All posts by mamont

Change Log: 0.8

  • Height contours in the contours shapefile. For the example planets, the contours were made at intervals of 200 and 250 meters; mlevel parameter specifies the height of a contour. For the present, contours are without smoothing.
  • The hlevel parameter value of the lakes is rounded down.

The Nishke II and Shkay II planets were updated.

The height contours are obtained without the use of Digital elevation model data (DEM), therefore they have their own distinctive properties.

The granularity of the contour data does not depend so much on position of the points. In the case of a traditional approach of getting contours from DEMs, their resolution vary significantly on the latitude. In the project, dependence on latitude and, possibly, on longitude is due to the method of cones to sphere mapping, but this dependence is not as great as in the former case.

The resolution of the contours is not limited by the generally accepted DEM formats, but only by the computational resolution of the planet, that is, by the parameters gn and ln.

If GeoTiffs are not required, but only height contours, then there is no need to produce a resource consuming DEMs generation operation.
 

Boundaries and nesting of land and water areas

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.

a lake island

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:

  • fid is feature id of the polygon.
  • pfid is the id of the polygon container in which the given polygon is directly contained.
  • zindex is the layer number in which a polygon is located.
  • area is 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.
  • geom is the field with geographic data in the Mercator projection (EPSG:3785).

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

area_idx and pfid_idx are the standard indices on area and pfid fields. 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.

The procedure 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.

fid pfid area zindex
2 2500 0
4 3 400 2
3 1 3600 1
1 10000 0

The 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.

Change log: 0.7.1

  • Correct lakes and lands borders. Excluding only small features which size comparable to the map resolution.
  • Fixed invalid river channel heights in riversz shapefile.
  • Water level of the lakes is stored in the new field ‘hlevel’ of lakes shapefile. For lowlakes water level is always 0 (water level of the ocean), so there is no need to store it.
  • Some minor bug fixes.

The Nishke and Shkay planets data are replaced by new ones generated under the same configuration parameters. The planets named Nishke II and Shkay II respectively.

This is an important correction, which gives way to some interesting applications.

 

Change Log: 0.7

  • The mechanism for combining two or more height meshes to obtain a great variety in planet relief forms. In particular, it is now possible to make continents with alternate mountainous and lowland areas.
  • Rivers are now more likely to flow into each other than flow in the same direction to the ocean.
  • River width variations fixed. Earlier, the value of the channel shift depended on the latitude.
  • The way of generating relief inside rhombs is again changed: sharp uplifts near the rivers are removed and there are no ubiquitous bulges.

The two new example planets added to the site.
 

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;