catch errors of postgis st_clip, and fix "Could not get band from working raster"

Question

I have millions of polygons, and one larger table of raster (tiled). I want to st_clip(raster, polygon), then apply st_summarystatsagg() on the clipped raster. Code looks like below.

with pieces as (
select d.pid,
       st_clip(r.rast, d.geom) as rast
       from tbl_raster as r
       inner join tbl_poly as d
       on st_intersects(r.rast, d.geom)
)
insert into tbl_out(pid, val1, val2, val3)
select pid,
       (stats1).mean as v1,
       (stats2).mean as v2,
       (stats3).mean as v3
from (
     select p.pid,
            st_summarystatsagg(p.rast, 1, true) as stats1,
            st_summarystatsagg(p.rast, 2, true) as stats2,
            st_summarystatsagg(p.rast, 3, true) as stats3
     from pieces as p
     group by p.pid
);

The code works for small subset of polygons (thousands), but when I pass millions, it comes back with error

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 want to get more diagnostic message of what polygon causes this trouble. Is there any way to do try-catch block type of thing then I am going to dump some info from there. Or let it give some null value or something and make the query to complete, then I look at output for the problem.

Version of postgresql is 9.6.1, postgis is 2.3.1, gdal 2.1.2, and running this on Linux using psql -f.

Thanks!

EDIT:

Following @Eelke's suggestion, I used for loop in PL/pgSQL and identified where the problem is.

It seems that if polygon touches raster with no shared interior points, ST_Clip fails, if there are multiple band calculation requested. As test code below shows, it seems to work as expected (by me) if I calculate union of raster that ST_Intersects with polygons. But this turned out to be costly. I am hoping to find alternative to ST_Intersects, which doesn't pick things up if no interior points are shared. But I couldn't find a function to do this and couldn't figure out expression... Is there such expression?

Below is an reproducible results I came up with :

--DROP SCHEMA IF EXISTS test_rast;
--CREATE SCHEMA test_rast;

SET search_path TO test_rast, public;

-- Create raster catalog of 3 bands
DROP TABLE IF EXISTS dummy_rast_3band;
CREATE TABLE dummy_rast_3band(rid int, rast raster);
INSERT INTO dummy_rast_3band(rid, rast)
VALUES
(1, ST_AddBand(
                ST_MakeEmptyRaster(240, 240, 82,   59.6, 6./3600., 6./3600., 0, 0, 4326),
                ARRAY[
                ROW(1, '8BUI'::text, 0, 255),
                ROW(2, '8BUI'::text, 0, 255),
                ROW(3, '8BUI'::text, 0, 255)
                ]::addbandarg[]
        )
),
(2, ST_AddBand(
                ST_MakeEmptyRaster(240, 240, 82.4, 59.6, 6./3600., 6./3600., 0, 0, 4326),
                ARRAY[
                ROW(1, '8BUI'::text, 0, 255),
                ROW(2, '8BUI'::text, 0, 255),
                ROW(3, '8BUI'::text, 0, 255)
                ]::addbandarg[]
        )
) ;

-- Almost identical raster catalogue, the only difference being there is only one  band
DROP TABLE IF EXISTS dummy_rast_1band;
CREATE TABLE dummy_rast_1band(rid int, rast raster);
INSERT INTO dummy_rast_1band(rid, rast)
VALUES
(1, ST_AddBand(
                ST_MakeEmptyRaster(240, 240, 82,   59.6, 6./3600., 6./3600., 0, 0, 4326),
                '8BUI'::text, 0, NULL
        )
),
(2, ST_AddBand(
                ST_MakeEmptyRaster(240, 240, 82.4, 59.6, 6./3600., 6./3600., 0, 0, 4326),
                '8BUI'::text, 0, NULL
        )
) ;

-- Prepare polygon which touches one raster
DROP TABLE IF EXISTS dummy_poly;
CREATE TABLE dummy_poly(pid int, geom geometry);
INSERT INTO dummy_poly(pid, geom)
VALUES
(1, ST_GeomFromText( 'MULTIPOLYGON(((82.4180000000001 59.8655000000001,82.4 59.8655000000001,82.4 59.8745000000001,82.4180000000001 59.8745000000001,82.4180000000001 59.8655000000001)))', 4326));

-- Four tests with ST_Clip()
-- This works, with 1-band raster
SELECT ST_Summary(ST_Clip(r.rast, d.geom))
FROM dummy_rast_1band AS r INNER JOIN dummy_poly AS d
ON ST_Intersects(r.rast, d.geom);

-- This fails with
--   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
SELECT ST_Summary(ST_Clip(r.rast, d.geom))
FROM dummy_rast_3band AS r INNER JOIN dummy_poly AS d
ON ST_Intersects(r.rast, d.geom);

-- If ask for one band, that works
SELECT ST_Summary(ST_Clip(r.rast, 2, d.geom, true))
FROM dummy_rast_3band AS r INNER JOIN dummy_poly AS d
ON ST_Intersects(r.rast, d.geom);

-- If I union raster into bigger piece, it works
SELECT ST_Summary(ST_Clip(ST_Union(r.rast), d.geom))
FROM dummy_rast_3band AS r INNER JOIN dummy_poly AS d
ON ST_Intersects(r.rast, d.geom)
GROUP BY d.geom;

I think that if you would write a pgsql function that retrieves the data and processes it row by row you could catch the error and report on which row it occured. — Eelke, Feb 12 '17 at 06:43

score 1 · Answer 1 · answered Mar 07 '17 at 17:53

As I said in the EDIT above, I think the problem comes from unexpected behavior of ST_Clip when acting on particular condition: (1) raster has multiband and (2) polygon touches raster with no shared interior. My workaround it to operate each band separately. It appears to run slightly slower than the original, but it's not a catastrophe.

with pieces as (
select d.pid,
       st_clip(r.rast, 1, d.geom) as rast1
       st_clip(r.rast, 2, d.geom) as rast2
       st_clip(r.rast, 3, d.geom) as rast3
       from tbl_raster as r
       inner join tbl_poly as d
       on st_intersects(r.rast, d.geom)
)
insert into tbl_out(pid, val1, val2, val3)
select pid,
       (stats1).mean as v1,
       (stats2).mean as v2,
       (stats3).mean as v3
from (
     select p.pid,
            st_summarystatsagg(p.rast1, 1, true) as stats1,
            st_summarystatsagg(p.rast2, 1, true) as stats2,
            st_summarystatsagg(p.rast3, 1, true) as stats3
     from pieces as p
     group by p.pid
);

catch errors of postgis st_clip, and fix "Could not get band from working raster"

1 Answers1