Version 8 (modified by 12 years ago) ( diff ) | ,
---|
Another example showing how to simplify (generalize) a multipolygon layer using PostGIS Topology
This article is a sort of "follow-up" for this one: UsersWikiSimplifyPreserveTopology, showing another method to build generalized polygons keeping their spatial relationship, based on spatial operators and SQL queries.
Another take is on http://strk.keybit.net/blog/2012/04/13/simplifying-a-map-layer-using-postgis-topology/
This example uses the PostGIS Topology module (http://postgis.org/documentation/manual-2.0/Topology.html) and aims at producing a simplified polygonal layer, with initial polygons' attributes correctly associated to simplified polygons
(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
- 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);
- 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;
- 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;
- 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, departement's codes are correcty re-associated:
Notes
Find a way to associate faces to initial polygon using topology, not using spatial operations. see strk mail
Attachments (1)
- SPT_simple_good_topo.png (78.3 KB ) - added by 12 years ago.
Download all attachments as: .zip