Creating raster with cell value counts from another raster in buffer around each cell using PostGIS

Question

• Main objective

My main objective is to compute 1) the proportion of a given type of land cover around every raster cell in a radius of 500 m at the scale of an island (e.g.: there is 40 % of agricultural areas at a radius of 500 m around this point, if the selected land cover type is agricultural area), 2) the majority land cover type at a similar scale around each cell. (Such values are needed to extrapolate a Species Distribution Model).

• Source layer

The source layer for this calculation is a categorical PostGIS raster layer with 12 different cell values corresponding to land cover types (named 'land_cover') within a schema named 'schema'. It has around 18*10^6 cells, which necessarily makes the calculation computationally intense. Agricultural areas correspond to value 1 in 'schema.land_cover'

• First tries

I have tried several PostGIS queries. The rationale followed so far involves:

• Getting the coordinates of each raster cell in 'land_cover' (ST_PixelAsCentroids(rast,1)) and creating a circular buffer around that point centroid geometry (ST_Buffer(geom, 500)).
• Clipping the initial raster 'land_cover' with that buffer and merging the result to have a unique result value even if the buffer intersects several 'rid' of 'land_cover' (ST_Clip(ST_Union(rast), geom))
• Using ST_ValueCount(rast, 1) on the clipped circular raster obtained
• Selecting the 'count' of values corresponding to agricultural areas (WHERE value = 1)
• Giving this value 'count' and the pixel centroid coordinates as input to ST_SetValues(ST_AddBand(ST_MakeEmptyRaster(rast), '32BF'), 1, ARRAY_AGG(ROW(geom, count)::geomval)::geomval[] ) by taking the initial 'land_cover' raster as a reference for creating the new raster grid.

Applying this while paying attention to GROUP BY conditions leads to the following query, which seems to be correct (checked with 'Explain' in pgAdmin) but it is still running after one complete day, and I need to run several of these queries:

CREATE TABLE schema.nb_agri_500 AS 
 SELECT r.rid, 
 ST_SetValues(ST_AddBand(ST_MakeEmptyRaster(rast), '32BF'), 1, 
 ARRAY_AGG(ROW(geom, count)::geomval)::geomval[] ) rast 
 FROM ( 
 SELECT * FROM ( 
 SELECT rid, geom, (ST_ValueCount(rast, 1)).* 
 FROM ( 
 SELECT b.rid, b.x, b.y, b.geom, ST_CLIP(ST_UNION(r.rast), b.geombuf) AS rast 
 FROM ( 
 SELECT rid, 
 (ST_PixelAsCentroids(rast, 1)).*, 
 ST_Buffer((ST_PixelAsCentroids(rast, 1)).geom, 500) geombuf 
 FROM 
 schema.land_cover ) 
 AS b, 
 schema.land_cover 
 AS r 
 WHERE ST_Intersects(b.geombuf, r.rast) 
 GROUP BY b.rid, b.x, b.y, b.geom, b.geombuf 
 ) AS a 
 ) t 
 WHERE value = 1 
 ) t, 
 schema.land_cover AS r 
 WHERE t.rid = r.rid 
 GROUP BY r.rid, r.rast

Then, for getting the proportion of a given value, I need the total number of cells around each cell. I therefore adapt this query to obtain the total number of cells in a radius around a cell, basically removing the WHERE value = 1 and several SELECT * FROM() and changing ST_ValueCount() for ST_Count().

For getting the majority cell value in a given area, I also use a similar query with an additional operation ranking cell values according to their frequency and retaining only the first line of the resulting table. It leads to the following query, which has GROUP BY troubles:

CREATE TABLE schema.land_cover_max AS 
 SELECT rid, 
 ST_SetValues(ST_AddBand(ST_MakeEmptyRaster(rast), '8BUI'), 1, 
 ARRAY_AGG(ROW(geom, value)::geomval)::geomval[] ) rast 
 FROM ( 
 SELECT * FROM ( 
 SELECT *, rank() OVER (PARTITION BY rid, x, y ORDER BY count DESC) AS rang 
 FROM ( 
 SELECT *, (ST_ValueCount(rast, 1)).* 
 FROM ( 
 SELECT b.rid, b.x, b.y, b.geom, ST_CLIP(ST_UNION(r.rast), b.geombuf) AS rast 
 FROM ( 
 SELECT rid, (ST_PixelAsCentroids(rast, 1)).*, 
 ST_Buffer((ST_PixelAsCentroids(rast, 1)).geom, 500) geombuf 
 FROM schema.land_cover ) 
 AS b, 
 schema.land_cover 
 AS r 
 WHERE ST_Intersects(b.geombuf, r.rast) 
 GROUP BY b.rid, b.x, b.y, b.geom, b.geombuf 
 ) AS a 
 ) t 
 ) t 
 WHERE rang = 1 
 ) t, 
 schema.land_cover AS r 
 WHERE t.rid = r.rid 
 GROUP BY r.rid, r.rast

How do I fasten these queries and make them functional?

I have dug into PostGIS add ons but 1) I was helped for these first tries and am not sure to handle PostGIS sufficiently well to use correctly these functions and 2) did not find functions performing exactly what I am trying to do here.

Answer 1

So far, I have not found more efficient solutions for my calculations on PostGIS. Due to memory limits, even the queries searching for total number of cells of a given value could not be run entirely. A solution would have consisted in splitting the raster extent into several parts, running the queries for these separate parts and merging the results (by taking into account edge effects).

Yet, the function r.neighbors on GRASS (launched from QGIS) exactly does the calculations I was looking for. I created binary rasters with cell values 0 (background) and 1 (agricultural area) and computed the sum of cell values in a circular zone around every raster cell defined by its radius. The proportion could then be obtained by dividing the result by the number of cells in a circular buffer.

The option 'mode' in r.neighbors function makes it possible to get the majority cell value in the circular area around every raster cell.

Also useful for Species Distribution Modelling, the function r.grow.distance computes the distance from each raster cell to raster cells in another raster.

https://grass.osgeo.org/grass78/manuals/r.neighbors.html https://grass.osgeo.org/grass79/manuals/r.grow.distance.html

So, I did not find an optimal solution on PostGIS but GRASS made the job.