官术网_书友最值得收藏!

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:

  1. Converting polygons to linestrings.
  2. Converting linestrings back to polygons.
  3. Finding center points of resultant polygons.
  4. Use resultant points to query tabular relationships.

Converting polygons to linestrings

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
;

Converting linestrings back to polygons

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);

Finding center points of resultant polygons

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);

Using resultant points to query tabular relationships

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);
主站蜘蛛池模板: 顺平县| 崇文区| 桂平市| 东方市| 马鞍山市| 佛教| 临海市| 宾阳县| 深泽县| 星子县| 宽城| 辽中县| 阿鲁科尔沁旗| 贵定县| 灯塔市| 华容县| 且末县| 淮安市| 西峡县| 昭通市| 英吉沙县| 九龙城区| 西宁市| 城固县| 兰州市| 汉阴县| 德阳市| 安图县| 太谷县| 辽宁省| 马公市| 沙湾县| 凯里市| 新乡县| 龙江县| 察隅县| 通城县| 柏乡县| 新津县| 大石桥市| 新源县|