[GIS] How to create a valid global polygon grid in PostGIS

postgisvector-grid

I want to create a grid spawning the whole globe. The grid is used to measure a density value by counting the number of certain spatial objects in each grid cell.

To create the grid, I have created multiline polygons, both using a script of my own, and the function suggested here: How to create a regular polygon grid in PostGIS? . Both algorithms loop from -90 .. 90 , -180 .. 180 in a fixed step size creating the multiline polygons.

I have two problems.

First, adjacent cells share a separating line, making objects on it counting towards two or more cells. I could subtract a very small value from corners, but this risks loosing objects when they come to lay exactly on the original line. The question is thus if there is a constant for the absolute minimum step value that I could use to guarantee that one point exactly lays next to another with no space inbetween.

Second, it seems the algorithm generates invalid geometries. e.g when I compute the square meters of a grid (for the density)

SELECT ST_Area(ST_GeographyFromText('POLYGON((-15 0,-15 5,-10 5,-10 0,-15 0))'))

throws an error ptarray_area_spheroid: cannot handle ptarray that crosses equator
If I avoid the geography type, and use a geometry there is no global projection that outputs the correct size of all grids. e.g

Select ST_Transform( ST_GeomFromEWKT('SRID=4326;POINT(-180 90 0)'), 32630 )

throws couldn't project point (-180 90 0): latitude or longitude exceeded limits (-14)

How could I generate a grid of valid geometries that covers the whole globe?

Best Answer

As for your first problem, what would you want to happen with a point that falls exactly on the boundary between two grid cells?

As for the area calculation:

geodata=# SELECT ST_Area(ST_GeogFromText('POLYGON((-10 0,-10 5,-5 5,-5 0,-10 0))'));

st_area

308911036269.806 (1 row)

geodata=# select postgis_full_version();
                                                           postgis_full_version                                                            
-------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="1.5.3" GEOS="3.3.2-CAPI-1.7.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.6" USE_STATS (procs from 1.5 r5385 need upgrade)
(1 row)

Seems to work for me. What version are you using??

Searching in Google I found this maillist exchange regarding AT_Area and polys touching the equator: postgis-users

Perhaps, to better understand, could you run the same query but crossing the equator. Such as:

geodata=# SELECT ST_Area(ST_GeographyFromText('POLYGON((-2 -1,-2 1,-1 1,-1 -1,-2 -1))'));
     st_area      
------------------
 24728063597.0354
(1 row)

Related Solutions

[GIS] Creating regular polygon grid in PostGIS

Here is a set returning function ST_CreateFishnet that creates a 2D grid of polygon geometries:

CREATE OR REPLACE FUNCTION ST_CreateFishnet(
        nrow integer, ncol integer,
        xsize float8, ysize float8,
        x0 float8 DEFAULT 0, y0 float8 DEFAULT 0,
        OUT "row" integer, OUT col integer,
        OUT geom geometry)
    RETURNS SETOF record AS
$$
SELECT i + 1 AS row, j + 1 AS col, ST_Translate(cell, j * $3 + $5, i * $4 + $6) AS geom
FROM generate_series(0, $1 - 1) AS i,
     generate_series(0, $2 - 1) AS j,
(
SELECT ('POLYGON((0 0, 0 '||$4||', '||$3||' '||$4||', '||$3||' 0,0 0))')::geometry AS cell
) AS foo;
$$ LANGUAGE sql IMMUTABLE STRICT;

where nrow and ncol are the number of rows and columns, xsize and ysize are the lengths of the cell size, and optional x0 and y0 are coordinates for the bottom-left corner.

The result is row and col numbers, starting from 1 at the bottom-left corner, and geom rectangular polygons for each cell. So for example:

SELECT *
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;
 row | col |         geom
-----+-----+--------------------------------
   1 |   1 | 0103000000010000000500000000...
   2 |   1 | 0103000000010000000500000000...
   3 |   1 | 0103000000010000000500000000...
   4 |   1 | 0103000000010000000500000000...
   1 |   2 | 0103000000010000000500000000...
   2 |   2 | 0103000000010000000500000000...
   ...
   3 |   6 | 0103000000010000000500000000...
   4 |   6 | 0103000000010000000500000000...
(24 rows)

Or to make a single geometry collection for the full grid:

SELECT ST_Collect(cells.geom)
FROM ST_CreateFishnet(4, 6, 10, 10) AS cells;

4x6 grid

You can add the x0 / y0 origin offsets (these defaulted to zero).

[GIS] Converting polygon values to a vector grid using PostGIS

There are two approaches: vector-based and raster-based. The former is more accurate, and the later is generally faster but much more complicated.

Vector approach

Use the fraction of area overlap between each polygon and grid cell to weight the values. Say you have a grid (e.g. created this way), and an irregular poly with values in val. To determine val on grid:

SELECT g.gid, g.row, g.col, count(distinct p.gid) as num_poly,
  sum(ST_Area(ST_Intersection(g.geom, p.geom)) / ST_Area(g.geom) * p.val) AS grid_val
FROM grid g
LEFT JOIN poly p ON ST_Intersects(g.geom, p.geom)
GROUP BY g.gid, g.row, g.col;

This query will take a while there are many tens of thousands of grid cells (i.e. 100 × 100, or more).

Raster approach

The idea behind this method is if you rasterise values from an irregular polygon, you get the values. However, the rasterisation methods are limited, as they don't consider the area fraction of polygons within each pixel. So you will need to modify this approach by rasterising at a finer resolution than desired, then generalise the refined raster to the desired resolution. The accuracy of this method increases with a finer initial raster resolution.

I can add more details to this method, but it is much more complicated. I did this approach first with gdal_rasterize with files in a shell, then I moved the logic to GDAL/Python using gdal.RasterizeLayer.