I wanted to play with contours for some mapping. Ie: a shapefile of where the 8000′ or higher terrain is across the US. Turns out to be pretty simple to generate something basically working from SRTM data. The gdal_contour program does the work. The OpenStreetMap folks have some nice notes on how to use it with SRTM data. The details, however, are much harder.


$ gdal_contour -a elevation N40W108.hgt N40W108.shp -i 304.8 -snodata 32767

That generates 1000′ contours, very fast, maybe 150ms with the data in cache. The output file size varies wildly: for that hilly part of Colorado it’s 1 meg. Generating 100′ contours instead takes about 6x as long and results in 10x the resulting file size. (gdal_contour is dumb about files: if the output file already exists, the error you get is “ERROR 1: N40W108.shp is not a directory”)

North America has 2415 data files so even on a slow computer 1000′ contours can be generated in about an hour. Once you have the contour shapefiles, you can import them into PostGIS. It goes pretty easily, but is slower than generating the shape files themselves. Should probably create an index, too.

$ shp2pgsql -p N20W076.hgt.shp srtm_contours | psql contours
$ shp2pgsql -a N20W076.hgt.shp srtm_contours | psql contours
$ shp2pgsql -a N20W075.hgt.shp srtm_contours | psql contours
$ shp2pgsql -a N19W075.hgt.shp srtm_contours | psql contours
$ shp2pgsql -a N19W076.hgt.shp srtm_contours | psql contours

The import of all of North America took a couple of hours and the table is 14 gigs on disk. 32,865,353 rows. Some quick efforts to make maps with TileMill suggest it’s too much data. Need indexing on elevation and the geometry, also would really like to simplify those contour lines.

Thinking in feet

The elevations are in meters, stored as a double, but I think in feet and it’d be nice to simplify. So let’s convert them to integer feet and index on elevation.

alter table srtm_contours add column elevation_feet integer;
update srtm_contours set elevation_feet = round(3.28084*elevation);

This is very slow, over an hour. PostgreSQL is completely copying the whole table to make rows with the new column. Maybe it would have been better to update the value in place.

Indexing the data

Let’s create a basic spatial index!

create index idx_srtm_contours_the_geom on srtm_contours using gist(the_geom);  -- about 30 minutes, 2 gigs
vacuum analyze srtm_contours; -- 30+ minutes, at least before database tuning

Subsetting for TileMill

Even with the index the contours are just way, way too much for TileMill. I’m not sure how smart TileMill is about querying only data it needs. A different solution is to create a special table with just rows for the 8000′ elevation level (the one we care about). It also seems to help to simplify the contours some.

create table contour_8000_simple (
  gid integer,
  id integer,
  the_geom geometry,
  constraint contour_8000_pk primary key (gid)

insert into contour_8000_simple
select gid, id, ST_SimplifyPreserveTopology(the_geom, 0.001) from srtm_contours where elevation_feet = 8000;

(SQL is approximate, not cut-and-paste ready). The 0.001 parameter means “simplify to 0.001 degree”, in this case.

Next Steps

After all the selecting and simplifying the resulting data works OK in TileMill. It renders at an acceptable speed, at least, and if you draw big fat lines you can schmudge over the awkward too-detailed contour. But this work has gotten pretty hacky. I really need polygons, not linestrings, so I can fill the 8000’+ range. But because the contours were calculated on 1degree square tiles in the first place, closing the polygons takes a bit of work. Also I haven’t dealt properly with NODATA regions. Need to sit down and re-do it with more time.

Here’s what it looks like, for a region of mountains SE of Salt Lake City, UT.