当几何图形几乎不相交时，光栅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=“版本4.9.32016年8月15日”GDAL=“GDAL 2.1.2，发布日期20 16/10/24”LIBXML=“2.9.4”LIBJSON=“0.12.1”RASTER

谢谢！

Answer 1

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

-- 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