6

I am looking to speed up a postgis query, regarding zonal statistics. I am using postgis 3.0 with postgreSQL 12 and I was hoping that I could apply the below script to speed up processing using multiple CPUs;

// Parallel processing has no effect 
SET max_parallel_workers = 8;
SET max_parallel_workers_per_gather = 4;

// Zonal statistics
DROP TABLE IF EXISTS schema.zonal_stats;
CREATE TABLE schema.zonal_stats AS 
SELECT id,geom,
(ST_SummaryStats(ST_Union(ST_Clip(raster,geom)))).sum 
FROM schema.raster, schema.ply
GROUP BY id,geom; 

I could use the rasterstats and multiprocessing packages in Python, which does provide a computation speedup. However, I would like to apply parallel processing natively in postgis, as my data is stored in a postgreSQL db.

  • 1
    you should be using a join condition where they intersect – ziggy Feb 26 at 17:56
  • @ziggy , can you provide an example? – tastatham Feb 26 at 17:59
  • can you add the explain analyse output to the question – Ian Turton Feb 27 at 8:28
3

i would add a intersects join clause on this query like this

DROP TABLE IF EXISTS schema.zonal_stats;
CREATE TABLE schema.zonal_stats AS 
SELECT id,geom,
(ST_SummaryStats(ST_Union(ST_Clip(raster,geom)))).sum 
FROM schema.raster join schema.ply on st_intersects(raster,geom)
GROUP BY id,geom; 

this will run more efficiently because it will use an index and only clip the raster sections that intersect the vector.

In terms of the parrallel processing I am not well versed in what happens under the hood but it seems like you followed what the docs and other tutorials have said about setting the parralel workers etc..

Based on your comments that you pre-clipped the raster to the vector I would not even use the geometry in the query

I would try

SELECT id,rast,
    (ST_SummaryStats(rast)).sum 
    FROM schema.raster

or if you are hell bent on the st_union

SELECT id,rast,
    (ST_SummaryStats(st_union(rast))).sum 
    FROM schema.raster
group by id,rast
| improve this answer | |
  • thanks, that is certainly useful (I did clip my raster to the geom previously) but doesn't answer my question. – tastatham Feb 26 at 18:08
  • if you clipped the raster to the geom already then why would you need to do it again in the query you posted? – ziggy Feb 26 at 18:09
  • Sorry, misread that (confusion between ST_Clip and ST_intersects). Anyway, I clipped the raster previously but using the bounding box of the geom as a rough guide to decrease the raster size, when loading the data into postgreSQL. – tastatham Feb 26 at 18:18
  • 1
    okay but what I am saying is you might not need to use st_clip or even st_union in this process. why not just try putting the raster straight into ST_SummaryStats – ziggy Feb 26 at 18:22
  • 1
    I tiled the raster (-t 100x100), so the union is to merge the rasters together. I'm not too familiar with tiling and postgis to be honest. The union makes sense to me when using tiles, is that correct? I am using the GPW population database and applying this to Argentina. – tastatham Feb 26 at 18:40
1

This doesn't answer the question 100% but does answer my problem, which is to speed up zonal statistics.

I am using the following exactextract command line tool for fast zonal statistics. This also takes into account partial overlaps, by weighting the zonal statistic for each raster value by its coverage fraction. Many zonal statistic tools do not account for this.

If someone can produce an example of multithreaded zonal statistics in postgis, I will substitute this as the marked answer with another.

| improve this answer | |

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.