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 – Mike T Jul 26 '12 at 21:45
  • 2
    For a truly spherical analysis, a better question could be how to create a geodesic grid – Mike T Jul 26 '12 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))'));


308911036269.806 (1 row)

geodata=# select 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))'));
(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. – kurt-hectic Jul 26 '12 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. – kurt-hectic Jul 26 '12 at 9:05
  • I did a copy/paste of your query into PostGIS, and it worked OK. (shrug...) – Micha Jul 26 '12 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. – Micha Jul 26 '12 at 11:17

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


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

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? – Micha Jul 28 '12 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.