Making river polygons with PostGIS

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;

Forged Maps

Procedural Planet Map Generation

Making river polygons with PostGIS