How to create a valid global polygon grid in PostGIS?

Question

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?

Answer 1

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)

Answer 2

Going back to the issue of points exactly on the boundary between two grid cells: PostGIS has a function

 `ST_ContainsProperly()` 

which finds points that intersect the interior of a polygon, but not the boundary. This will mean that those points exactly on a boundary will not be counted at all.

UPDATE cells SET point_count=(
  SELECT count(p.id) 
  FROM points AS p, cells AS c
  WHERE ST_ContainsProperly(c.geom, p.geom) AND
  c.id=cells.id; 

Next, you could possibly find points that intersect the boundary of each grid cell (those that got skipped by ST_ContainsProperly) and add 1/2 of them to the count. So points on the boundary would be weighted 1/2 for each of the adjacent grid cells.

UPDATE cells SET point_count=(
  SELECT COUNT(p.id)/2
  FROM points AS p, cells AS c
  WHERE ST_Intersects(ST_Boundary(c.geom), p.geom) AND
    c.id=cells.id) + point_count;

Based on something like that, you could update each cells point count, adding 1/2 of the boundary points.