Postgres never ceases to amaze me


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:

heightmap

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:

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

Which produces:

pgadmin

Using a CTE, we could easily break down the line geometry representing our embankment lines with the PostGIS function “ST_DumpPoints”:

pgadmin2

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:

pgadmin3

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.

pgadmin4

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[1] 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:

pgadmin5

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

Happy GISing :-), Shawty
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s