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;