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