Edit:
Per a comment, here is some more context and the basics of what I'm looking to do, that I think would better enable someone to offer input from a first principles perspective. Its a lot of detail because my ideal can seem a bit off from what my original post was about, so wanted to show how I got there, and ultimately I'd be happy to either get to my ideal results, or the adaptation I describe and that my original post was about.
My ideal analysis would actually be a "moving window analysis" from the raster landcover dataset I mentioned in my original post. The query I shared is all adapted to be based on a vector data model because I wasn't seeing how I could do what I want with a raster in PostGIS, and my attempts in a couple of other tools weren't fruitful. But I made progress with the below so figured I'd see how I could improve it.
So for my ideal moving window analysis the output would be a new raster with more the same resolution and dimensions of the original, in which the value of every pixel would represent:
- the proportion of land area within a circular area around that pixel with the land cover delineated as either one of two classes (classes "1" or "2" from the 8-class landcover dataset).
The original raster is 6 inch resolution, and the circular area around each pixel would have a radius of 656 feet so each window would have a large number of pixels (~5.5M). A decision would need to be made around how values around edges of the raster are calculated - either the resultant raster would not have any results values (i.e., it would have no data values) within 656 feet of the edge, or it would have values calculated for pixels up to the edge based on the data that exist (so the moving window for those outer areas would be smaller and not actually circular); it would be reasonable to go either direction on that. Likewise, whether the results are calculated based on pixel centroids that intersect the radius of interest or a calculation based on cutting the pixels would be a decision point, but either way would ultimately be okay. (I believe st_intersection would cut the pixels for the calculations, which would be fine for this.)
I don't really need the results for every 6" pixel, so I was adapting this to calculate the results for every point within a coarser grid to see if I could make progress and get something that is in line with what I need, even if its not my ideal, above. In this adaptation, per the query I shared in the original post, I was looking to calculate:
- For every point in the grid, what is the proportion of area that is delineated as class 1 or 2 in the landcover dataset within a 656' radius.
- Note 1 - in the query I shared, it was convenient, and simple/low-cost to get the proportion of the area that was either of the two focal land cover classes as well as both, so I was including that in the query.
- Note 2 - In some past work involving spatial intersections with the landcover data, I found things ran faster with the vectorized version of the landcover data, so my assumption was that would be the case here, though for something like this, its also more intuitive for me in terms of coding up the SQL.
Original Post:
I'm running a query, where, for a lot of points, I'm looking to calculate % or proportion of specific landcover types within a circular buffer. The landcover is vectorized polygon data from a raster. Both the landcover polygons and point data have gist style spatial indexes built for them.
I've run this for a sample area of interest, with ~800,000 points, and and it takes up around 90 hours. For how many points this is, it maybe isn't too bad, but am wondering if there's a way to improve performance, as I'd like to do this for a much larger area and number of points. The query, below.
The query includes a lateral join, which in my understanding doesn't leverage parallel processing at all (and doesn't appear to from looking at processor usage). Thus, I've tried to incorporate a join with another polygon layer (some simple administrative boundaries) that has multiple units in the focal area, such that I could include the join and group by outside of the lateral join and get some parallel processing going, but no success.
I'm also open to doing this type of thing as a moving window raster analysis, if there would be some performance boosts, but the best way to do that in PostGIS wasn't very clear to me, and attempts in R (using terra) have not been successful for sample areas with the size of the area I'm looking at in the window.
Below are the query, PostGIS version info, and the Query Plan. I think one of the aggregate steps is the main bottleneck.
Here's the query
explain (analyze, buffers, verbose, format json)
create table exploration.pct_veg_bk_32_34_localdata as
select
su.pgid,
su.geom_2263,
(sum(cliparea) filter (where lc_value=1)/sum(cliparea) filter (where lc_value is not null)) as pct_canopy,
(sum(cliparea) filter (where lc_value=2)/sum(cliparea) filter (where lc_value is not null)) as pct_grassshrub,
(sum(cliparea) filter (where lc_value in (1,2))/sum(cliparea) filter (where lc_value is not null)) as _veg
from exploration.points_bk32_34_nobldg as su
cross join lateral
(select lc_value, sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, 656.168)))) as cliparea
from exploration.landcover6in_polygons_bk32_34 as a where st_dwithin(su.geom_2263, a.geom_2263, 656.168) group by a.lc_value) foo
group by su.pgid,
su.geom_2263;
PostGIS Version Info (running on Windows 10 x64)
POSTGIS="3.1.1 3.1.1" [EXTENSION] PGSQL="130" GEOS="3.9.0-CAPI-1.14.1" SFCGAL="1.3.8" PROJ="7.1.1" GDAL="GDAL 3.2.1, released 2020/12/29" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER
Here's the executed query plan (pardon the format - output it as a json to see if pgmustard would offer any tips but there was nothing substantial that I could see)
[
{
"Plan": {
"Node Type": "Aggregate",
"Strategy": "Sorted",
"Partial Mode": "Simple",
"Parallel Aware": false,
"Startup Cost": 469.43,
"Total Cost": 1000851265.44,
"Plan Rows": 786877,
"Plan Width": 60,
"Actual Startup Time": 802.957,
"Actual Total Time": 343221227.984,
"Actual Rows": 787860,
"Actual Loops": 1,
"Output": ["su.pgid", "su.geom_2263", "(sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value = '1'::double precision)) / sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value IS NOT NULL)))", "(sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value = '2'::double precision)) / sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value IS NOT NULL)))", "(sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value = ANY ('{1,2}'::double precision[]))) / sum((sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))) FILTER (WHERE (a.lc_value IS NOT NULL)))"],
"Group Key": ["su.pgid"],
"Shared Hit Blocks": 345877949,
"Shared Read Blocks": 92926,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0,
"Plans": [
{
"Node Type": "Nested Loop",
"Parent Relationship": "Outer",
"Parallel Aware": false,
"Join Type": "Inner",
"Startup Cost": 469.43,
"Total Cost": 1000699791.62,
"Plan Rows": 5508139,
"Plan Width": 52,
"Actual Startup Time": 529.355,
"Actual Total Time": 343207780.279,
"Actual Rows": 4547009,
"Actual Loops": 1,
"Output": ["su.pgid", "su.geom_2263", "(sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision))))", "a.lc_value"],
"Inner Unique": false,
"Shared Hit Blocks": 345877949,
"Shared Read Blocks": 92926,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0,
"Plans": [
{
"Node Type": "Index Scan",
"Parent Relationship": "Outer",
"Parallel Aware": false,
"Scan Direction": "Forward",
"Index Name": "points_bk32_34_nobldg_pk",
"Relation Name": "points_bk32_34_nobldg",
"Schema": "exploration",
"Alias": "su",
"Startup Cost": 0.42,
"Total Cost": 315520.82,
"Plan Rows": 786877,
"Plan Width": 36,
"Actual Startup Time": 24.813,
"Actual Total Time": 231611.089,
"Actual Rows": 787860,
"Actual Loops": 1,
"Output": ["su.pgid", "su.geom_2263", "su.geom_2263_buff656"],
"Shared Hit Blocks": 632162,
"Shared Read Blocks": 70054,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0
},
{
"Node Type": "Aggregate",
"Strategy": "Sorted",
"Partial Mode": "Simple",
"Parent Relationship": "Inner",
"Parallel Aware": false,
"Startup Cost": 469.00,
"Total Cost": 1271.19,
"Plan Rows": 7,
"Plan Width": 16,
"Actual Startup Time": 65.122,
"Actual Total Time": 435.317,
"Actual Rows": 6,
"Actual Loops": 787860,
"Output": ["a.lc_value", "sum(st_area(st_intersection(a.geom_2263, st_buffer(su.geom_2263, '656.168'::double precision, ''::text), '-1'::double precision)))"],
"Group Key": ["a.lc_value"],
"Shared Hit Blocks": 345245787,
"Shared Read Blocks": 22872,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0,
"Plans": [
{
"Node Type": "Sort",
"Parent Relationship": "Outer",
"Parallel Aware": false,
"Startup Cost": 469.00,
"Total Cost": 469.05,
"Plan Rows": 16,
"Plan Width": 1829,
"Actual Startup Time": 7.866,
"Actual Total Time": 8.293,
"Actual Rows": 1890,
"Actual Loops": 787860,
"Output": ["a.lc_value", "a.geom_2263"],
"Sort Key": ["a.lc_value"],
"Sort Method": "quicksort",
"Sort Space Used": 2615,
"Sort Space Type": "Memory",
"Shared Hit Blocks": 345245787,
"Shared Read Blocks": 22872,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0,
"Plans": [
{
"Node Type": "Index Scan",
"Parent Relationship": "Outer",
"Parallel Aware": false,
"Scan Direction": "NoMovement",
"Index Name": "landcover6in_polygons_bk32_43_geom_idx",
"Relation Name": "landcover6in_polygons_bk32_34",
"Schema": "exploration",
"Alias": "a",
"Startup Cost": 0.41,
"Total Cost": 468.69,
"Plan Rows": 16,
"Plan Width": 1829,
"Actual Startup Time": 0.088,
"Actual Total Time": 6.735,
"Actual Rows": 1890,
"Actual Loops": 787860,
"Output": ["a.lc_value", "a.geom_2263"],
"Index Cond": "(a.geom_2263 && st_expand(su.geom_2263, '656.168'::double precision))",
"Rows Removed by Index Recheck": 0,
"Filter": "st_dwithin(su.geom_2263, a.geom_2263, '656.168'::double precision)",
"Rows Removed by Filter": 504,
"Shared Hit Blocks": 345245784,
"Shared Read Blocks": 22872,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0
}
]
}
]
}
]
}
]
},
"Planning": {
"Shared Hit Blocks": 199,
"Shared Read Blocks": 51,
"Shared Dirtied Blocks": 0,
"Shared Written Blocks": 0,
"Local Hit Blocks": 0,
"Local Read Blocks": 0,
"Local Dirtied Blocks": 0,
"Local Written Blocks": 0,
"Temp Read Blocks": 0,
"Temp Written Blocks": 0
},
"Planning Time": 158.742,
"Triggers": [
],
"Execution Time": 343230478.411
}
]
