Getting the median pixel value for a raster/polygon intersection in PostGIS

Question

I followed the instructions of this article to obtain zonal statistics for the intersection of a raster and polygon. For context, I'm wanting to obtain zonal statistics for every plot within a field experiment. These are the polygons (.gpkg) overlaid on the raster (.tif) in qGIS:

This is my implementation of the SQL code:

CREATE TABLE geospatial.ortho06_21_2022_sub_stats AS
SELECT
    CAST (geo.id AS INTEGER) AS geom_id,
    geo.geom AS geometry,
    COUNT(1) AS n_tiles,
    (ST_SummaryStatsAgg(ST_Clip(raster.rast, geo.geom, true), 1, true)).*
FROM
    geospatial.ortho06_21_2022_sub_ndvi AS raster
JOIN geospatial.allbuffers0724_23 AS geo
    ON ST_Intersects(geo.geom, raster.rast)
GROUP BY
    geom_id, geometry;

This code seems to produce my desired output, providing summary stats for each of the 1800 polygons, as shown in this table:

geom_id	geometry	n_tiles	count	sum	mean	stddev	min	max
1	0103000020647F00000100000...	8	7251	5921.466987863183	0.8166414270946328	0.16005520481227006	0.1760474294424057	0.9414851069450378
2	0103000020647F00000100000...	4	7245	6093.110233664513	0.8410090039564545	0.10651537882310805	0.20441795885562897	0.9480032324790955
3	0103000020647F00000100000...	4	7248	5920.850203022361	0.8168943436841005	0.156718020606931	0.22542105615139008	1

However, the ST_SummaryStatsAgg() function does not produce the median, and I can't find any PostGIS functions that can deal with this. Other online information seems pretty scarce as well. Does anyone know how to obtain the median in this context?

Answer 1

For each polygon:

• Clip the raster
• Convert the pixels to points
• Extract the raster values from all points
• Then calculate mean value per polygon id with percentile_cont:

SELECT polyid, PERCENTILE_CONT(0.5) WITHIN GROUP(ORDER BY rastervalue asc) AS median
FROM
    (SELECT poly.id as polyid, (st_pixelaspoints(ST_CLIP(raster.rast, poly.geom))).val as rastervalue
    FROM test.grid_nh_72_7 as poly
    JOIN test.nh_72_7_tiled AS raster
    ON ST_INTERSECTS(poly.geom, raster.rast)) sub
GROUP BY polyid