In the practice of the procedural generation of maps, I found it convenient to use the facilities introduced by PostGIS extender for PostgreSQL. One of the problems that arises, the creation of river polygons, has a fairly simple solution by means of PostGIS.
Here, I will give almost all code that is needed for making river polygons. Suppose, that we already have somehow obtained the rivers line segments of different streamflow (discharge) value. All river lines consist of such segments.
Also, we can assume that the segments with the same streamflow can be one or more for any river. More than one segments are needed to create variations in the width of the river channel (see the previous post).
If the database is new, then do not forget to add the PostGIS extension.
CREATE EXTENSION IF NOT EXISTS postgis;
There are two tables in the database. One contains river segments and other contains resulting river polygons.
CREATE TABLE riversegs ( rid bigint NOT NULL, streamflow integer NOT NULL, geom_wgs84 geometry(LINESTRING, 4326) NOT NULL ); CREATE INDEX riversegs_idx ON riversegs (rid); CREATE TABLE rivers ( rid bigint PRIMARY KEY, geom geometry(MULTIPOLYGON, 3785) ); CREATE INDEX rivers_geom_idx ON rivers USING GIST (geom);
We specify different projections for fields of geometry type. EPSG:4326 or WGS 84 (World Geodetic System) projection is convenient for the initial creation of geographic objects during procedural generation, because its unit of measurement is usual degrees of longitude and latitude. EPSG:3785 (Web Mercator) projection is more convenient when it is necessary to make various transformations of geometries, and its unit of measurement is meter. Since the projections are different, at the certain stage we will have to make a transformation from one projection to another.
This code also includes the creation of two indices, the second of which is the spatial index. Such indices are a good tool for speed up many operations with geospatial data; see details in the PostGIS documentation. In this example, rivers_geom_idx
index is not particularly useful, but if you intend to use the resulting rivers data further, then including it more than desirable. On the other hand, we do not create any spatial index with the geom_wgs84
field, since we do not need bounding box property of its values in our procedures.
The river line segments are of LINESTRING
spatial type, and resulting rivers are of MULTIPOLYGON
spatial type. We do not specify the POLYGON
spatial type in favor of the more general type, since the rivers can consist of separate pieces, for example, if a lake divides the river into parts.
For contrived, but simple, example we write down the data of the river consisting of three segments.
INSERT INTO riversegs VALUES (1,1,ST_GeomFromText('LINESTRING(1 0,2 0,2 1,3 1)',4326)); INSERT INTO riversegs VALUES (1,2,ST_GeomFromText('LINESTRING(3 1,3 2,3 2,3 4)',4326)); INSERT INTO riversegs VALUES (1,3,ST_GeomFromText('LINESTRING(3 4,4 4,5 4,6 4)',4326));
Below is the main procedure in which we collect all possible rivers and write them into the rivers
table.
CREATE OR REPLACE FUNCTION makerivers( halfw double precision ) RETURNS void AS $BODY$ DECLARE curs_rivers CURSOR FOR SELECT rid FROM riversegs GROUP BY rid; target_geom Geometry; BEGIN <<river_loop>> FOR river IN curs_rivers LOOP SELECT INTO target_geom ( SELECT ST_Multi(ST_Union(segs.geom)) FROM (SELECT makeBuffer(geom_wgs84, halfw, sqrt(streamflow)) as geom FROM riversegs WHERE rid=river.rid ) AS segs); target_geom = validateMultiPolygon(target_geom); INSERT INTO rivers VALUES (river.rid, target_geom); END LOOP river_loop; END; $BODY$ LANGUAGE plpgsql;
The procedure takes one argument, the half value of the minimum river width in meters. Further, in the cycle for all possible rivers, we select all segments of each river and make the combined geometry with ST_Union
function. We use the variant of ST_Union
function which accepts a set of spatial objects as a single argument. With ST_Multi
function we ensure that the geometry will be returned as a MULTI* geometry.
makeBuffer
function is defined as following.
CREATE OR REPLACE FUNCTION makeBuffer( geom Geometry, halfw double precision, sqrt_sf double precision ) RETURNS Geometry AS $BODY$ DECLARE buffed_geom Geometry; w double precision; BEGIN w = sqrt_sf * halfw; buffed_geom = ST_Buffer( ST_Transform(geom, 3785), w, 'endcap=round join=round'); RETURN buffed_geom; END; $BODY$ LANGUAGE plpgsql;
In this procedure we turn the river line segment into the ‘thick’ segment, which is a polygon. The function ST_Buffer
helps us in this. 'endcap=round join=round'
argument forces the function to use this rounding method for line joinings and endings. As the buffer width we take the product of half of the minimum width and the square root of the river current streamflow.
In accordance with the project’s agreement of semi-circular profile of the river bed, we take that the river width is proportional to the square root of the current streamflow. For the sake of justice, it should be noted that there are other ways of calculating the width of rivers depending on the streamflow.
Also, note how the function ST_Transform
transforms our original geometry projection into another projection, in which the unit of measurement is meter.
The validateMultiPolygon
procedure does a validation check and, if necessary, changes the geometry type to MULTIPOLYGON
. If resulting geometry is invalid the procedure returns NULL
. We try to get a valid object two times, firstly by ST_MakeValid
procedure, and then by ST_Buffer
procedure, which also is used to obtain valid geometries.
CREATE OR REPLACE FUNCTION validateMultiPolygon(geom Geometry) RETURNS Geometry AS $BODY$ DECLARE new_geom Geometry; BEGIN IF ST_IsValid(geom)=false THEN new_geom = findSubGeom(ST_MakeValid(geom)); IF (new_geom IS NULL) OR ST_IsEmpty(new_geom)=true THEN BEGIN new_geom = findSubGeom(ST_Buffer(geom, 1)); EXCEPTION WHEN SQLSTATE 'XX000' THEN new_geom = NULL; END; END IF; ELSE IF ST_GeometryType(geom)<>'ST_MultiPolygon' THEN new_geom = findSubGeom(geom); ELSE new_geom = geom; END IF; END IF; RETURN new_geom; END; $BODY$ LANGUAGE plpgsql;
The last procedure tries to make MULTIPOLYGON
geometry from input geometry; returns NULL
if this can not be done.
CREATE OR REPLACE FUNCTION findSubGeom(org_geom Geometry) RETURNS Geometry AS $BODY$ DECLARE new_geom Geometry; BEGIN new_geom = NULL; SELECT ST_Multi(ST_Collect(a.gdump)) as b INTO new_geom FROM (SELECT (ST_Dump(org_geom)).geom AS gdump) AS a WHERE ST_GeometryType(a.gdump)='ST_Polygon'; RETURN new_geom; END; $BODY$ LANGUAGE plpgsql;
That’s all. As you can see, with four small procedures we were able to perform a difficult job. To start the process with minimum river width 100
, just run the command
SELECT makerivers(50);
To view resulting data
SELECT ST_AsText(geom) FROM rivers;