当几何几乎不相交时，光栅 ST_Clip 失败

Question

我正在尝试使用 postgis 执行空间统计。有时我会遇到 ST_Clip 崩溃并停止查询。我认为当多边形几乎不与栅格相交时会发生这种情况。请参阅下面的示例。

SELECT ST_Summary(
        ST_Clip(
                ST_AddBand(
                        ST_MakeEmptyRaster(16, 16, 0, 0, 1, 1, 0, 0),
                        ARRAY[
                        ROW(1, '8BUI'::text, 0, 255),
                        ROW(2, '8BUI'::text, 0, 255),
                        ROW(3, '8BUI'::text, 0, 255)
                        ]::addbandarg[]
                )
                -- this works
                --, ST_GeomFromText('POLYGON((15.999999 15.999999, 15.999999 17, 17 17, 17 15.999999, 15.999999 15.999999))')
                -- this fails
                , ST_GeomFromText('POLYGON((15.9999999 15.9999999, 15.9999999 17, 17 17, 17 15.9999999, 15.9999999 15.9999999))')
        )
);

通过上述查询，我​​收到以下错误。

psql:demo_clip_fail_barelyintersects.sql:16: ERROR:  RASTER_clip: Could not get band from working raster 
CONTEXT:  PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN

我希望不会返回任何记录，或者某种空栅格。在我的生产代码中，ST_Intersects(r.rast, p.geom) 在多边形表和栅格之间找到了几何/栅格对。我想过为栅格制作边界框的一种方法，它比栅格的范围略小，但这很丑……

我的postgres和postgis版本是

• x86_64-pc-linux-gnu 上的 PostgreSQL 9.6.1，由 gcc (GCC) 4.9.1 编译， 64 位
• POSTGIS="2.3.1 r15264" GEOS="3.6.0-CAPI-1.10.0 r0" PROJ="Rel. 4.9.3，2016 年 8 月 15 日" GDAL="GDAL 2.1.2，于 2024 年 16 月 10 日 20 日发布" LIBXML="2.9.4" LIBJSON="0.12.1" RASTER

谢谢！

Answer 1

我的暂定解决方案，用 begin/exception/end 块包装，让异常部分返回空栅格。性能受到影响（〜两次）。它会产生假阴性，但不确定要寻找什么......

-- function to work around bug in st_clip (fails when polygon barely intersects with raster)
-- not sure how much damage this has on performance
create or replace function st_clip_fuzzy(
        rast raster, nband integer[],
        geom geometry,
        nodataval double precision[] DEFAULT NULL, crop boolean DEFAULT TRUE
)
        returns raster
        as $$
        declare
        rec record;
        g geometry;
        begin
                return st_clip($1, $2, $3, $4, $5);
        exception
        when others then
                select st_intersection(st_envelope(rast), geom) into g;
                raise warning 'st_clip_fuzzy: intersection %', st_astext(g);
                raise warning 'st_clip_fuzzy: area intersection %', st_area(g);
                raise warning 'st_clip_fuzzy: area pixel %', abs(ST_ScaleX(rast) * ST_ScaleY(rast));
                raise warning 'st_clip_fuzzy: area ratio %', st_area(g) / abs(ST_ScaleX(rast) * ST_ScaleY(rast));

                return ST_MakeEmptyRaster(0, 0, ST_UpperLeftX(rast), ST_UpperLeftY(rast), ST_ScaleX(rast), ST_ScaleY(rast), ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast));
        end;
        $$ language 'plpgsql' immutable;

CREATE OR REPLACE FUNCTION st_clip_fuzzy(
        rast raster, nband integer,
        geom geometry,
        nodataval double precision, crop boolean DEFAULT TRUE
)
-- four more interfaces with different set of arguments