Category Archives: All

Planet surface tessellation with irregular tiles. Part 2

Ocean segments

Actually, the ocean segments are almost the same polygons, about which we spoke earlier, only without the subtracted areas of the continents and small polygons merged with neighbors.

Some of the parts left after lands subtracting can be very small. Also, small polygons can be formed near the border of the ‘world’ polygon. We should add them to the neighboring sectors, so as not to spoil the overall picture.

Red lines in the figure indicate the desired merging of polygons.

desired merging

The choice of rectangles for merging with can be done in different ways. The project implemented two methods: by the longest border and by the nearest center. In the following image, the second method is applied.

ocean sectors

In general, we achieved what we wanted. The area of the sectors is approximately the same (on the sphere) and sectors are mainly convex with small exceptions.

Land segments

On the continents we can create a more interesting subdivision. It based, from one side, on generated small Voronoi polygons (base polygons), and, from the other side, on natural geographical objects such as rivers and watersheds. The obtained structure is very similar to the world’s subdivision into the states and provinces in them. That is, it is the process of the procedural generation of the world political map.

Unlike the ocean, convexity is unnecessary, and we should consider natural boundaries, other than land borders, as possible boundaries between sectors.

Watersheds in the project are not yet build, but there are already rivers!

As the first step, we create the Voronoi diagram of small area polygons for the selected land. In the project, each division of the global landmass is marked with the identifier aid (area Id) and can consist of a continent and nearest islands or only small islands located close by each other. The availability of areas and their identifiers can be determined from the content of the lands shapefile.

The average size of the base polygons is made as small as possible (it increases the computation time!). For the current example, we choose the average sector size of 40,000 km2 and the average size of the base polygons at 5000 km2. This is a rather large size chosen so only for clearness of the example.

In the picture there is a part of Voronoi diagram inside a land.

small Voronoi polygons

In practice, it is better to choose the average size of the base polygons in the interval 500–1000 km2. Smaller sizes are also permissible (up to several m2), but take into account the computation time and use a more suitable database configuration.

To merge base polygons into sectors, we randomly choose from them as much as we want to make separate sectors, and then gradually join remaining adjacent base polygons to resulting sectors.

In the picture an example of what we can get.

raw sectors

Now it is time to consider the rivers. We use the representation of rivers as two-dimensional linestrings, by which we can divide the initial sectors into parts (with ST_Split function). The river linestrings can be taken from riversz shapefile. Since we do not need the z-coordinate stored in that shapefile, we cut it off when importing shapefile to PostGIS enabled database.

Division of sectors by rivers can be performed depending on a number of conditions: the final streamflow of the river, the area of the separated sector part, etc.

After this division we will get many separate parts.

rivers cut sectors

Small parts of sectors (cuts) need to be attached to large sectors. But we must try not to attach parts to the same sectors from which they were cut off.

After joining, we will see the following picture.

land sectors

You might think that we have done everything that is required. But in fact it is not so, because there are still small islands. Those of them that have an area close to the average sector area are made as individual sectors right away. But with those whose area is much smaller than the average sector, we act in two ways.

1. If such islands are located relatively far from existing sectors, then we make them separate sectors. “Farness” here depends on the defined average sector area: the smaller this area the more likely that the island will become a separate sector.

2. If the island is close to other already defined sectors, then it merges with one of them. With the one that has the largest area in some polygon obtained with ST_Buffer of the island polygon.

On the next picture, four islands that have merged into two sectors.
2 sector archipelago

With an increase in the defined average sector size, these islands will fall into single sector.

Planet surface tessellation with irregular tiles. Part 1

Spherical tiling

Usually, the problem of dividing the sphere into tiles (tessellation) is solved seeking to obtain as much as possible a regular structure on the sphere. If there are few tiles, then a spherical polyhedra can be used. 20 is the maximum number of spherical divisions which can be achieved by a spherical icosahedron. And hosohedra is more often far from the desired.

There are derived methods.
1. Further fission of divisions can produce tessellation with a large number of tiles.
2. Using of semiregular polyhedra.
But the tessellation obtained by such methods is not always convenient in practice. It can completely consist of triangles, or it can consist of dissimilar elements. Imagine a strategic game with a tile map made up of triangles!

Now let’s take a look at this from the other side. What if the tiles must match some defined structure on the sphere?

We are dealing with whole planets and we want some tiles to be parts of the ocean and others to be parts of the continents. Usually, when creating tile maps, tiles are divided into the ocean and land ones, according to some principle, thereby determining the map content. And the resulting map has a resolution defining by the the number of tiles in the world map.

But, what if the division into the ocean and continents is already defined on the planet surface, and with a resolution far exceeding possible tile resolution? There may come in handy irregular tiles. Despite the name, irregular tiles may have imposed requirements. For example, there may be a convexity of tiles or an uniformity of size within certain limits.

On the plane such tessellation is, usually, made by using a Voronoi diagram. Here, we will see how to generalize this method to the surface of the sphere. And how the generated Voronoi polygons can be used to make ocean and land sectors.

We use the term ‘sector’ instead of ’tile’, since tiles are generally acting as elementary objects. But inside our sectors there are always an internal structure (terrain and geographic objects). Sectors can also act as elementary objects, serving as a tile map on sphere. For this purpose, a graph of possible transitions between sectors is also build.

SQL script constructing planet surface sectors. See also the tutorial on its use.

Polygons on planet sphere

Here, I give a description of how to partition the planet sphere into polygons using the arsenal of PostGIS tools. The only imposed requirements for polygons: convexity and uniformity of their spherical area with possible deviations of the order of the predetermined average area. Thus obtained polygons will be used to create sectors in the ocean and in the continents.

In the case of a flat map, one could break it into a Voronoi diagram of convex polygons, and could further use Lloyd’s algorithm to bring the diagram to ultimate form (Centroidal Voronoi tessellation). The essence of Lloyd’s algorithm is to repeat the creation of the Voronoi diagram, where at each iteration the centers of polygons obtained in the last iteration are taken as reference points for the new iteration.

In PostGIS there is a function ST_VoronoiPolygons, which builds the Voronoi diagram in a square area. Let’s see how we can use it for our purposes.

What will happen if we will try to use a straightforward approach? Any rectangular projection of the planet can be converted to a square by the coordinate transformation, then diagram is build, and then we can apply reverse coordinate transformation. But, the rectangles constructed will be stretched in one direction, which would be undesirable. And, if we will try to directly use Lloyd’s algorithm, polygons close to the poles of the sphere will have substantially smaller spherical area than those close to the equator.

Therefore, I offer another adapted method, which has not these flaws.

The random reference points for Voronoi diagram are chosen uniformly distributed on the sphere. In a Mercator projection, this means that they should appear less frequently towards the poles with probability proportional to the Cosine of the latitude. We construct the Voronoi diagram in the ‘world’ polygon, which is the entire planet projection or part of the projection of the planet. The ST_VoronoiPolygons function itself completes the rectangle to a square, we just need to subtract added space from the resulting diagram polygons.

The sample of the Voronoi diagram obtained by the proposed method is on the picture. Here is a part of the test planet map from the equator at the top to 73 degrees of south latitude at the bottom edge. And the lands are already subtracted from the Voronoi polygons.

Voronoi polygons

Apparently, in the Mercator projection, the polygons are usually have larger area when approaching the poles, but as spherical polygons they have the same distribution anywhere, this is what we wanted. But at the same time, the polygons have a very large dispersion in area and, because of this, the overall picture is unattractive.

Let’s try to apply the partial Lloyd’s algorithm. As new reference points for the subsequent iterations we select the centers of the polygons obtained by clipping with ‘world’ polygon. And, to prevent too much area reduction of polygons near the poles, we do only a small number of iterations (3 or so).

After applying such a modified algorithm, we obtain the following picture, which can already be acceptable.

Lloyd's algorithm

In the next post we will see how to make ocean sectors and, more interestingly, land sectors.

Change Log: 0.9

  • Possible tiling of the planet sphere with tiles of various polygonal form. Ocean sectors on the planet map (ocean_sectors shapefile).
  • Division of the landmasses into sectors taking into account rivers. Sectors can act as provinces of possible states (land_sectors shapefile). The procedural generation of the political map.
  • Graph of possible transitions between neighboring sectors (ocean_sectors_graph.dbf and land_sectors_graph.dbf files respectively).
  • The height contours of the example planets are moved to a separate archive.
  • Fixed bug in generation of LineStringZ rivers.

The Nishke II and Shkay II planets were updated.

The next few posts and tutorials will be devoted to the spherical tiling theory and practice.

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.


    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.

    RETURNS Void AS $proc$
    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;
    SET enable_seqscan TO off;
    FOR feat IN curs_feats LOOP
        pfidV = NULL;
        IF feat.geom IS NOT NULL THEN
            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;
$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$
    feat record;
    curs refcursor;
    IF pfid IS NULL THEN
        OPEN curs NO SCROLL FOR
            SELECT fid FROM polygons AS f
            WHERE f.pfid IS NULL;
        OPEN curs NO SCROLL FOR
            SELECT fid FROM polygons AS f
            WHERE f.pfid=setZIndex.pfid;
    END IF;
        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;
$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.