Simon Willison's Weblog: geospatial

Coding agents for data analysis

2026-03-16T20:12:32+00:00

Here's the handout I prepared for my NICAR 2026 workshop "Coding agents for data analysis" - a three hour session aimed at data journalists demonstrating ways that tools like Claude Code and OpenAI Codex can be used to explore, analyze and clean data.

Here's the table of contents:

Coding agents

Warmup: ChatGPT and Claude

Setup Claude Code and Codex

Asking questions against a database

Exploring data with agents

Cleaning data: decoding neighborhood codes

Creating visualizations with agents

Scraping data with agents

I ran the workshop using GitHub Codespaces and OpenAI Codex, since it was easy (and inexpensive) to distribute a budget-restricted API key for Codex that attendees could use during the class. Participants ended up burning $23 of Codex tokens.

The exercises all used Python and SQLite and some of them used Datasette.

One highlight of the workshop was when we started running Datasette such that it served static content from a viz/ folder, then had Claude Code start vibe coding new interactive visualizations directly in that folder. Here's a heat map it created for my trees database using Leaflet and Leaflet.heat, source code here.

I designed the handout to also be useful for people who weren't able to attend the session in person. As is usually the case, material aimed at data journalists is equally applicable to anyone else with data to explore.

Tags: data-journalism, geospatial, python, speaking, sqlite, ai, datasette, generative-ai, llms, github-codespaces, nicar, coding-agents, claude-code, codex, leaflet

Let the LLM Write the Prompts: An Intro to DSPy in Compound Al Pipelines

2025-10-04T22:48:59+00:00

Let the LLM Write the Prompts: An Intro to DSPy in Compound Al Pipelines

I've had trouble getting my head around DSPy in the past. This half hour talk by Drew Breunig at the recent Databricks Data + AI Summit is the clearest explanation I've seen yet of the kinds of problems it can help solve.

Here's Drew's written version of the talk.

Drew works on Overture Maps, which combines Point Of Interest data from numerous providers to create a single unified POI database. This is an example of conflation, a notoriously difficult task in GIS where multiple datasets are deduped and merged together.

Drew uses an inexpensive local model, Qwen3-0.6B, to compare 70 million addresses and identity matches, for example between Place(address="3359 FOOTHILL BLVD", name="RESTAURANT LOS ARCOS") and Place(address="3359 FOOTHILL BLVD", name="Los Arcos Taqueria"').

DSPy's role is to optimize the prompt used for that smaller model. Drew used GPT-4.1 and the dspy.MIPROv2 optimizer, producing a 700 token prompt that increased the score from 60.7% to 82%.

Why bother? Drew points out that having a prompt optimization pipeline makes it trivial to evaluate and switch to other models if they can score higher with a custom optimized prompt - without needing to execute that trial-and-error optimization by hand.

Tags: geospatial, ai, prompt-engineering, generative-ai, llms, drew-breunig, overture, dspy

Spatial Joins in DuckDB

2025-08-23T21:21:02+00:00

Spatial Joins in DuckDB

Extremely detailed overview by Max Gabrielsson of DuckDB's new spatial join optimizations.

Consider the following query, which counts the number of NYC Citi Bike Trips for each of the neighborhoods defined by the NYC Neighborhood Tabulation Areas polygons and returns the top three:

SELECT neighborhood,
  count(*) AS num_rides
FROM rides
JOIN hoods ON ST_Intersects(
  rides.start_geom, hoods.geom
)
GROUP BY neighborhood
ORDER BY num_rides DESC
LIMIT 3;

The rides table contains 58,033,724 rows. The hoods table has polygons for 310 neighborhoods.

Without an optimized spatial joins this query requires a nested loop join, executing that expensive ST_Intersects() operation 58m * 310 ~= 18 billion times. This took around 30 minutes on the 36GB MacBook M3 Pro used for the benchmark.

The first optimization described - implemented from DuckDB 1.2.0 onwards - uses a "piecewise merge join". This takes advantage of the fact that a bounding box intersection is a whole lot faster to calculate, especially if you pre-cache the bounding box (aka the minimum bounding rectangle or MBR) in the stored binary GEOMETRY representation.

Rewriting the query to use a fast bounding box intersection and then only running the more expensive ST_Intersects() filters on those matches drops the runtime from 1800 seconds to 107 seconds.

The second optimization, added in DuckDB 1.3.0 in May 2025 using the new SPATIAL_JOIN operator, is significantly more sophisticated.

DuckDB can now identify when a spatial join is working against large volumes of data and automatically build an in-memory R-Tree of bounding boxes for the larger of the two tables being joined.

This new R-Tree further accelerates the bounding box intersection part of the join, and drops the runtime down to just 30 seconds.

Via @mackaszechno.bsky.social

Tags: geospatial, sql, duckdb

I counted all of the yurts in Mongolia using machine learning

2025-06-18T14:11:33+00:00

I counted all of the yurts in Mongolia using machine learning

Fascinating, detailed account by Monroe Clinton of a geospatial machine learning project. Monroe wanted to count visible yurts in Mongolia using Google Maps satellite view. The resulting project incorporates mercantile for tile calculations, Label Studio for help label the first 10,000 examples, a model trained on top of YOLO11 and a bunch of clever custom Python code to co-ordinate a brute force search across 120 CPU workers running the model.

Via Hacker News

Tags: geospatial, machine-learning, python, ai

DuckDB is Probably the Most Important Geospatial Software of the Last Decade

2025-05-04T00:28:35+00:00

DuckDB is Probably the Most Important Geospatial Software of the Last Decade

Drew Breunig argues that the ease of installation of DuckDB is opening up geospatial analysis to a whole new set of developers.

This inspired a comment on Hacker News from DuckDB Labs geospatial engineer Max Gabrielsson which helps explain why the drop in friction introduced by DuckDB is so significant:

I think a big part is that duckdbs spatial extension provides a SQL interface to a whole suite of standard foss gis packages by statically bundling everything (including inlining the default PROJ database of coordinate projection systems into the binary) and providing it for multiple platforms (including WASM). I.E there are no transitive dependencies except libc.

[...] the fact that you can e.g. convert too and from a myriad of different geospatial formats by utilizing GDAL, transforming through SQL, or pulling down the latest overture dump without having the whole workflow break just cause you updated QGIS has probably been the main killer feature for a lot of the early adopters.

I've lost count of the time I've spent fiddling with dependencies like GDAL trying to get various geospatial tools to work in the past. Bundling difficult dependencies statically is an under-appreciated trick!

If the bold claim in the headline inspires you to provide a counter-example, bear in mind that a decade ago is 2015, and most of the key technologies In the modern geospatial stack - QGIS, PostGIS, geopandas, SpatiaLite - predate that by quite a bit.

Tags: geospatial, sql, duckdb, drew-breunig

The GeoGuessr StreetView meta-game

2025-04-26T16:56:59+00:00

My post on o3 guessing locations from photos made it to Hacker News and by far the most interesting comments are from SamPatt, a self-described competitive GeoGuessr player.

In a thread about meta-knowledge of the StreetView card uses in different regions:

The photography matters a great deal - they're categorized into "Generations" of coverage. Gen 2 is low resolution, Gen 3 is pretty good but has a distinct car blur, Gen 4 is highest quality. Each country tends to have only one or two categories of coverage, and some are so distinct you can immediately know a location based solely on that (India is the best example here). [...]

Nigeria and Tunisia have follow cars. Senegal, Montenegro and Albania have large rifts in the sky where the panorama stitching software did a poor job. Some parts of Russia had recent forest fires and are very smokey. One road in Turkey is in absurdly thick fog. The list is endless, which is why it's so fun!

Sam also has his own custom Obsidian flashcard deck "with hundreds of entries to help me remember road lines, power poles, bollards, architecture, license plates, etc".

I asked Sam how closely the GeoGuessr community track updates to street view imagery, and unsurprisingly those are a big deal. Sam pointed me to this 10 minute video review by zi8gzag of the latest big update from three weeks ago:

This is one of the biggest updates in years in my opinion. It could be the biggest update since the 2022 update that gave Gen 4 to Nigeria, Senegal, and Rwanda. It's definitely on the same level as the Kazakhstan update or the Germany update in my opinion.

Tags: geospatial, hacker-news, streetview, geoguessing

A5

2025-04-22T14:37:36+00:00

A5 is a new "global, equal-area, millimeter-accurate geospatial index" by Felix Palmer:

It is the pentagonal equivalent of other DGGSs, like S2 or H3, but with higher accuracy and lower distortion.

Effectively it's a way of dividing the entire world into pentagons where each one covers the same physical area (to within a 2% threshold) - like Uber's H3 but a bit weirder and more fun. An A5 reference implementation written in TypeScript is available on GitHub.

This interactive demo helps show how it works:

Why pentagons? Here's what the A5 docs say:

A5 is unique in that it uses a pentagonal tiling of a dodecahedron. [...] The benefit of choosing a dodecahedron is that it is the platonic solid with the lowest vertex curvature, and by this measure it is the most spherical of all the platonic solids. This is key for minimizing cell distortion as the process of projecting a platonic solid onto a sphere involves warping the cell geometry to force the vertex curvature to approach zero. Thus, the lower the original vertex curvature, the less distortion will be introduced by the projection.

I had to look up platonic solids on Wikipedia. There are only five: Tetrahedron, Cube, Octahedron, Dodecahedron and Icosahedron and they can be made using squares, triangles or (in the case of the Dodecahedron) pentagons, making the pentagon the most circle-like option.

Via Tom MacWright

Tags: geospatial

OpenTimes

2025-03-17T22:49:59+00:00

OpenTimes

Spectacular new open geospatial project by Dan Snow:

OpenTimes is a database of pre-computed, point-to-point travel times between United States Census geographies. It lets you download bulk travel time data for free and with no limits.

Here's what I get for travel times by car from El Granada, California:

The technical details are fascinating:

The entire OpenTimes backend is just static Parquet files on Cloudflare's R2. There's no RDBMS or running service, just files and a CDN. The whole thing costs about $10/month to host and costs nothing to serve. In my opinion, this is a great way to serve infrequently updated, large public datasets at low cost (as long as you partition the files correctly).

Sure enough, R2 pricing charges "based on the total volume of data stored" - $0.015 / GB-month for standard storage, then $0.36 / million requests for "Class B" operations which include reads. They charge nothing for outbound bandwidth.

All travel times were calculated by pre-building the inputs (OSM, OSRM networks) and then distributing the compute over hundreds of GitHub Actions jobs. This worked shockingly well for this specific workload (and was also completely free).

Here's a GitHub Actions run of the calculate-times.yaml workflow which uses a matrix to run 255 jobs!

Relevant YAML:

  matrix:
    year: ${{ fromJSON(needs.setup-jobs.outputs.years) }}
    state: ${{ fromJSON(needs.setup-jobs.outputs.states) }}

Where those JSON files were created by the previous step, which reads in the year and state values from this params.yaml file.

The query layer uses a single DuckDB database file with views that point to static Parquet files via HTTP. This lets you query a table with hundreds of billions of records after downloading just the ~5MB pointer file.

This is a really creative use of DuckDB's feature that lets you run queries against large data from a laptop using HTTP range queries to avoid downloading the whole thing.

The README shows how to use that from R and Python - I got this working in the duckdb client (brew install duckdb):

INSTALL httpfs;
LOAD httpfs;
ATTACH 'https://data.opentimes.org/databases/0.0.1.duckdb' AS opentimes;

SELECT origin_id, destination_id, duration_sec
  FROM opentimes.public.times
  WHERE version = '0.0.1'
      AND mode = 'car'
      AND year = '2024'
      AND geography = 'tract'
      AND state = '17'
      AND origin_id LIKE '17031%' limit 10;

In answer to a question about adding public transit times Dan said:

In the next year or so maybe. The biggest obstacles to adding public transit are:

Collecting all the necessary scheduling data (e.g. GTFS feeds) for every transit system in the county. Not insurmountable since there are services that do this currently.

Finding a routing engine that can compute nation-scale travel time matrices quickly. Currently, the two fastest open-source engines I've tried (OSRM and Valhalla) don't support public transit for matrix calculations and the engines that do support public transit (R5, OpenTripPlanner, etc.) are too slow.

GTFS is a popular CSV-based format for sharing transit schedules - here's an official list of available feed directories.

This whole project feels to me like a great example of the baked data architectural pattern in action.

Via Hacker News

Tags: census, geospatial, open-data, openstreetmap, cloudflare, parquet, github-actions, baked-data, duckdb, http-range-requests

wolf-h3-viewer.glitch.me

2025-03-09T14:51:55+00:00

wolf-h3-viewer.glitch.me

Neat interactive visualization of Uber's H3 hexagonal geographical indexing mechanism.

Here's the source code.

Why does H3 use hexagons? Because Hexagons are the Bestagons:

When hexagons come together, they form three-sided joints 120 degrees apart. This, for the least material, is the most mechanically stable arrangement.

Only triangles, squares, and hexagons can tile a plane without gaps, and of those three shapes hexagons offer the best ratio of perimeter to area.

Via Hacker News comment

Tags: geospatial, javascript, leaflet, uber

OpenStreetMap embed URL

2024-11-25T19:29:16+00:00

OpenStreetMap embed URL

I just found out OpenStreetMap have a "share" button which produces HTML for an iframe targetting https://www.openstreetmap.org/export/embed.html, making it easy to drop an OpenStreetMap map onto any web page that allows iframes.

As far as I can tell the supported parameters are:

bbox= then min longitude, min latitude, max longitude, max latitude
marker= optional latitude, longitude coordinate for a marker (only a single marker is supported)
layer=mapnik - other values I've found that work are cyclosm, cyclemap, transportmap and hot (for humanitarian)

Here's HTML for embedding this on a page using a sandboxed iframe - the allow-scripts is necessary for the map to display.

<iframe
  sandbox="allow-scripts"
  style="border: none; width: 100%; height: 20em;"
  src="https://www.openstreetmap.org/export/embed.html?bbox=-122.613%2C37.431%2C-122.382%2C37.559&amp;layer=mapnik&amp;marker=37.495%2C-122.497"
></iframe>

Thanks to this post I learned that iframes are rendered correctly in NetNewsWire, NewsExplorer, NewsBlur and Feedly on Android.

Tags: geospatial, iframes, openstreetmap, sandboxing

Foursquare Open Source Places: A new foundational dataset for the geospatial community

2024-11-20T05:52:38+00:00

Foursquare Open Source Places: A new foundational dataset for the geospatial community

I did not expect this!

[...] we are announcing today the general availability of a foundational open data set, Foursquare Open Source Places ("FSQ OS Places"). This base layer of 100mm+ global places of interest ("POI") includes 22 core attributes (see schema here) that will be updated monthly and available for commercial use under the Apache 2.0 license framework.

The data is available as Parquet files hosted on Amazon S3.

Here's how to list the available files:

aws s3 ls s3://fsq-os-places-us-east-1/release/dt=2024-11-19/places/parquet/

I got back places-00000.snappy.parquet through places-00024.snappy.parquet, each file around 455MB for a total of 10.6GB of data.

I ran duckdb and then used DuckDB's ability to remotely query Parquet on S3 to explore the data a bit more without downloading it to my laptop first:

select count(*) from 's3://fsq-os-places-us-east-1/release/dt=2024-11-19/places/parquet/places-00000.snappy.parquet';

This got back 4,180,424 - that number is similar for each file, suggesting around 104,000,000 records total.

Update: DuckDB can use wildcards in S3 paths (thanks, Paul) so this query provides an exact count:

select count(*) from 's3://fsq-os-places-us-east-1/release/dt=2024-11-19/places/parquet/places-*.snappy.parquet';

That returned 104,511,073 - and Activity Monitor on my Mac confirmed that DuckDB only needed to fetch 1.2MB of data to answer that query.

I ran this query to retrieve 1,000 places from that first file as newline-delimited JSON:

copy (
    select * from 's3://fsq-os-places-us-east-1/release/dt=2024-11-19/places/parquet/places-00000.snappy.parquet'
    limit 1000
) to '/tmp/places.json';

Here's that places.json file, and here it is imported into Datasette Lite.

Finally, I got ChatGPT Code Interpreter to convert that file to GeoJSON and pasted the result into this Gist, giving me a map of those thousand places (because Gists automatically render GeoJSON):

Via Andy Baio

Tags: geospatial, open-source, foursquare, geojson, parquet, duckdb, datasette-lite, ai-assisted-programming, code-interpreter, coding-agents

OpenStreetMap vector tiles demo

2024-11-19T23:39:18+00:00

OpenStreetMap vector tiles demo

Long-time OpenStreetMap developer Paul Norman has been working on adding vector tile support to OpenStreetMap for quite a while. Paul recently announced that vector.openstreetmap.org is now serving vector tiles (in Mapbox Vector Tiles (MVT) format) - here's his interactive demo for seeing what they look like.

Via Mark Litwintschik

Tags: geospatial, mapping, openstreetmap

Visualizing local election results with Datasette, Observable and MapLibre GL

2024-11-09T23:32:06+00:00

Alex Garcia and myself hosted the first Datasette Open Office Hours on Friday - a live-streamed video session where we hacked on a project together and took questions and tips from community members on Discord.

We didn't record this one (surprisingly not a feature that Discord offers) but we hope to do more of these and record them in the future.

This post is a detailed write-up of what we built during the session.

San Mateo County election results

I live in El Granada, a tiny town just north of Half Moon Bay in San Mateo County, California.

Every county appears to handle counting and publishing election results differently. For San Mateo County the results are published on this page, and detailed per-precinct and per-candidate breakdowns are made available as a CSV file.

(I optimistically set up a Git scraper for these results in simonw/scrape-san-mateo-county-election-results-2024 only to learn that the CSV is updated just once a day, not continually as the ballots are counted.)

I'm particularly invested in the results of the Granada Community Services District board member elections. Our little town of El Granada is in "unincorporated San Mateo County" which means we don't have a mayor or any local officials, so the closest we get to hyper-local government is the officials that run our local sewage and parks organization! My partner Natalie ran the candidate forum event (effectively the debate) featuring three of the four candidates running for the two open places on the board.

Let's explore the data for that race using Datasette.

Importing CSV data into Datasette

I ran my part of the demo using Datasette Cloud, the beta of my new hosted Datasette service.

I started by using the pre-configured datasette-import plugin to import the data from the CSV file into a fresh table:

Modifying the schema

The table imported cleanly, but all of the columns from the CSV were still being treated as text. I used the datasette-edit-schema plugin to switch the relevant columns to integers so that we could run sums and sorts against them.

(I also noted that I really should add a "detect column types" feature to that plugin!)

The resulting 15,589 rows represent counts from individual precincts around the county for each of the races and measures on the ballot, with a row per precinct per candidate/choice per race.

Faceting and filtering the table

Since I'm interested in the Granada Community Services District election, I applied a facet on "Contest_title" and then used that to select that specific race.

I applied additional facets on "candidate_name" and "Precinct name".

This looks right to me: we have 7 precincts and 4 candidates for 28 rows in total.

Importing geospatial precinct shapes

Those precinct names are pretty non-descriptive! What does 33001 mean?

To answer that question, I added a new table.

San Mateo County offers precinct maps in the form of 23 PDF files. Our precincts are in the "Unincorporated Coastside" file:

Thankfully the county also makes that data available as geospatial data, hosted using Socrata with an option to export as GeoJSON.

The datasette-import plugin can handle JSON files... and if a JSON file contains a top-level object with a key that is an array of objects, it will import those objects as a table.

Dragging that file into Datasette is enough to import it as a table with a properties JSON column containing properties and a geometry JSON columnn with the GeoJSON geometry.

Here's where another plugin kicks in: datasette-leaflet-geojson looks for columns that contain valid GeoJSON geometries and... draws them on a map!

So now we can see the shape of the individual geometries.

Enriching that data to extract the precinct IDs

The precinctid is present in the data, but it's tucked away in a JSON object in that properties JSON blob. It would be more convenient if it was a top-level column.

Datasette's enrichments feature provides tools for running operations against every row in a table and adding new columns based on the results.

My Datasette Cloud instance was missing the datasette-enrichments-quickjs plugin that would let me run JavaScript code against the data. I used my privileged access on Datasette Cloud to add that plugin to my requirements and restarted the instance to install it.

I used that to run this JavaScript code against every row in the table and saved the output in a new precinct_id column:

function enrich(row) {
  return JSON.parse(row.properties).precinctid;
}

This took less than a second to run, adding and populating a new precinct_id column for the table.

Running a join

I demonstrated how to run a join between the election results and the precincts table using the Datasette SQL query editor.

I tried a few different things, but the most interesting query was this one:

select
  Precinct_name,
  precincts.geometry,
  total_ballots,
  json_group_object(
    candidate_name,
    total_votes
  ) as votes_by_candidate
from
  election_results 
  join precincts on election_results.Precinct_name = precincts.precinct_id
where 
  Contest_title = "Granada Community Services District Members, Board of Directors"
group by 
  Precinct_name,
  precincts.geometry,
  total_ballots;

Creating an API token to access the data

I was nearly ready to hand over to Alex for the second half of our demo, where he would use Observable Notebooks to build some custom visualizations on top of the data.

A great pattern for this is to host the data in Datasette and then fetch it into Observable via the Datasette JSON API.

Since Datasette Cloud instances are private by default we would need to create an API token that could do this.

I used this interface (from the datasette-auth-tokens plugin) to create a new token with read-only access to all databases and tables in the instance:

Since we're running a dedicated instance just for Datasette Public Office Hours there's no reason not to distribute that read-only token in publically accessible code.

Getting CORS working

Embarrassingly, I had forgotten that we would need CORS headers in order to access the data from an Observable notebook. Thankfully we have another plugin for that: datasette-cors. I installed that quickly and we confirmed that it granted access to the API from Observable as intended.

I handed over to Alex for the next section of the demo.

Working with Datasette in Observable

Alex started by running a SQL query from client-side JavaScript to pull in the joined data for our specific El Granada race:

sql = `
select
  Precinct_name,
  precincts.geometry,
  Split_name,
  Reporting_flag,
  Update_count,
  Pct_Id,
  Pct_seq_nbr,
  Reg_voters,
  Turn_Out,
  Contest_Id,
  Contest_seq_nbr,
  Contest_title,
  Contest_party_name,
  Selectable_Options,
  candidate_id,
  candidate_name,
  Candidate_Type,
  cand_seq_nbr,
  Party_Code,
  total_ballots,
  total_votes,
  total_under_votes,
  total_over_votes,
  [Vote Centers_ballots],
  [Vote Centers_votes],
  [Vote Centers_under_votes],
  [Vote Centers_over_votes],
  [Vote by Mail_ballots],
  [Vote by Mail_votes],
  [Vote by Mail_under_votes],
  [Vote by Mail_over_votes]
from
  election_results join precincts on election_results.Precinct_name = precincts.precinct_id
where "Contest_title" = "Granada Community Services District Members, Board of Directors"
limit 101;`

And in the next cell:

raw_data = fetch(
  `https://datasette-public-office-hours.datasette.cloud/data/-/query.json?_shape=array&sql=${encodeURIComponent(
    sql
  )}`,
  {
    headers: {
      Authorization: `Bearer ${secret}`
    }
  }
).then((r) => r.json())

Note the ?_shape=array parameter there, which causes Datasette to output the results directly as a JSON array of objects.

That's all it takes to get the data into Observable. Adding another cell like this confirms that the data is now available:

Inputs.table(raw_data)

Visualizing those with MapLibre GL

There are plenty of good options for visualizing GeoJSON data using JavaScript in an Observable notebook.

Alex started with MapLibre GL, using the excellent OpenFreeMap 3D tiles:

viewof map = {
  const container = html`<div style="height:800px;">`;
  yield container;
  const map = (container.value = new maplibregl.Map({
    container,
    zoom: 2,
    //style: "https://basemaps.cartocdn.com/gl/voyager-gl-style/style.json",
    style: "https://tiles.openfreemap.org/styles/liberty",
    scrollZoom: true
  }));
  yield container;

  map.on("load", function () {
    map.fitBounds(d3.geoBounds(data), { duration: 0 });
    map.addSource("precincts", {
      type: "geojson",
      data: data
    });
    map.addLayer({
      id: "precincts",
      type: "fill",
      source: "precincts",
      paint: {
        "fill-opacity": 0.4,
        "fill-color": [
          "case",
          ["==", ["get", "ratio"], null], "#000000",
          [
            "interpolate",
            ["linear"],
            ["get", "ratio"],
            0.0, "#0000ff",
            0.5, "#d3d3d3",
            1.0, "#ff0000"
          ]
        ]
      }
    });
    map.on("click", "precincts", (e) => {
      const { precinct, ratio } = e.features[0].properties;
      const description = JSON.stringify();

      new maplibregl.Popup()
        .setLngLat(e.lngLat)
        .setHTML(description)
        .addTo(map);
    });
  });
  invalidation.then(() => map.remove());
}

(This is just one of several iterations, I didn't capture detailed notes of every change Alex made to the code.)

Observable Plot

Observable notebooks come pre-loaded with the excellent Observable Plot charting library - Mike Bostock's high-level charting tool built on top of D3.

Alex used that to first render the shapes of the precincts directly, without even needing a tiled basemap:

Plot.plot({
  width,
  height: 600,
  legend: true,
  projection: {
    type: "conic-conformal",
    parallels: [37 + 4 / 60, 38 + 26 / 60],
    rotate: [120 + 30 / 60, 0],
    domain: data
  },
  marks: [
    Plot.geo(data, {
      strokeOpacity: 0.1,
      fill: "total_votes",
      title: (d) => JSON.stringify(d.properties),
      tip: true
    })
  ]
})

The parallels and rotate options there come from the handy veltman/d3-stateplane repo, which lists recommended settings for the State Plane Coordinate System used with projections in D3. Those values are for California Zone 3.

Bringing it all together

For the grand finale, Alex combined everything learned so far to build an interactive map allowing a user to select any of the 110 races on the ballot and see a heatmap of results for any selected candidate and option:

You can try this out in Alex's notebook. Here's the relevant code (Observable cells are divided by // --- comments). Note that Observable notebooks are reactive and allow variables to be referenced out of order.

// Select the contest
viewof contest = Inputs.select(contests, { label: "Choose a contest" })

// ---

// And the candidate
viewof candidate = Inputs.radio(
  candidates,
  {
    label: "Choose a candidate",
    value: candidates[0]
  }
)

// ---

// Show the map itself
Plot.plot({
  width,
  height: 600,
  legend: true,
  color: { scheme: "blues", legend: true },
  projection: {
    type: "mercator",
    domain: data2
  },
  marks: [
    Plot.geo(data2, {
      strokeOpacity: 0.1,
      fill: "ratio",
      tip: true
    })
  ]
})
# ---
data2 = ({
  type: "FeatureCollection",
  features: raw_data2.map((d) => ({
    type: "Feature",
    properties: {
      precinct: d.Precinct_name,
      total_ballots: d.total_ballots,
      ratio: JSON.parse(d.votes_by_candidate)[candidate] / d.total_ballots
    },
    geometry: JSON.parse(d.geometry)
  }))
})

// ---

raw_data2 = query(
  `select
  Precinct_name,
  precincts.geometry,
  total_ballots,
  json_grop_object(
    candidate_name,
    total_votes
  ) as votes_by_candidate
from
  election_results 
  join precincts on election_results.Precinct_name = precincts.precinct_id
where Contest_title = :contest
group by 
  Precinct_name,
  precincts.geometry,
  total_ballots;`,
  { contest }
)

// ---

raw_data2 = query(
  `select
  Precinct_name,
  precincts.geometry,
  total_ballots,
  json_group_object(
    candidate_name,
    total_votes
  ) as votes_by_candidate
from
  election_results 
  join precincts on election_results.Precinct_name = precincts.precinct_id
where Contest_title = :contest
group by 
  Precinct_name,
  precincts.geometry,
  total_ballots;`,
  { contest }
)

// ---

// Fetch the available contests
contests = query("select distinct Contest_title from election_results").then(
  (d) => d.map((d) => d.Contest_title)
)

// ---

// Extract available candidates for selected contest

candidates = Object.keys(
  JSON.parse(raw_data2[0].votes_by_candidate)
)

// ---

function query(sql, params = {}) {
  return fetch(
    `https://datasette-public-office-hours.datasette.cloud/data/-/query.json?${new URLSearchParams(
      { sql, _shape: "array", ...params }
    ).toString()}`,
    {
      headers: {
        Authorization: `Bearer ${secret}`
      }
    }
  ).then((r) => r.json());
}

We'll be doing this again

This was our first time trying something like this and I think it worked really well. We're already thinking about ways to improve it next time:

I want to record these sessions and make them available on YouTube for people who couldn't be there live
It would be fun to mix up the format. I'm particularly keen on getting more people involved giving demos - maybe having 5-10 minute lightning demo slots so we can see what other people are working on

Keep an eye on this blog or on the Datasette Discord for news about future sessions.

Tags: geospatial, mapping, politics, projects, datasette, datasette-cloud, alex-garcia, datasette-public-office-hours, leaflet

Conflating Overture Places Using DuckDB, Ollama, Embeddings, and More

2024-09-30T17:24:03+00:00

Conflating Overture Places Using DuckDB, Ollama, Embeddings, and More

Drew Breunig's detailed tutorial on "conflation" - combining different geospatial data sources by de-duplicating address strings such as RESTAURANT LOS ARCOS,3359 FOOTHILL BLVD,OAKLAND,94601 and LOS ARCOS TAQUERIA,3359 FOOTHILL BLVD,OAKLAND,94601.

Drew uses an entirely offline stack based around Python, DuckDB and Ollama and finds that a combination of H3 geospatial tiles and mxbai-embed-large embeddings (though other embedding models should work equally well) gets really good results.

Tags: geospatial, python, ai, duckdb, embeddings, drew-breunig, overture, ollama

OpenFreeMap

2024-09-28T21:41:15+00:00

OpenFreeMap

New free map tile hosting service from Zsolt Ero:

OpenFreeMap lets you display custom maps on your website and apps for free. […] Using our public instance is completely free: there are no limits on the number of map views or requests. There’s no registration, no user database, no API keys, and no cookies. We aim to cover the running costs of our public instance through donations.

The site serves static vector tiles that work with MapLibre GL. It deliberately doesn’t offer any other services such as search or routing.

From the project README looks like it’s hosted on two Hetzner machines. I don’t think the public server is behind a CDN.

Part of the trick to serving the tiles efficiently is the way it takes advantage of Btrfs:

Production-quality hosting of 300 million tiny files is hard. The average file size is just 450 byte. Dozens of tile servers have been written to tackle this problem, but they all have their limitations.

The original idea of this project is to avoid using tile servers altogether. Instead, the tiles are directly served from Btrfs partition images + hard links using an optimised nginx config.

The self-hosting guide describes the scripts that are provided for downloading their pre-built tiles (needing a fresh Ubuntu server with 300GB of SSD and 4GB of RAM) or building the tiles yourself using Planetiler (needs 500GB of disk and 64GB of RAM).

Getting started is delightfully straightforward:

const map = new maplibregl.Map({
  style: 'https://tiles.openfreemap.org/styles/liberty',
  center: [13.388, 52.517],
  zoom: 9.5,
  container: 'map',
})

I got Claude to help build this demo showing a thousand random markers dotted around San Francisco. The 3D tiles even include building shapes!

Zsolt built OpenFreeMap based on his experience running MapHub over the last 9 years. Here’s a 2018 interview about that project.

It’s pretty incredible that the OpenStreetMap and open geospatial stack has evolved to the point now where it’s economically feasible for an individual to offer a service like this. I hope this turns out to be sustainable. Hetzner charge just €1 per TB for bandwidth (S3 can cost $90/TB) which should help a lot.

Via Tim Bray

Tags: geospatial, maps, openstreetmap, ai-assisted-programming

Towards Standardizing Place

2024-08-01T23:14:49+00:00

Towards Standardizing Place

Overture Maps announced General Availability of its global maps datasets last week, covering places, buildings, divisions, and base layers.

Drew Breunig demonstrates how this can be accessed using both the Overture Explorer tool and DuckDB, and talks about Overture's GERS IDs - reminiscent of Who's On First IDs - which provide stable IDs for all kinds of geographic places.

Tags: geospatial, whosonfirst, drew-breunig, overture

The many lives of Null Island

2024-07-28T17:44:54+00:00

The many lives of Null Island

Stamen's custom basemaps have long harbored an Easter egg: zoom all the way in on 0, 0 to see the outline of the mystical "null island", the place where GIS glitches and data bugs accumulate, in the Gulf of Guinea south of Ghana.

Stamen's Alan McConchie provides a detailed history of the Easter egg - first introduced by Mike Migurski in 2010 - along with a definitive guide to the GIS jokes and traditions that surround it.

Here's Null Island on Stamen's Toner map. The shape (also available as GeoJSON) is an homage to the island from 1993's Myst, hence the outline of a large docked ship at the bottom.

Alan recently gave a talk about Stamen's updated custom maps at State of the Map US 2024 (video, slides) - their Toner and Terrain maps are now available as vector tiles served by Stadia Maps (here's the announcement), but their iconic watercolor style is yet to be updated to vectors, due to the weird array of raster tricks it used to achieve the effect.

In researching this post I searched for null island on Google Maps and was delighted to learn that a bunch of entrepreneurs in Western Africa have tapped into the meme for their own businesses:

Via Andy Baio

Tags: geospatial, maps, michal-migurski, stamen-design

Searching an aerial photo with text queries

2024-07-12T18:07:48+00:00

Searching an aerial photo with text queries

Robin Wilson built a demo that lets you search a large aerial photograph of Southampton for things like "roundabout" or "tennis court". He explains how it works in detail: he used the SkyCLIP model, which is trained on "5.2 million remote sensing image-text pairs in total, covering more than 29K distinct semantic tags" to generate embeddings for 200x200 image segments (with 100px of overlap), then stored them in Pinecone.

Tags: geospatial, embeddings, clip

tiny-world-map

2024-04-21T22:11:14+00:00

tiny-world-map

I love this project. It’s a JavaScript file (694K uncompressed, 283KB compressed) which can be used with the Leaflet mapping library and provides a SVG base map of the world with country borders and labels for every world city with a population more than 48,000—10,000 cities total.

This means you can bundle an offline map of the world as part of any application that doesn’t need a higher level of detail. A lot of smaller island nations are missing entirely though, so this may not be right for every project.

It even includes a service worker to help implement offline mapping support, plus several variants of the map with less cities that are even smaller.

Via Hacker News

Tags: geospatial, javascript, mapping, svg, service-workers, leaflet

A POI Database in One Line

2024-04-19T02:44:58+00:00

A POI Database in One Line

Overture maps offer an extraordinarily useful freely licensed databases of POI (point of interest) listings, principally derived from partners such as Facebook and including restaurants, shops, museums and other locations from all around the world.

Their new "overturemaps" Python CLI utility makes it easy to quickly pull subsets of their data... but requires you to provide a bounding box to do so.

Drew Breunig came up with this delightful recipe for fetching data using LLM and gpt-3.5-turbo to fill in those bounding boxes:

overturemaps download --bbox=$(llm 'Give me a bounding box for Alameda, California expressed as only four numbers delineated by commas, with no spaces, longitude preceding latitude.') -f geojsonseq --type=place | geojson-to-sqlite alameda.db places - --nl --pk=id

Via @dbreunig

Tags: cli, geospatial, geojson, llm, drew-breunig, overture

mapshaper.org

2024-03-23T03:44:22+00:00

mapshaper.org

It turns out the mapshaper CLI tool for manipulating geospatial data—including converting shapefiles to GeoJSON and back again—also has a web UI that runs the conversions entirely in your browser. If you need to convert between those (and other) formats it’s hard to imagine a more convenient option.

Via @adamrpearce

Tags: cli, geospatial, javascript, shapefiles, geojson

Claude and ChatGPT for ad-hoc sidequests

2024-03-22T19:44:12+00:00

Here is a short, illustrative example of one of the ways in which I use Claude and ChatGPT on a daily basis.

I recently learned that the Adirondack Park is the single largest park in the contiguous United States, taking up a fifth of the state of New York.

Naturally, my first thought was that it would be neat to have a GeoJSON file representing the boundary of the park.

A quick search landed me on the Adirondack Park Agency GIS data page, which offered me a shapefile of the "Outer boundary of the New York State Adirondack Park as described in Section 9-0101 of the New York Environmental Conservation Law". Sounds good!

I knew there were tools for converting shapefiles to GeoJSON, but I couldn't remember what they were. Since I had a terminal window open already, I typed the following:

llm -m opus -c 'give me options on macOS for CLI tools to turn a shapefile into GeoJSON'

Here I am using my LLM tool (and llm-claude-3 plugin) to run a prompt through the new Claude 3 Opus, my current favorite language model.

It replied with a couple of options, but the first was this:

ogr2ogr -f GeoJSON output.geojson input.shp

So I ran that against the shapefile, and then pasted the resulting GeoJSON into geojson.io to check if it worked... and nothing displayed. Then I looked at the GeoJSON and spotted this:

"coordinates": [ [ -8358911.527799999341369, 5379193.197800002992153 ] ...

That didn't look right. Those co-ordinates aren't the correct scale for latitude and longitude values.

So I sent a follow-up prompt to the model (the -c option means "continue previous conversation"):

llm -c 'i tried using ogr2ogr but it gave me back GeoJSON with a weird coordinate system that was not lat/lon that i am used to'

It suggested this new command:

ogr2ogr -f GeoJSON -t_srs EPSG:4326 output.geojson input.shp

This time it worked! The shapefile has now been converted to GeoJSON.

Time elapsed so far: 2.5 minutes (I can tell from my LLM logs).

I pasted it into Datasette (with datasette-paste and datasette-leaflet-geojson) to take a look at it more closely, and got this:

That's not a single polygon! That's 106 line segments... and they are fascinating. Look at those descriptions:

thence westerly along the northern line of lots 204 and 203 to the midpoint of the northern line of lot 203

This is utterly delightful. The shapefile description did say "as described in Section 9-0101 of the New York Environmental Conservation Law", so I guess this is how you write geographically boundaries into law!

But it's not what I wanted. I want a single polygon of the whole park, not 106 separate lines.

I decided to switch models. ChatGPT has access to Code Interpreter, and I happen to know that Code Interpreter is quite effective at processing GeoJSON.

I opened a new ChatGPT (with GPT-4) browser tab, uploaded my GeoJSON file and prompted it:

This GeoJSON file is full of line segments. Use them to create me a single shape that is a Polygon

OK, so it wrote some Python code and ran it. But did it work?

I happen to know that Code Interpreter can save files to disk and provide links to download them, so I told it to do that:

Save it to a GeoJSON file for me to download

I pasted that into geojson.io, and it was clearly wrong:

So I told it to try again. I didn't think very hard about this prompt, I basically went with a version of "do better":

that doesn't look right to me, check that it has all of the lines in it

It gave me a new file, optimistically named complete_polygon.geojson. Here's what that one looked like:

This is getting a lot closer! Note how the right hand boundary of the park looks correct, but the rest of the image is scrambled.

I had a hunch about the fix. I pasted in a screenshot of where we were so far and added my hunch about the solution:

That almost works but you need to sort the line segments first, it looked like this:

Honestly, pasting in the screenshot probably wasn't necessary here, but it amused me.

... and ChatGPT churned away again ...

sorted_polygon.geojson is spot on! Here's what it looks like:

Total time spent in ChatGPT: 3 minutes and 35 seconds. Plus 2.5 minutes with Claude 3 earlier, so an overall total of just over 6 minutes.

Here's the full Claude transcript and the full transcript from ChatGPT.

This isn't notable

The most notable thing about this example is how completely not notable it is.

I get results like this from these tools several times a day. I'm not at all surprised that this worked, in fact, I would've been mildly surprised if it had not.

Could I have done this without LLM assistance? Yes, but not nearly as quickly. And this was not a task on my critical path for the day - it was a sidequest at best and honestly more of a distraction.

So, without LLM tools, I would likely have given this one up at the first hurdle.

A year ago I wrote about how AI-enhanced development makes me more ambitious with my projects. They are now so firmly baked into my daily work that they influence not just side projects but tiny sidequests like this one as well.

This certainly wasn't simple

Something else I like about this example is that it illustrates quite how much depth there is to getting great results out of these systems.

In those few minutes I used two different interfaces to call two different models. I sent multiple follow-up prompts. I triggered Code Interpreter, took advantage of GPT-4 Vision and mixed in external tools like geojson.io and Datasette as well.

I leaned a lot on my existing knowledge and experience:

I knew that tools existed for commandline processing of shapefiles and GeoJSON
I instinctively knew that Claude 3 Opus was likely to correctly answer my initial prompt
I knew the capabilities of Code Interpreter, including that it has libraries that can process geometries, what to say to get it to kick into action and how to get it to give me files to download
My limited GIS knowledge was strong enough to spot a likely coordinate system problem, and I guessed the fix for the jumbled lines
My prompting intuition is developed to the point that I didn't have to think very hard about what to say to get the best results

If you have the right combination of domain knowledge and hard-won experience driving LLMs, you can fly with these things.

Isn't this a bit trivial?

Yes it is, and that's the point. This was a five minute sidequest. Writing about it here took ten times longer than the exercise itself.

I take on LLM-assisted sidequests like this one dozens of times a week. Many of them are substantially larger and more useful. They are having a very material impact on my work: I can get more done and solve much more interesting problems, because I'm not wasting valuable cycles figuring out ogr2ogr invocations or mucking around with polygon libraries.

Not to mention that I find working this way fun! It feels like science fiction every time I do it. Our AI-assisted future is here right now and I'm still finding it weird, fascinating and deeply entertaining.

LLMs are useful

There are many legitimate criticisms of LLMs. The copyright issues involved in their training, their enormous power consumption and the risks of people trusting them when they shouldn't (considering both accuracy and bias) are three that I think about a lot.

The one criticism I wont accept is that they aren't useful.

One of the greatest misconceptions concerning LLMs is the idea that they are easy to use. They really aren't: getting great results out of them requires a great deal of experience and hard-fought intuition, combined with deep domain knowledge of the problem you are applying them to.

I use these things every day. They help me take on much more interesting and ambitious problems than I could otherwise. I would miss them terribly if they were no longer available to me.

Tags: geospatial, shapefiles, geojson, ai, openai, generative-ai, chatgpt, llms, ai-assisted-programming, anthropic, claude, code-interpreter, coding-agents, leaflet

How to make self-hosted maps that work everywhere and cost next to nothing

2024-02-24T04:19:05+00:00

How to make self-hosted maps that work everywhere and cost next to nothing

Chris Amico provides a detailed roundup of the state of web mapping in 2024. It’s never been easier to entirely host your own mapping infrastructure, thanks to OpenStreetMap, Overture, MBTiles, PMTiles, Maplibre and a whole ecosystem of other fine open source projects.

I like Protomaps creator Brandon Liu’s description of this: “post-scarcity web mapping”.

Tags: geospatial, maps, openstreetmap, chris-amico, overture

Database generated columns: GeoDjango & PostGIS

2023-12-11T19:14:28+00:00

Database generated columns: GeoDjango & PostGIS

Paolo Melchiorre advocated for the inclusion of generated columns, one of the biggest features in Django 5.0. Here he provides a detailed tutorial showing how they can be used with PostGIS to create database tables that offer columns such as geohash that are automatically calculated from other columns in the table.

Tags: django, geospatial, postgresql

Geospatial SQL queries in SQLite using TG, sqlite-tg and datasette-sqlite-tg

2023-09-25T19:45:03+00:00

Geospatial SQL queries in SQLite using TG, sqlite-tg and datasette-sqlite-tg

Alex Garcia built sqlite-tg—a SQLite extension that uses the brand new TG geospatial library to provide a whole suite of custom SQL functions for working with geospatial data.

Here are my notes on trying out his initial alpha releases. The extension already provides tools for converting between GeoJSON, WKT and WKB, plus the all important tg_intersects() function for testing if a polygon or point overlap each other.

It’s pretty useful already. Without any geospatial indexing at all I was still able to get 700ms replies to a brute-force point-in-polygon query against 150MB of GeoJSON timezone boundaries stored as JSON text in a table.

Tags: geospatial, sqlite, geojson, datasette, alex-garcia, tg

TG: Polygon indexing

2023-09-23T04:32:14+00:00

TG: Polygon indexing

TG is a brand new geospatial library by Josh Baker, author of the Tile38 in-memory spatial server (kind of a geospatial Redis). TG is written in pure C and delivered as a single C file, reminiscent of the SQLite amalgamation.

TG looks really interesting. It implements almost the exact subset of geospatial functionality that I find most useful: point-in-polygon, intersect, WKT, WKB, and GeoJSON—all with no additional dependencies.

The most interesting thing about it is the way it handles indexing. In this documentation Josh describes two approaches he uses to speeding up point-in-polygon and intersection using a novel approach that goes beyond the usual RTree implementation.

I think this could make the basis of a really useful SQLite extension—a lighter-weight alternative to SpatiaLite.

Via Hacker News

Tags: c, geospatial, spatialite, sqlite, geojson, tg

Overture Maps Foundation Releases Its First World-Wide Open Map Dataset

2023-07-27T16:45:09+00:00

Overture Maps Foundation Releases Its First World-Wide Open Map Dataset

The Overture Maps Foundation is a collaboration lead by Amazon, Meta, Microsoft and TomTom dedicated to producing “reliable, easy-to-use, and interoperable open map data”.

Yesterday they put out their first release and it’s pretty astonishing: four different layers of geodata, covering Places of Interest (shops, restaurants, attractions etc), administrative boundaries, building outlines and transportation networks.

The data is available as Parquet. I just downloaded the 8GB places dataset and can confirm that it contains 59 million listings from around the world—I filtered to just places in my local town and a spot check showed that recently opened businesses (last 12 months) were present and the details all looked accurate.

The places data is licensed under “Community Data License Agreement – Permissive” which looks like the only restriction is that you have to include that license when you further share the data.

Tags: geospatial, open-data, parquet, meta, overture

GPSJam

2022-07-30T19:51:31+00:00

GPSJam

John Wiseman’s “Daily maps of GPS interference” —a beautiful interactive globe (powered by Mapbox GL) which you can use to see points of heaviest GPS interference over a 24 hour period, using data collected from commercial airline radios by ADS-B Exchange. “From what I can tell the most common reason for aircraft GPS systems to have degraded accuracy is jamming by military systems. At least, the vast majority of aircraft that I see with bad GPS accuracy are flying near conflict zones where GPS jamming is known to occur.”

Via @lemonodor

Tags: geospatial, gps, mapping, john-wiseman

A tiny web app to create images from OpenStreetMap maps

2022-06-12T05:49:35+00:00

Earlier today I found myself wanting to programmatically generate some images of maps.

I wanted to create a map centered around a location, at a specific zoom level, and with a marker in a specific place.

Some cursory searches failed to turn up exactly what I wanted, so I decided to build a tiny project to solve the problem, taking advantage of my shot-scraper tool for automating screenshots of web pages.

The result is map.simonwillison.net - hosted on GitHub Pages from my simonw/url-map repository.

Here's how to generate a map image of Washington DC:

shot-scraper 'https://map.simonwillison.net/?q=washington+dc' \
  --retina --width 600 --height 400 --wait 3000

That command generates a PNG 1200x800 image that's a retina screenshot of the map displayed at https://map.simonwillison.net/?q=washington+dc - after waiting three seconds to esure all of the tiles have fully loaded.

The website itself is documented here. It displays a map with no visible controls, though you can use gestures to zoom in and pan around - and the URL bar will update to reflect your navigation, so you can bookmark or share the URL once you've got it to the right spot.

You can also use query string parameters to specify the map that should be initially displayed:

https://map.simonwillison.net/?center=51.49,0&zoom=8 displays a map at zoom level 8 centered on the specified latitude, longitude coordinate pair.
https://map.simonwillison.net/?q=islington+london geocodes the ?q= text using OpenStreetMap Nominatim and zooms to the level that best fits the bounding box of the first returned result.
https://map.simonwillison.net/?q=islington+london&zoom=12 does that but zooms to level 12 instead of using the best fit for the bounding box
https://map.simonwillison.net/?center=51.49,0&zoom=8&marker=51.49,0&marker=51.3,0.2 adds two blue markers to the specified map. You can pass &marker=lat,lon as many times as you like to add multiple markers.

Annotated source code

The entire mapping application is contained in a single 68 line index.html file that mixes HTML and JavaScript. It's built using the fantastic Leaflet open source mapping library.

Since the code is so short, I'll enclude the entire thing here with some additional annotating comments.

It started out as a copy of the first example in the Leaflet quick start guide.

<!DOCTYPE html>
<!-- Regular HTML boilerplate -->
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>map.simonwillison.net</title>
<!--
  Leaflet's CSS and JS are loaded from the unpgk.com CDN, with the
  Subresource Integrity (SRI) integrity="sha512..." attribute to ensure
  that the exact expected code is served by the CDN.
-->
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.8.0/dist/leaflet.css" integrity="sha512-hoalWLoI8r4UszCkZ5kL8vayOGVae1oxXe/2A4AO6J9+580uKHDO3JdHb7NzwwzK5xr/Fs0W40kiNHxM9vyTtQ==" crossorigin=""/>
<script src="https://unpkg.com/leaflet@1.8.0/dist/leaflet.js" integrity="sha512-BB3hKbKWOc9Ez/TAwyWxNXeoV9c1v6FIeYiBieIWkpLjauysF18NzgR1MBNBXf8/KABdlkX68nAhlwcDFLGPCQ==" crossorigin=""></script>
<!-- I want the map to occupy the entire browser window with no margins -->
<style>
html, body {
  height: 100%;
  margin: 0;
}
</style>
</head>
<body>
<!-- The Leaflet map renders in this 100% high/wide div -->
<div id="map" style="width: 100%; height: 100%;"></div>
<script>
function toPoint(s) {
  // Convert "51.5,2.1" into [51.5, 2.1]
  return s.split(",").map(parseFloat);
}
// An async function so we can 'await fetch(...)' later on
async function load() {
  // URLSearchParams is a fantastic browser API - it makes it easy to both read
  // query string parameters from the URL and later to generate new ones
  let params = new URLSearchParams(location.search);
  // If the starting URL is /?center=51,32&zoom=3 this will pull those values out
  let center = params.get('center') || '0,0';
  let initialZoom = params.get('zoom');
  let zoom = parseInt(initialZoom || '2', 10);
  let q = params.get('q');
  // .getAll() turns &marker=51.49,0&marker=51.3,0.2 into ['51.49,0', '51.3,0.2']
  let markers = params.getAll('marker');
  // zoomControl: false turns off the visible +/- zoom buttons in Leaflet
  let map = L.map('map', { zoomControl: false }).setView(toPoint(center), zoom);
  L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 19,
    attribution: '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>',
    // This option means retina-capable devices will get double-resolution tiles:
    detectRetina: true
  }).addTo(map);
  // We only pay attention to ?q= if ?center= was not provided:
  if (q && !params.get('center')) {
    // We use fetch to pass ?q= to the Nominatim API and get back JSON
    let response = await fetch(
      `https://nominatim.openstreetmap.org/search.php?q=${encodeURIComponent(q)}&format=jsonv2`
    )
    let data = await response.json();
    // data[0] is the first result - it has a boundingbox array of four floats
    // which we can convert into a Leaflet-compatible bounding box like this:
    let bounds = [
      [data[0].boundingbox[0],data[0].boundingbox[2]],
      [data[0].boundingbox[1],data[0].boundingbox[3]]
    ];
    // This sets both the map center and zooms to the correct level for the bbox:
    map.fitBounds(bounds);
    // User-provided zoom over-rides this
    if (initialZoom) {
      map.setZoom(parseInt(initialZoom));
    }
  }
  // This is the code that updates the URL as the user pans or zooms around.
  // You can subscribe to both the moveend and zoomend Leaflet events in one go:
  map.on('moveend zoomend', () => {
    // Update URL bar with current location
    let newZoom = map.getZoom();
    let center = map.getCenter();
    // This time we use URLSearchParams to construct a center...=&zoom=... URL
    let u = new URLSearchParams();
    // Copy across ?marker=x&marker=y from existing URL, if they were set:
    markers.forEach(s => u.append('marker', s));
    u.append('center', `${center.lat},${center.lng}`);
    u.append('zoom', newZoom);
    // replaceState() is a weird API - the third argument is the one we care about:
    history.replaceState(null, null, '?' + u.toString());
  });
  // This bit adds Leaflet markers to the map for ?marker= query string arguments:
  markers.forEach(s => {
    L.marker(toPoint(s)).addTo(map);
  });
}
load();
</script>
</body>
</html>
<!-- See https://github.com/simonw/url-map for documentation -->

Tags: geospatial, maps, openstreetmap, projects, shot-scraper, leaflet

Datasette for geospatial analysis

2022-04-13T00:48:14+00:00

Datasette for geospatial analysis

I added a new page to the Datasette website describing how Datasette can be used for geospatial analysis, pulling together several of the relevant plugins and tools from the Datasette ecosystem.

Via @simonw

Tags: geospatial, plugins, datasette