George MacKerron: code blog

GIS, software development, and other snippets

Polygons from PostGIS to Processing

There are plenty of ways to get spatial data from a PostGIS database into a Processing sketch.

You can export to CSV or SVG and load it from there; you can query the database directly; or, depending on context, you might choose to generate Processing commands directly, which is the route I went to display a background map of the UK in a recent visualization project.

First, I load outlines of the 4 constituent countries of the UK into my PostGIS database.

shp2pgsql -d -D -I -s 27700 "/pgdata/geodata/Nations/england.shp" eng \
  | sudo -u postgres psql mappiness
shp2pgsql -d -D -I -s 27700 "/pgdata/geodata/Nations/wales.shp" wal \
  | sudo -u postgres psql mappiness
shp2pgsql -d -D -I -s 27700 "/pgdata/geodata/Nations/scotland.shp" sco \
  | sudo -u postgres psql mappiness
shp2pgsql -d -D -I -s 27700 "/pgdata/geodata/Nations/nireland_ol_2001.shp" ni \
  | sudo -u postgres psql mappiness

Next, I combine these into one table, using st_simplifypreservetopology to discard a lot of unnecessary detail and st_dump to dissolve some English MULTIPOLYGONs into simple POLYGONS.

create table uk1km as (
        select st_simplifypreservetopology(the_geom, 1000) as the_geom from 
          (select (st_dump(the_geom)).geom as the_geom from eng) as eng_single
  union select st_simplifypreservetopology(the_geom, 1000) as the_geom from wal
  union select st_simplifypreservetopology(the_geom, 1000) as the_geom from sco
  union select st_simplifypreservetopology(the_geom, 1000) as the_geom from ni
);
alter table uk1km add column id serial primary key;
drop table eng; drop table wal; drop table sco; drop table ni;

Finally, I use string_agg to produce the vertex() commands required to reproduce my polygons in Processing. I change the sign on the y-coordinate (from northings to ‘southings’, to match the Processing coordinate space). I scale the OSGB36 coordinates down from metres to kilometres (I could otherwise have done this with a scale() command in Processing). And I discard a lot of small polygons (comprising 20 POINTs or fewer) I don’t need in my output.

select 'beginShape(); ' 
    || string_agg('vertex(' || st_x(geom) || ', ' || -st_y(geom) || '); ', '')
    || 'endShape(CLOSE); ' 
from (
  select id, (st_dumppoints(st_exteriorring(st_scale(the_geom, 0.001, 0.001)))).geom as geom 
  from uk1km 
  where st_numpoints(st_exteriorring(the_geom)) > 20
) as d
group by id;

Then, to draw the UK in a sketch I do this.

size(350, 620);
translate(10, height);
scale(0.5);
 
// paste in the output of the previous query

Which has this result (with thanks to EDINA for the data):

The UK as drawn by this method

Share

Written by George

November 17th, 2011 at 12:28 pm

Posted in GIS,PostGIS,Processing