捕獲postgis st_clip的錯誤，並修復「無法從工作柵格獲得帶子」

Question

我擁有數百萬個多邊形，以及一個更大的柵格（平鋪）表。我想要st_clip（光柵，多邊形），然後在裁剪的光柵上應用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塊類型的東西，然後我要從那裏轉儲一些信息。或者讓它給一些空值或什麼，並使查詢完成，然後我看問題的輸出。

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

謝謝！

編輯：

繼@ Eelke的建議下，我用PL/pgSQL的環路，並確定問題的所在。

似乎如果多邊形接觸沒有共享內部點的柵格，ST_Clip會失敗，如果有多個頻段計算請求。如下面的測試代碼所示，如​​果我計算ST_與多邊形交互的柵格聯合，它似乎按預期工作（由我）。但結果是代價很高。我希望找到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 
);