- 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.
