Simon Willison's Weblog: shapefiles

mapshaper.org

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

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

geoBoundaries

2022-03-24T14:03:56+00:00

geoBoundaries

This looks useful: “The world’s largest open, free and research-ready database of political administrative boundaries.” Founded by the geoLab at William & Mary university, and released under a Creative Commons Attribution license that includes a requirement for a citation. File formats offered include shapefiles, GeoJSON and TopoJSON.

Tags: geospatial, shapefiles, geojson

Weeknotes: datasette-export-notebook, PyInstaller packaged Datasette, CBSAs

2021-01-10T00:26:13+00:00

What a terrible week. I've found it hard to concentrate on anything substantial. In a mostly futile attempt to distract myself from doomscrolling I've mainly been building some experimental output plugins, fiddling with PyInstaller and messing around with shapefiles.

Packaged Datasette with PyInstaller

A long running goal for Datasette has been to make it as easy to install as possible - something that's not particularly straight-forward for applications written in Python, at least in comparison to toolchains like Rust, Go or Deno.

Back in November 2017 Raffaele Messuti suggested using PyInstaller for this. I revisited that issue while looking through open Datasette issues ordered by least recently updated and decided to try it out again - and it worked! Here's the resulting TIL, and I've attached a bundled datasette macOS binary file to the 0.53 release on GitHub.

There's one catch: the binary isn't signed, which means it shows security warnings that have to be worked around if you try to run it on macOS. I've started looking into the signing process - I'm going to need an Apple Developer account and to jump through a bunch of different hoops, but it looks like I should be able to get that working. Here's my issue for that: GitHub Actions workflow to build and sign macOS binary executables - it looks like gon is the missing automation piece I need.

One thing that really impressed me about PyInstaller is the size of the resulting file. On both macOS and Linux it was able to create a roughly 8MB file containing Datasette, all of its dependencies AND a working Python environment. It's pretty magic!

datasette-css-properties

I wrote this up in detail a couple of days ago: datasette-css-properties is an amusingly weird output plugin that turns the results of a SQL query into CSS custom property definitions which can then be used to style or insert content into the current page.

sqlite-utils 3.2

The big new feature in this release is cached table counts using triggers, which I described last week. Full release notes here.

I've opened an issue to take advantage of this optimization in Datasette itself.

datasette-export-notebook

This is an idea I've been bouncing around for a while, and during a bout of attempted-coup-induced insomnia I decided to sketch out an initial version.

datasette-export-notebook is a plugin that adds an export-to-notebook option to any table or query.

This provides a page of documentation with copy-and-paste examples for loading data from that table or query into a Jupyter or Observable notebook.

Here's a live demo. The interface currently looks like this:

As often happens when building even simple plugins like this I identified some small improvements I can make to Datasette.

cbsa-datasette

Core-based statistical areas are a US government concept used for various statistical purposes. They are essentially metropolitan areas, based on central cities and the commuting area that they sit inside.

I built cbsa.datasettes.com this week to provide an API for looking up a CBSA based on a latitude and longitude point. Here's a location within San Francisco for example:

https://cbsa.datasettes.com/core/by_lat_lon?longitude=-122.51&latitude=37.78

This returns San Francisco-Oakland-Berkeley, CA. Add .json and &_shape=array to the above URL to get a JSON API version.

The data comes from a shapefile published by the Bureau of Transportation Stastics. I'm using shapefile-to-sqlite to import it into a SpatiaLite database, then publishing it to Cloud Run using this GitHub Actions workfow. Full details in the README.

I built this mainly to act as a simple updated example of how to use Datasette and SpatiaLite to provide an API against data from a shapefile. I published a tutorial about doing this for timezones three years ago, but shapefile-to-sqlite makes it much easier.

Releases this week

datasette-export-notebook: 0.1.1 - 2021-01-09
Datasette plugin providing instructions for exporting data to Jupyter or Observable
datasette-css-properties: 0.2 - 2021-01-07
Experimental Datasette output plugin using CSS properties
sqlite-utils: 3.2 - 2021-01-03
Python CLI utility and library for manipulating SQLite databases

TIL this week

Packaging a Python app as a standalone binary with PyInstaller

Tags: geospatial, projects, python, shapefiles, datasette, weeknotes, sqlite-utils

California Protected Areas Database in Datasette

2020-08-21T23:15:58+00:00

California Protected Areas Database in Datasette

I built this yesterday: it’s a Datasette interface on top of the CPAD 2020 GIS database of protected areas in California maintained by GreenInfo Network. This was a useful excuse to build a GitHub Actions flow that builds a SpatiaLite database using my shapefile-to-sqlite tool, and I fixed a few bugs in my datasette-leaflet-geojson plugin as well.

Via calands-datasette on GitHub

Tags: geospatial, projects, shapefiles, spatialite, california, datasette, github-actions, leaflet

Things I learned about shapefiles building shapefile-to-sqlite

2020-02-19T05:25:58+00:00

The latest in my series of x-to-sqlite tools is shapefile-to-sqlite. I learned a whole bunch of things about the ESRI shapefile format while building it.

Governments really love ESRI shapefiles. There is a huge amount of interesting geospatial data made available in the format - 4,614 on Data.gov!

shapefile-to-sqlite

shapefile-to-sqlite loads the data from these files into a SQLite database, turning geometry properties into database columns and the geometry itself into a blob of GeoJSON. Let's try it out on a shapefile containing the boundaries of US national parks.

$ wget http://nrdata.nps.gov/programs/lands/nps_boundary.zip
...
Saving to: ‘nps_boundary.zip’
nps_boundary.zip                           100%[=====================================================================================>]  12.61M   705KB/s    in 22s     
2020-02-18 19:59:22 (597 KB/s) - ‘nps_boundary.zip’ saved [13227561/13227561]

$ unzip nps_boundary.zip 
Archive:  nps_boundary.zip
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.xml  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.dbf  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.prj  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shx

$ shapefile-to-sqlite nps.db temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
[####################################]  100%

$ datasette nps.db
Serve! files=('nps.db',) (immutables=()) on port 8003
INFO:     Started server process [33534]
INFO:     Waiting for application startup.
INFO:     Application startup complete.
INFO:     Uvicorn running on http://127.0.0.1:8001 (Press CTRL+C to quit)

I recommend installing the datasette-leaflet-geojson plugin, which will turn any column containing GeoJSON into a Leaflet map.

If you've installed SpatiaLite (installation instructions here) you can use the --spatialite option to instead store the geometry in a SpatiaLite column, unlocking a bewildering array of SQL geometry functions.

$ shapefile-to-sqlite nps.db temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp --spatialite --table=nps-spatialite
temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
[##################################--]   94%  00:00:00

I deployed a copy of the resulting database using Cloud Run:

$ datasette publish cloudrun nps.db \
    --service national-parks \
    --title "National Parks" \
    --source_url="https://catalog.data.gov/dataset/national-parks" \
    --source="data.gov" \
    --spatialite \
    --install=datasette-leaflet-geojson \
    --install=datasette-render-binary \
    --extra-options="--config max_returned_rows:5"

I used max_returned_rows:5 there because these geometrries are pretty big - without it a page with 100 rows on it can return over 90MB of HTML!

You can browse the GeoJSON version of the table here and the SpatiaLite version here.

The SpatiaLite version defaults to rendering each geometry as an ugly binary blob. You can convert them to GeoJSON for compatibility with datasette-leaflet-geojson using the SpatiaLite AsGeoJSON() function:

select id, UNIT_NAME, AsGeoJSON(geometry)
from [nps-spatialite]

Here's the result of that query running against the demo.

Understanding shapefiles

The most confusing thing about shapefiles is that they aren't a single file. A shapefile comes as a minimum of three files: foo.shp containing geometries, foo.shx containing an index into those geometries (really more of an implementation detail) and foo.dbf contains key/value properties for each geometry.

They often come bundled with other files too. foo.prj is a WKT projection for the data for example. Wikipedia lists a whole bunch of other possibilities.

As a result, shapefiles are usually distributed as a zip file. Some shapefile libraries can even read directly from a zip.

The GeoJSON format was designed as a modern alternative to shapefiles, so understanding GeoJSON really helps in understanding shapefiles. In particular the GeoJSON geometry types: Point, LineString, MultiLineString, Polygon and MultiPolygon match how shapefile geometries work.

An important detail in shapefiles is that data in the .shp and .dbf files is matched by array index - so the first geometry can be considered as having ID=0, the second ID=1 and so on.

You can read the properties from the .dbf file using the dbfread Python module like this:

$ ipython
In [1]: import dbfread
In [2]: db = dbfread.DBF("temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.dbf")
In [3]: next(iter(db))
Out[3]: 
OrderedDict([('UNIT_TYPE', 'Park'),
            ('STATE', ''),
            ('REGION', 'NC'),
            ('UNIT_CODE', 'NACC'),
            ('UNIT_NAME', 'West Potomac Park'),
            ('DATE_EDIT', None),
            ('GIS_NOTES', ''),
            ('CREATED_BY', 'Legacy'),
            ('METADATA', ''),
            ('PARKNAME', '')])

Reading shapefiles in Python

I'm a big fan of the Shapely Python library, so I was delighted to see that Sean Gillies, creator of Shapely, also created a library for reading and writing shapefiles: Fiona.

GIS with Python, Shapely, and Fiona by Tom MacWright was particularly useful for figuring this out. I like how he wrote that post in 2012 but added a note in 2017 that it's still his recommended way of getting started with GIS in Python.

Projections

The trickiest part of working with any GIS data is always figuring out how to deal with projections.

GeoJSON attempts to standardize on WGS 84, otherwise known as the latitude/longitude model used by GPS. But... shapefiles frequently use something else. The Santa Clara county parks shapefiles for example use EPSG:2227, also known as California zone 3.

(Fun fact: ESPG stands for European Petroleum Survey Group, a now defunct oil industry group that today lives on only as a database of projected coordinate systems.)

I spent quite a while thinking about how to best handle projections. In the end I decided that I'd follow GeoJSON's lead and attempt to convert everything to WGS 84, but allow users to skip that behaviour using --crs=keep or to specify an alternative projection to convert to with --crs=epsg:2227 or similar.

SpatiaLite creates its geometry columns with a baked in SRID (a code which usually maps to the EPSG identifier). You can see which SRID was used for a specific geometry using the srid() function:

select srid(geometry) from "nps-spatialite" limit 1

SpatiaLite can also convert to another projection using the Transform() function:

select ':' || AsGeoJSON(Transform(geometry, 2227)) from "nps-spatialite" limit 1

(I'm using ':' || AsGeoJSON(...) here to disable the datasette-leaflet-geojson plugin, since it can't correctly render data that has been transformed to a non-WGS-84 proection.)

Pulling it all together

I now have two tools for imorting geospatial data into SQLite (or SpatiaLite) databases: shapefile-to-sqlite and geojson-to-sqlite.

I'm excited about Datasette's potential as a tool for GIS. I started exploring this back in 2017 when I used it to build a location to timezone API - but adding easy shapefile imports to the toolchain should unlock all kinds of interesting new geospatial projects.

Tags: geospatial, projects, shapefiles, spatialite, sqlite, geojson, weeknotes, tom-macwright, leaflet

Flickr Shapefiles Public Dataset 1.0

2009-05-22T18:12:10+00:00

Flickr Shapefiles Public Dataset 1.0

Another awesome Geo dataset from the Yahoo! stable—this time it’s Flickr releasing shapefiles (geometrical shapes) for hundreds of thousands of places around the world, under the CC0 license which makes them essentially public domain. The shapes themselves have been crowdsourced from geocoded photos uploaded to Flickr, where users can “correct” the textual location assigned to each photo. Combine this with the GeoPlanet WOE data and you get a huge, free dataset describing the human geography of the world.

Tags: creativecommons, crowdsourcing, flickr, geoplanet, geospatial, maps, shapefiles, yahoo