June 7, 2010

Glacier Contours

I consider the USGS 7.5' quads to be among the best topographic maps around. That said, working around their shortcomings is the primary cartographic motivation for Closed Contour. Right now I feel like I have a ways to go before I can claim that my maps are 'better' (how ever you define it) than the USGS quads. Really this is a long-winded way of saying that I'm going to blatantly steal ideas and styling from the USGS quads until they don't have anything more to offer!

Here's a prime example: contour lines on glaciers (or, in California, permanent snow). Take a look at the area around North Palisade:

And now the same area as rendered in the v1 SPS Map tiles:

See the problem? On the USGS map, contour lines on the Palisade Glacier are light blue, as is the outline. Even though the glacier is white on the Closed Contour map you can barely see it. I decided that my maps needed this feature, both for usability and aesthetics.

Another motivation for making this map was to learn some new technologies. Sure, I'd written a little SQL before and I'm pretty familiar with geospatial software and techniques, but I hadn't ever really used PostgreSQL and/or PostGIS. This would be a great opportunity to try to do something useful purely in SQL.

A little background: at this point all of the map data (except the hillshade raster) are stored in a haphazardly organized PostgreSQL with PostGIS database; each data source with it's own table. Actually, in the case of contours I have them separated into 3 tables, one for 40', one for 200', and one for 1000' intervals. This is an artifact of the contour generation/import process and something I'll hopefully clean up. There's also a polygon land-cover table.

So all I need to do is to iterate over the snow/glacier land-cover polygons, find all contours inside them, slice them to be exactly inside the snow and mark them as such, and then subtract out the contours from the original data set. No problem, right? Here's how I accomplished that in SQL:

CREATE OR REPLACE FUNCTION find_snow_contours()
    RETURNS  void AS $$
DECLARE
    contour_line RECORD;
    snow_poly RECORD;
    g geometry;
BEGIN
    -- Holding place for contours that are purely contained by snow/glacier
    DROP TABLE IF EXISTS temp_snow_contours;
    CREATE TABLE temp_snow_contours (like contours_40);

    -- Iterate over all snow polygons ('snow' is a temporary table with just snow polygons')
    FOR snow_poly IN SELECT * FROM snow LOOP
    -- Iterate over all contour lines that intersect this snow polygon
    FOR contour_line IN SELECT * from contours_40 WHERE ST_Intersects(the_geom, snow_poly.the_geom) LOOP

        -- Subtract off the snow section and save resulting lines as non-snow contours
        SELECT INTO g ST_Difference(contour_line.the_geom, snow_poly.the_geom);

        -- Delete the original contour line and replace it with the remaining segments
        -- It's common for there to be many segments after the difference, hence the 'ST_Dump'
        DELETE FROM contours_40 WHERE gid = contour_line.gid;
        INSERT INTO contours_40 (elev, is_snow, the_geom)
            SELECT contour_line.elev, 0, (ST_Dump(g)).geom;

        -- Add the intersection as snow contours
        SELECT INTO g ST_Intersection(contour_line.the_geom, snow_poly.the_geom);
        INSERT INTO temp_snow_contours (gid, elev, is_snow, the_geom) VALUES (contour_line.gid, contour_line.elev, 1, g);

    END LOOP;
    END LOOP;

    -- Finally put the snow contours back in the main table.
    INSERT INTO contours_40 (elev, the_geom, is_snow)
    SELECT  elev, (ST_Dump(the_geom)).geom, is_snow FROM temp_snow_contours;
    DROP TABLE temp_snow_contours;
    RETURN;
END;
$$
    

I imagine that a PostGIS SQL-guru would look at this function and shudder... forgive me, I'm a newb. Anyway, this technique worked out nicely and I'm pretty happy with the result:

Not sure when I'll generate them but I'll definitely be rolling this feature into v2 of the tiles.

Back