- PostGIS Cookbook
- Paolo Corti Thomas J. Kraft Stephen Vincent Mather Bborie Park
- 921字
- 2021-07-19 18:29:44
Normalizing internal overlays
Data from an external source can have not just table structure issues, but also topological issues endemic to the geospatial data itself. Take, for example, the problem of data with overlapping polygons. If our dataset has polygons that overlap with internal overlays, queries for area, perimeter, and other metrics may not produce predictable or consistent results.
There are a few approaches that can solve the problem of polygon datasets with internal overlays. The general approach presented here was originally proposed by Kevin Neufeld of Refractions Research.
Over the course of writing our query, we will also produce a solution for converting polygons to linestrings.
Getting ready
First, we'll load our dataset using the following command:
shp2pgsql -s 3734 -d -i -I -W LATIN1 -g the_geom cm_usearea_polygon chp02.use_area | psql -U me -d postgis_cookbook
How to do it...
Now that the data is loaded into a table in the database, we can leverage PostGIS to flatten and get the union of the polygons, such that we have a normalized dataset. The first step in doing so using this approach will be to convert the polygons to linestrings. We can then node those linestrings and convert them back to polygons, representing the union of all the polygon inputs. We will perform the following tasks:
- Converting polygons to linestrings.
- Converting linestrings back to polygons.
- Finding center points of resultant polygons.
- Use resultant points to query tabular relationships.
To accomplish this, we'll need to extract just the portions of the polygons we want using ST_ExteriorRing
, convert those parts to points using ST_DumpPoints
, and then connect those points back into lines like a "connect-the-dots" coloring book using ST_MakeLine
.
Breaking it down further, ST_ExteriorRing (the_geom)
will grab just the outer boundary of our polygons. But ST_ExteriorRing
returns polygons, so we need to take that output and create a line from it. The easiest way to do this is to convert it to points using ST_DumpPoints
and then connect those points. By default, the Dump
function returns an object called a geometry_dump
, which is not just simple geometry but the geometry in combination with an array of integers. The easiest way to return the geometry alone is to leverage the object notation to extract just the geometry portion of geometry_dump
as follows:
(ST_DumpPoints(geom)).geom
Piecing the geometry back together with ST_ExteriorRing
is done using the following query:
SELECT (ST_DumpPoints(ST_ExteriorRing(geom))).geom
This should give us a listing of points in order from the exterior rings of all the points from which we want to construct our lines using ST_MakeLine, as shown in the following query:
SELECT ST_MakeLine(geom) FROM ( SELECT (ST_DumpPoints(ST_ExteriorRing(geom))).geom ) AS linpoints
Since the preceding approach is a process we may want to use in many other places, it might be prudent to create a function from this using the following query:
CREATE OR REPLACE FUNCTION chp02.polygon_to_line(geometry) RETURNS geometry AS $BODY$ SELECT ST_MakeLine(geom) FROM ( SELECT (ST_DumpPoints(ST_ExteriorRing( (ST_Dump($1)).geom ))).geom ) AS linpoints $BODY$ LANGUAGE sql VOLATILE; ALTER FUNCTION chp02.polygon_to_line(geometry) OWNER TO me;
Now that we have the polygon_to_line
function, we still need to force the noding of overlapping lines in our particular use. The ST_Union
function will aid in this as shown in the following query:
SELECT ST_Union(geom) AS geom FROM ( SELECT chp02.polygon_to_line(geom) AS geom FROM chp02.use_area ) AS unioned ;
Now we can polygonize this result using ST_Polygonize
, as shown in the following query:
SELECT ST_Polygonize(geom) AS geom FROM ( SELECT ST_Union(geom) AS geom FROM ( SELECT chp02.polygon_to_line(geom) AS geom FROM chp02.use_area ) AS unioned ) as polygonized;
The ST_Polygonize
function will create a single multi polygon, so we need to explode this into multiple single polygon geometries if we are to do anything useful with it. While we are at it, we might as well do the following within a CREATE TABLE
statement:
CREATE TABLE chp02.use_area_alt AS ( SELECT (ST_Dump(the_geom)).geom AS the_geom FROM ( SELECT ST_Polygonize(the_geom) AS the_geom FROM ( SELECT ST_Union(the_geom) AS the_geom FROM ( SELECT chp02.polygon_to_line(the_geom) AS the_geom FROM chp02.use_area ) AS unioned ) as polygonized ) AS exploded );
We will be performing spatial queries against this geometry, so we should create an index in order to ensure our query performs well, as shown in the following query:
CREATE INDEX chp02_use_area_alt_the_geom_gist ON chp02.use_area_alt USING gist(the_geom);
In order to extract the appropriate table information from the original geometry and apply that back to our resultant geometries, we will perform a point-in-polygon query. For that, we first need to calculate centroids on the resultant geometry:
CREATE TABLE chp02.use_area_alt_p AS SELECT ST_SetSRID(ST_PointOnSurface(the_geom), 3734) AS the_geom FROM chp02.use_area_alt; ALTER TABLE chp02.use_area_alt_p ADD COLUMN gid serial; ALTER TABLE chp02.use_area_alt_p ADD PRIMARY KEY (gid);
And, as always, create a spatial index using the following query:
CREATE INDEX chp02_use_area_alt_p_the_geom_gist ON chp02.use_area_alt_p USING gist(the_geom);
The centroids then structure our point-in-polygon (ST_Intersects
) relationship between the original tabular information and resultant polygons, using the following query:
CREATE TABLE chp02.use_area_alt_relation AS SELECT points.gid, cu.location FROM chp02.use_area_alt_p AS points, chp02.use_area AS cu WHERE ST_Intersects(points.the_geom, cu.the_geom);
How it works...
Our essential approach here is to look at the underlying topology of the geometry and reconstruct a topology that is nonoverlapping, and then use the centroids of that new geometry to construct a query that establishes the relationship to the original data.
There's more...
At this stage, we can optionally establish a framework for referential integrity using a foreign key as follows:
ALTER TABLE chp02.use_area_alt_relation ADD FOREIGN KEY (gid) REFERENCES chp02.use_area_alt_p (gid);
- Spring 5.0 Microservices(Second Edition)
- Angular UI Development with PrimeNG
- Rake Task Management Essentials
- 零基礎學Java程序設計
- Backbone.js Blueprints
- 精通Python設計模式(第2版)
- 劍指大數據:企業級數據倉庫項目實戰(在線教育版)
- UML2面向對象分析與設計(第2版)
- 物聯網系統架構設計與邊緣計算(原書第2版)
- iOS開發項目化入門教程
- HTML5移動Web開發
- INSTANT PLC Programming with RSLogix 5000
- 微信公眾平臺開發最佳實踐
- Java無難事:詳解Java編程核心思想與技術
- WCF 4.5 Multi-Layer Services Development with Entity Framework(Third Edition)