What is the best way of converting a value from irregularly shaped poylgons to a polygon grid using PostGIS? I have a set of polygons that look like this:

enter image description here

and I would like to overlay a grid on these polygons, and assign values from the polygon to the grid value, proportional to the area of the polygon that intersects the grid. For example, where a grid intersects more than one polygon, I want to use the fraction of each polygon area that intersects with the gridcell to proportionally calculate the value.

enter image description here

My approach for doing this is as follows. Assuming each polygon has a value X_i and an area A_i, I calculate the density using D_i = X_i / A_i. Then, for each gridcell, I multiply the fraction of the polygon area that intersects that gridcell by the appropriate density value. I think this approach is satisfactory but I have not figured out a way of abstracting it so that it is easier for me to do this regularly.

I am curious to learn if there is a better way to do this, if there is a term for this process and if there is a function that exists to perform a calculation of this type.

  • I had a similar question gis.stackexchange.com/questions/31108/… which I solved using whuber's suggestion
    – Mike T
    Jan 10 '13 at 18:54
  • hmm, interesting. I've done this before using rasters, but I feel like the vector-intersection approach is more accurate. I feel like a vector 'gridfying' function could be a worthy endeavor.
    – djq
    Jan 10 '13 at 19:42

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.

  • That works perfectly, thanks. I grouped by gid and geom. I'm trying to understand the logic behind the LEFT JOIN poly p ON ST_Intersects(g.geom, p.geom) though. Is it like a true/false condition for relating the geometries?
    – djq
    Jan 10 '13 at 22:24
  • 1
    yup, it finds each polygon (1 or more) that intersect with each grid cell
    – Mike T
    Jan 10 '13 at 22:25
  • I forgot to say, it does return this warning message when running NOTICE: geometry_gist_joinsel called with incorrect join type. I can't tell if this is a significant or not.
    – djq
    Jan 10 '13 at 22:46
  • you can probably [safely] ignore it. It happens when there are no overlaps on the grid. If an INNER JOIN is used, the notice should go away, but you will see fewer rows (where there are no overlaps on the grid).
    – Mike T
    Jan 10 '13 at 23:47

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.