postgresql - 捕获 postgis st_clip 的错误，并修复 "Could not get band from working raster"

Question

我有数百万个多边形和一个更大的光栅表(平铺)。我想 st_clip(raster, polygon)，然后在裁剪的栅格上应用 st_summarystatsagg()。代码如下所示。

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

该代码适用于一小部分多边形(数千个)，但是当我超过数百万个时，它会返回错误

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

我想获得更多有关多边形导致此问题的诊断消息。有什么方法可以做 try-catch block 类型的事情然后我将从那里转储一些信息。或者让它给出一些空值或其他东西并使查询完成，然后我查看问题的输出。

postgresql 的版本是 9.6.1，postgis 是 2.3.1，gdal 是 2.1.2，在 Linux 上使用 psql -f 运行它。

谢谢!

编辑:

按照@Eelke 的建议，我在PL/pgSQL 中使用for 循环并确定了问题所在。

似乎如果多边形接触没有共享内部点的栅格，ST_Clip 会失败，如果请求了多个波段计算。正如下面的测试代码所示，如​​果我计算 ST_Intersects 与多边形的栅格并集，它似乎(按我的预期)工作。但事实证明这是代价高昂的。我希望找到 ST_Intersects 的替代方案，如果没有共享内部点，它不会拾取东西。但是我找不到执行此操作的功能，也无法弄清楚表达方式……有这样的表达方式吗？

下面是我想出的可重现的结果:

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

Answer 1

正如我在上面的编辑中所说，我认为问题来自 ST_Clip 在特定条件下的意外行为:(1) 栅格具有多波段和 (2) 多边形接触没有共享内部的栅格。我的解决方法是分别操作每个波段。它似乎比原来的运行速度稍慢，但这并不是灾难。

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