Version 1 (modified by 12 years ago) ( diff ) | ,
---|
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
- 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:
Notes
Find a way to associate faces to initial polygon using topology, not using spatial operations.
Attachments (1)
- SPT_simple_good_topo.png (78.3 KB ) - added by 12 years ago.
Download all attachments as: .zip