2

I am trying to perform spatial statistics using postgis. Once in a while I have ST_Clip crushes and halt the query. I figure that this occurs when polygon barely intersects with raster. Please see the sample below.

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))')
        )
);

With the above query I am getting following error.

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

I am hoping to getting no record returned instead, or some kind of empty raster. In my production code, the geometry/raster pair was found by ST_Intersects(r.rast, p.geom) between table of polygons and raster. One way I thought about making bounding box for raster which is slightly smaller than the extent of raster, but this is pretty ugly...

My version of postgres and postgis are

  • PostgreSQL 9.6.1 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 4.9.1, 64-bit
  • POSTGIS="2.3.1 r15264" GEOS="3.6.0-CAPI-1.10.0 r0" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.1.2, released 20 16/10/24" LIBXML="2.9.4" LIBJSON="0.12.1" RASTER

Thanks!

yosukesabai
  • 6,184
  • 4
  • 30
  • 42

1 Answers1

1

My tentative solution, wrap with begin/exception/end block and let the exception part to return empty raster. Performance suffers (~ two times). It will create false negative, but not sure what to look for...

-- 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
yosukesabai
  • 6,184
  • 4
  • 30
  • 42