Another example showing how to simplify (generalize) a multipolygon layer using PostGIS Topology ¶

This article is a kind of "follower" for this one: UsersWikiSimplifyPreserveTopology, showing another method to build generalized polygons keeping their spatial relationship.

It uses the PostGIS Topology module (​http://postgis.org/documentation/manual-2.0/Topology.html) and aims at producing

A big thank to Brent Wood for its Topology Example and to Sandro Santilli for the topology module

What we want ¶

Simplifying this layer:	into this:	avoiding that:

The data ¶

French administrative subdivisions, called "départements", will be used. Data can be downloaded here: ​http://professionnels.ign.fr/DISPLAY/000/528/175/5281750/GEOFLADept_FR_Corse_AV_L93.zip

Data projection is Lambert-93, EPSG:2154

Loading the data ¶

shp2pgsql -IiD -g geom -s 2154 DEPARTEMENT.SHP departement | psql

Principle of simplification ¶

• The ​| Postgis Topology module is used to build a topology from input data
• A new topology is then created, using a simplified (generalized) version of the first one.
• Node connections seem to be preserved doing so,
• Then, simplified polygons are matched against initial ones using spatial operations

Steps ¶

1. Creates the target table that will contain simplified geometry. Initial departement polygons are dumped to handle simple objects.

create table new_dept as (
       select gid, code_dept, (st_dump(geom)).*
       from departement
);
create index new_dept_geom_gist on new_dept using gist(geom);

-- adds the new geom column that will contain simplified geoms
alter table new_dept add column simple_geom geometry(POLYGON, 2154);

2. Creates a topology from the departements

-- create new empty topology structure
select CreateTopology('topo1',2154,0);

-- add all departements polygons to topology in one operation as a collection
select ST_CreateTopoGeo('topo1',ST_Collect(geom))
from departement;

3. Create a new topology based on the simplification of existing one. (should not be the right way to do it, but calling ST_ChangeEdgeGeom)

select CreateTopology('topo2',2154,0);

select ST_CreateTopoGeo('topo2', geom)
from (
       select ST_Collect(st_simplifyPreserveTopology(geom, 10000)) as geom
       from topo1.edge_data
) as foo;

4. Retrieves polygons by comparing surfaces (pip is not enough for odd-shaped polygons)

with simple_face as (
       select st_getFaceGeometry('topo2', face_id) as geom
       from topo2.face
       where face_id > 0
) update new_dept d set simple_geom = sf.geom
from simple_face sf
where st_intersects(d.geom, sf.geom)
and st_area(st_intersection(sf.geom, d.geom))/st_area(sf.geom) > 0.5;

Results ¶

This method produces a clean result:

Notes ¶

Find a way to associate faces to initial polygon using topology, not using spatial operations.