I’ve been using Postgresql for years now, those who know me, know fine well that it’s my absolute favorite database in the whole wide world. Those who know me may also recall I do a lot of GIS and mapping work using Spatial SQL among other things.
For a recent project I was working on, I was required to model a 3D landscape/terrain of a given location.
This landscape was then turned into a 3D model complete with buildings, navigation controls etc, and all wrapped up in a neat EXE file that allowed folks to walk around the landscape and look at the buildings. The project is far from complete however, and is likely to be ongoing for sometime as there are still lots of 3D modeling of assets to be completed
The interesting part so far however, was modeling the terrain.
First and foremost, getting the general lay of the land, was it has to be said quite easy, using the UK Ordnance Survey Opendata “Terrain 50” data set, I extracted the contours clipped to the geographic area I was interested in, leaving me with a nice neat table containing a bunch of line geometry’s with associated height values
I then used the Quantum GIS interpolation tools to create a smooth interpolated greyscale map of the area, with white being the high ground, and black being the lowground, giving me the following image:
This worked well for the overall landscape, giving me the valleys and the hills in a 10 square mile area, but there was one thing that didn’t quite get shaped as expected
If you look at the base map I’m digitizing, you’ll see there are a number of man made embankments on it. Most of these represent steep drop offs, such as embankments surrounding bridges, or cuttings through which railway lines are built.
Using the OS Contours the resolution is not fine enough to be able to capture these embankments easily, some of which are only half a meter in height, and since the object of modeling the terrain was to produce everything as accurately as possible, this meant I had to find a way to represent these embankments.
Performing a few measurements over the height map I’d already created, did show me that there was a height difference along most of the lines, but because very few of them actually cut through the contours, the height values at these locations where interpolated values over 50 meters or more. This got me thinking, that if I could create some extra contours along the edges of the man made embankments, then set the height value on those lines to a suitable height (In this case an average along the length of the line) I might be able to force the interpolation tool to create some extra ridges, the only question now was how to achieve it.
Postgres to the Rescue
Prior to doing this task, I’d been doing some work on writing my own map server, and I’d been learning how to use PostGIS raster. PostGIS Raster is part of the ever awesome spatial toolkit for PostgreSQL, and if you do any kind of geographic processing at all, it’s a tool that absolutely needs to be in your toolbox. I started to wonder, if I’d be able to load my height map into Postgres and use PostGIS to perform some maths on it, first though I needed to create some new contours
I set about using Quantum GIS to create 4 test contours around embankments I had on the base map:
Once I had these lines, I just needed to figure out how to get the points along the line. I’d already added my height map to Postgress using the “raster2pgsql” tool that comes with Postgres as follows:
raster2pgsql -t 200x200 -s 27700 -q -I -M heightmap.tif heightmaptable > heightmap.sql psql -U gisuser < heightmap.sql
What this does is to take the tif image of my height map as exported from Quantum GIS, then chop it up into 200 pixel square tiles, that are georeferenced to the UK “OSGB36” mapping coordinate system (EPSG code: 27700) , the sql file that’s produced is then run through the postgres command line client and executed on the server to create a table of small 200×200 raster images representing geographically placed portions of the height map.
You create things using tiles in this fashion, simply because it’s easier to scan over a bunch of heavily indexed smaller raster’s, than it is to scan one huge monster raster with no index.
The embankment contours meanwhile where inserted into a geospatial enabled vector table, that holds the line geometry, and an attribute (initially set to 0) that represents the height of the line, the table definitions for both, look like this:
The Embankment Vectors
CREATE TABLE embankmentcontours ( pkid serial NOT NULL, height numeric(6,2), geometry geometry(LineString,27700), CONSTRAINT embankmentcontours_pkey PRIMARY KEY (pkid) ) WITH ( OIDS=FALSE ); ALTER TABLE embankmentcontours OWNER TO gismapuser;
The Height map rasters
CREATE TABLE heightmap ( pkid integer serial NOT NULL, rast raster, CONSTRAINT heightmap_pkey PRIMARY KEY (pkid) ) WITH ( OIDS=FALSE ); ALTER TABLE heightmap OWNER TO gismapuser; CREATE INDEX heightmap_st_convexhull_idx ON heightmap USING gist (st_convexhull(rast));
Now I had all my raw materials in the database, all that was left was to figure out some postgres magic.
One thing that Postgres (IMHO anyway) does, far better than any other database out there, is a construct called a “Common Table Expression”, for those who’ve never seen “CTE’s” before, it’s a bit like defining a temporary table, but one that lasts only for the duration of the query, something like the following:
with mycte as ( select something from somewhere where something = somethingelse )
You can then perform normal queries over the result of that query:
with mycte as ( select something from somewhere where something = somethingelse ) select * from mycte
The best part is, you can chain cte’s to form more complex queries:
with mycte1 as ( select 'hello'::text as hello, 'world'::text as world ), mycte2 as ( select hello || ' ' || world as greeting from mycte1 ) select * from mycte2
Using a CTE, we could easily break down the line geometry representing our embankment lines with the PostGIS function “ST_DumpPoints”:
We can then, cross reference each of these points using a join on our raster table, to find out which raster tile they fall in, and extract the height of the point at that location:
Finally, we average out the values of the points, grouped by the ID of the embankment contour they came from, so we get a median value representing a leveled value accross that contour line.
Beacuse we did this using a chain of 3 CTE’s, we just need one simple update statement that uses the data derived from the last link in the CTE chain to update the height field of the original contour. Since the average height data also includes the original embankment contour ID, it’s as simple as just setting the height value on a contour line, where the id matches the currently accessed value, the full Postgres script to perform the entire task is thus:
with points as ( select contourid, (dp).path as pointindex, (dp).geom as pointgeometry from ( select pkid as contourid, ST_DumpPoints(geometry) as dp from embankmentcontours ) as pointlist ), pointswithrasterid as ( select pt.contourid, pt.pointindex, hm.pkid as rasterid, pt.pointgeometry, hm.rast as raster from points as pt left join heightmap as hm on ST_Intersects(hm.rast, pt.pointgeometry) ), contouraverages as ( select contourid, avg(ST_Value(raster, pointgeometry)) as averageheight from pointswithrasterid group by contourid ) update embankmentcontours cntrs set height = cte3.averageheight from (select * from contouraverages order by contourid) as cte3 where cntrs.pkid = cte3.contourid
If we now look at the data in our embankment contours table, we’ll see that each line now has an associated median height:
Will this actually work to help me shape the embankments, well that we’ll find out once I’ve sat and drawn all the contour lines in, then re-run the interpolation job in Qauntum, for now though, as the title states….
PostgreSQL never ceases to amaze me, the things you can do in a fairly simple single script like the one above put most other DBMS systems to shame