# How to create a valid global polygon grid in PostGIS?

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?

• In short, you need the `geometry` type for this particular grid, not `geography` Commented Jul 26, 2012 at 21:45
• For a truly spherical analysis, a better question could be how to create a geodesic grid Commented Jul 26, 2012 at 21:57

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)
``````
• As to version: I use "POSTGIS="2.0.0 r9605" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER", on windows xp. Regarding behaviour. I would want the point to fall into exactly one cell. What I could do is to use degrees, minutes, seconds (which are the reference system of my objects) to generate the cells. 1 sec would be the minimum step in this domain, but I am still interested how this can be solved generically in postgis. Commented Jul 26, 2012 at 8:51
• note that your query had different coordinates from mine, even though the shape is similar. I just ran yours and it works, but mine keeps failing. Commented Jul 26, 2012 at 9:05
• I did a copy/paste of your query into PostGIS, and it worked OK. (shrug...) Commented Jul 26, 2012 at 10:52
• Regarding the PostGIS version, I found a postgis-users maillist entry regarding that error with ST_Area(). See the additional query in the answer above. Commented Jul 26, 2012 at 11:17

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.

• This question of points falling precisely on the boundary piqued my curiosity. I did a quick experiment: created (in Spatialite for convenience) a polygon grid of 5 deg x 5 deg, with a spacing of 0.02 degrees (100 rows x 100 cols = 10,000 cells). THen I made a layer of 1,000,000 randomly located points in the same extent. I then did a query of how many points fell into each cell. Results were: Avg=100 (not surprising), min=68 and max=144. So I guess the spread was OK. Then I did a count of how many points fell on the boundary. Answer = 0. Maybe this is a non-issue? Commented Jul 28, 2012 at 10:31