One of the problems that arises in the process of drawing lands and lakes boundaries is the problem of determining the relative nesting of land and water areas. It’s a common thing, when there are lakes inside the continent, and, in turn, there can be islands in those lakes. In principle, infinite nesting is possible: there can be pools of water on the islands, in which there may be small islets, and so on.
When the boundaries between land and water are determined the role of these boundaries are not always clear. A boundary line may be the mainland coastline or the shore of a lake. The task here is to strictly define the boundary of which object is the given boundary line.
An example of in-lake island
You can accomplish this task using PostGIS extension, just as we used it to create rivers with variable widths.
As an example, we consider four simple geographical objects as in the picture: red squares are the continents, blue square is the lake, and yellow square is the island in this lake.
It is natural to consider that the lakes’ areas belong to the mainland, and the lake islands’ areas belong to the lakes. That is, all our objects are polygons without holes. In this case, the objects will be located either side by side (if they do not intersect) or one above the other, as if in layers.
Partial overlapping of objects are not allowed (in the project this is provided by geography objects creating procedures).
Let’s create a GIS database of those objects.
CREATE EXTENSION postgis; CREATE TABLE IF NOT EXISTS polygons ( fid Bigint NOT NULL, pfid Bigint, zindex Smallint, area Double Precision, geom Geometry(POLYGON, 3785) );
polygons table has following fields:
fidis feature id of the polygon.
pfidis the id of the polygon container in which the given polygon is directly contained.
zindexis the layer number in which a polygon is located.
areais the area of the polygon. In this example, we calculate the 2D Cartesian area of the polygon; if you need an exact area on the sphere, then you need to go to PostGIS Geography Type.
geomis the field with geographic data in the Mercator projection (
For efficiency we should create some indices.
CREATE INDEX area_idx ON polygons (area); CREATE INDEX pfid_idx ON polygons (pfid); CREATE INDEX polygons_geom_idx ON polygons USING GIST (geom);
pfid_idx are the standard indices on
polygons_geom_idx is the GIST index on
geom field, which speeds up the process of finding the nested polygons. You may prefer to create indices after filling out the table.
arrangePolygons() scans the
polygons table and finds the container polygon id. If there is no container for a polygon then NULL value remains.
CREATE OR REPLACE FUNCTION arrangePolygons() RETURNS Void AS $proc$ DECLARE curs_feats CURSOR FOR SELECT fid, area, geom FROM polygons ORDER BY area FOR UPDATE; curs_parents CURSOR ( feat_area double precision, feat_geom Geometry) FOR SELECT fid, geom FROM polygons WHERE area>feat_area AND geom && feat_geom ORDER BY area; pfidV Bigint; BEGIN SET enable_seqscan TO off; <<feats_loop>> FOR feat IN curs_feats LOOP pfidV = NULL; IF feat.geom IS NOT NULL THEN <<parents_loop>> FOR pfeat IN curs_parents ( feat_area:=feat.area, feat_geom:=feat.geom) LOOP IF _ST_Contains(pfeat.geom, feat.geom)=true THEN pfidV = pfeat.fid; EXIT parents_loop; END IF; END LOOP parents_loop; UPDATE polygons SET pfid = pfidV WHERE CURRENT OF curs_feats; END IF; END LOOP feats_loop; END; $proc$ LANGUAGE plpgsql;
The procedure is quite simple and consists of two nested loops. In
curs_parents cursor we use ordering on
area to find the object’s container: from those polygons that contain the given one the container will, of course, have the smallest area.
Note that we use a variant of the
_ST_Contains function with an underscore, since
ST_Contains could preliminary access the
polygons_geom_idx index, which we have already used with
&& operator on the bounding boxes.
The expression `
SET enable_seqscan TO off;‘ is here for forcing usage of defined indices.
The second procedure
setZIndex(Bigint, Smallint) sets
zindex of polygons and it is recursive. We pass the container id (
NULL if there is no container) and the current
zindex as arguments to it.
CREATE OR REPLACE FUNCTION setZIndex(pfid Bigint, zindex Smallint) RETURNS Void AS $proc$ DECLARE feat record; curs refcursor; BEGIN IF pfid IS NULL THEN OPEN curs NO SCROLL FOR SELECT fid FROM polygons AS f WHERE f.pfid IS NULL; ELSE OPEN curs NO SCROLL FOR SELECT fid FROM polygons AS f WHERE f.pfid=setZIndex.pfid; END IF; <<curs_loop>> WHILE TRUE LOOP FETCH NEXT FROM curs INTO feat; EXIT curs_loop WHEN NOT FOUND; PERFORM setZIndex(feat.fid, zindex+1::Smallint); UPDATE polygons SET zindex = setZIndex.zindex WHERE CURRENT OF curs; END LOOP curs_loop; END; $proc$ LANGUAGE plpgsql;
At the very beginning, one of the two applicable cursors is selected depending on the container id value. Inside the loop, the procedure is called again for those polygons whose container is current polygon.
zindex value is increased by one on every recursive call.
Let’s execute this code.
INSERT INTO polygons (fid, geom) VALUES (1, ST_GeomFromText( 'POLYGON((0 0,100 0,100 100,0 100,0 0))',3785)); INSERT INTO polygons (fid, geom) VALUES (2, ST_GeomFromText( 'POLYGON((150 0,200 0,200 50,150 50,150 0))',3785)); INSERT INTO polygons (fid, geom) VALUES (3, ST_GeomFromText( 'POLYGON((20 20,80 20,80 80,20 80,20 20))',3785)); INSERT INTO polygons (fid, geom) VALUES (4, ST_GeomFromText( 'POLYGON((40 40,60 40,60 60,40 60,40 40))',3785)); UPDATE polygons SET area=ST_Area(geom); SELECT arrangePolygons(); SELECT setZIndex(NULL::Bigint,0::Smallint); SELECT fid, pfid, area, zindex FROM polygons;
We should get following results.
pfid column gives us the container id, the value in the column
zindex shows to us the type of the polygon object: if the value is even then it is the land coastline, otherwise it is the lake shore.
Zero value corresponds to the isolated polygons: continents or oceanic islands.