ST_Clip, subpixel polygons, and related curiosities
Dian Fay
di at nmfay.com
Sun Jul 14 14:04:03 PDT 2024
When I went digging I found a TODO note "test polygon smaller than
pixel" in the ST_Clip source, so that goes a long way toward explaining
my initial issues. I started with 3.4 and upgraded to 3.5-dev (commit
2860cdd1, built from source) when I saw ST_Clip got the `touched`
argument. The new function signature was available immediately but
returned the same results until I came back after restarting Postgres
this morning, at which point it did exactly what I hoped for. There are
some other interesting things going on that I thought would be worth
writing up.
Here's a script that clips an all-999 raster to three numbered,
successively smaller triangles:
```
with r as (
select st_addband(
st_makeemptyraster(300, 251, 97.4999889, 52.716667838, .008333333, -.00833333, 0, 0, 4326),
1, '16BSI', 999, -32768
) as rast
), g as (
select 1 as rn, st_setsrid(st_geomfromtext('POLYGON((98 52,98 52.1,98.01 52,98 52))'), 4326) as geom
union
select 2 as rn, st_setsrid(st_geomfromtext('POLYGON((98 52,98 52.1,98.005 52,98 52))'), 4326) as geom
union
select 3 as rn, st_setsrid(st_geomfromtext('POLYGON((98 52,98 52.1,98.004 52,98 52))'), 4326) as geom
), c as (
select rast, geom, rn, st_clip(rast, 1, geom, -32768) as clipped
-- try using the `touched` arg in 3.5.0; if just upgrading from 3.4,
-- results are the same until the server restarts
-- select rast, geom, rn, st_clip(rast, 1, geom, -32768, true, true) as clipped
from r, g
)
select
rn,
st_summary(clipped) as clipped_summary,
(st_summarystats(clipped)).*,
st_value(clipped, st_centroid(geom)) as clipped_centroid_value
from c
order by rn;
```
Things I understand about the results:
- All clipped rasters have the same dimensions. The minimum width is 2px
whether the clip triangle's x-length is the high bound of 1/100° or
the low bound 1/250°.
- Somewhere between 1/200° and 1/250°, the length is too short for
ST_Clip with `touched=false` to populate even one pixel's value --
there is a point (if you will) at which it rounds down to nothingness,
indicated by ST_SummaryStats' `count` reaching 0.
- `clipped_centroid_value` becomes null _before_ ST_SummaryStats reaches
zero-count. My guess is that the data pixels are aligned vertically
(stacked) since the triangles are taller than they are long, and so
the centroid of triangle 2 falls just barely in the second column
which is all NODATA.
Something I don't yet understand: ST_SummaryStats prints the "All pixels
of band have the NODATA value" NOTICE message one time per row six
times, once for each attribute of the `summarystats` type. If you select
`st_summarystats(clipped)` as a tuple without the parentheses-asterisk
decomposition, it prints once. I found the message in
rt_band_get_summary_stats, but I'm not sure why that would execute
multiple times in the one case.
And I found a really odd behavior change while working on a fallback
strategy before I found the new `touched` arg. Add the following lines
to the final select list with ST_Clip `touched=false`:
```
-- st_astext((st_intersection(clipped, geom)).geom) as tile_intersection,
-- (st_valuecount(clipped)).*,
```
Uncomment ST_ValueCount. There's a new NOTICE, "Cannot count the values
for the band at index 1", and row 3 -- where the clipped raster is fully
NODATA -- is gone, with only two records returned.
Now uncomment the ST_Intersection. It still prints "Cannot count", _but_
it returns the third row! Both `value` and `count` are null as might be
expected. Intersecting the original raster also preserves the third row.
Is ST_ValueCount silently dropping all-NODATA rows? And why does an
unrelated invocation of ST_Intersection bring them back?
More information about the postgis-users
mailing list