[postgis-tickets] [PostGIS] #3483: Problem with ST_Union
PostGIS
trac at osgeo.org
Thu Mar 17 01:22:10 PDT 2016
#3483: Problem with ST_Union
-----------------------+-----------------------------------------
Reporter: krzysiek | Owner: dustymugs
Type: defect | Status: new
Priority: medium | Milestone: PostGIS 2.1.9
Component: raster | Version: 2.1.x
Resolution: | Keywords: ST_Union, raster, alignment
-----------------------+-----------------------------------------
Comment (by krzysiek):
A little update on the issue.
I've been struggling trying to workaround the problem and it seems that
90% of cases I managed to resolve using ST_MapAlgebraExpr and
incrementally merging rasters in a loop. It is extremely slow but at least
I managed to process some of the data (not sure whether it workaround all
of the problems).
Another issue (or maybe a clue) I've found processing data obtained from a
rasterized vector layer using ST_AsRaster. For some of pieces of data
ST_Union works fine but in some cases fails, below you can find the piece
of code that is failing due to lack of alignment :
{{{
with polys as
(
select ST_Union(rast,'LAST') as rast from
(
select
ST_Resample(ST_Clip(ST_AsRaster(ST_Intersection(c2d.rect,p_aoi),p_ref_raster,ARRAY['8BUI',
'8BUI', '8BUI'],
case
when c2d.value=10 then ARRAY[1,255,1]
when c2d.value=20 then ARRAY[255,255,1]
when c2d.value=30 then ARRAY[255,127,1]
else ARRAY[255,255,255]
end
, ARRAY[0,0,0]),p_aoi),p_ref_raster) as rast
from (
select * from
classified2d c2d1
where c2d1.rect && p_aoi
and c2d1.value<>40
and c2d1.mean_depth>0) as c2d
where
--
}}}
I've copied problematic data into table rasters.r and tested it using
following code:
{{{
-- Wykonanie zapytania:
DO $$
declare
l_buffer raster;
i integer;
j integer;
BEGIN
FOR i IN 1..192 LOOP
select ST_Union(rast) into l_buffer from
(select filename,rast
from rasters.r
-- order by filename
) as a
where filename::integer<i;
RAISE NOTICE '% % % ',i,ST_Width(l_buffer),ST_Height(l_buffer);
end loop;
END $$;
}}}
The result is as follows:
{{{
UWAGA: 1 <NULL> <NULL>
UWAGA: 2 73 78
UWAGA: 3 230 241
UWAGA: 4 593 586
UWAGA: 5 1165 1153
UWAGA: 6 1534 1379
UWAGA: 7 1534 1379
UWAGA: 8 1769 2342
UWAGA: 9 3009 2342
UWAGA: 10 3009 2374
UWAGA: 11 4542 3589
UWAGA: 12 4600 3589
UWAGA: 13 4600 3589
UWAGA: 14 4600 3589
UWAGA: 15 4600 3589
UWAGA: 16 4600 3589
UWAGA: 17 4600 3589
UWAGA: 18 4600 3589
UWAGA: 19 4600 3589
BŁĄD: rt_raster_from_two_rasters: The two rasters provided do not have
the same alignment
}}}
If I enable sorting :
{{{
DO $$
declare
l_buffer raster;
i integer;
j integer;
BEGIN
FOR i IN 1..192 LOOP
select ST_Union(rast) into l_buffer from
(select filename,rast
from rasters.r
order by filename
) as a
where filename::integer<i;
RAISE NOTICE '% % % ',i,ST_Width(l_buffer),ST_Height(l_buffer);
end loop;
END $$;
}}}
Everything works fine (at least for this case):
{{{
UWAGA: 1 <NULL> <NULL>
UWAGA: 2 73 78
UWAGA: 3 230 241
UWAGA: 4 593 586
UWAGA: 5 1165 1153
UWAGA: 6 1534 1379
UWAGA: 7 1534 1379
UWAGA: 8 1769 2342
UWAGA: 9 3009 2342
UWAGA: 10 3009 2374
UWAGA: 11 4542 3589
UWAGA: 12 4600 3589
UWAGA: 13 4600 3589
UWAGA: 14 4600 3589
UWAGA: 15 4600 3589
UWAGA: 16 4600 3589
UWAGA: 17 4600 3589
UWAGA: 18 4600 3589
UWAGA: 19 4600 3589
UWAGA: 20 4678 3589
UWAGA: 21 4796 3589
.
.
.
UWAGA: 177 8337 7601
UWAGA: 178 8337 7601
UWAGA: 179 8337 7601
UWAGA: 180 8337 7601
UWAGA: 181 8337 7601
UWAGA: 182 8337 7601
UWAGA: 183 8337 7601
UWAGA: 184 8337 7601
UWAGA: 185 8337 7601
UWAGA: 186 8337 7601
UWAGA: 187 8337 7601
UWAGA: 188 8337 7601
UWAGA: 189 8337 7601
UWAGA: 190 8337 7601
UWAGA: 191 8337 7601
UWAGA: 192 8337 7601
Query returned successfully with no result in 05:57 minutes.
}}}
I've been testing it with different data but it didnt' help. From other
tests I've done (mostly changing the order of the rasters returned by the
internal query) the impression I have is that the result, for some reason,
depends on sequence in which the data arrive to ST_Union.
Are there any chances that somebody will take a look at the issue in near
future ? Maybe I am doing something particularly wrong in this code:
{{{
CREATE OR REPLACE FUNCTION get_rasterized_features1(
p_aoi geometry, -- area of interest
p_ref_raster raster) -- raster that will be background for the
data returned by this function
RETURNS raster AS
$BODY$
declare
--
l_tmp_classified raster;
l_tmp_patches raster;
l_tmp_features raster;
l_tmp_skeletons raster;
--
l_res raster;
--
begin
...
with classified as
(
select ST_Union(rast,'LAST') as rast from
--select rast as rast,classified2d_id from
(
select
ST_MakeEmptyRaster(ST_AddBand(ST_AddBand(p_ref_raster,p_ref_raster),p_ref_raster))
as rast, -1 as classified2d_id
union all
--select
ST_Resample(ST_Clip(ST_AsRaster(ST_Intersection(c2d.rect,p_aoi),p_ref_raster,ARRAY['8BUI',
'8BUI', '8BUI'],
select ST_AsRaster(c2d.rect,p_ref_raster,ARRAY['8BUI',
'8BUI', '8BUI'],
case
when c2d.value=10 then ARRAY[1,255,1]
when c2d.value=20 then ARRAY[255,255,1]
when c2d.value=30 then ARRAY[255,127,1]
else ARRAY[255,255,255]
end
, ARRAY[0,0,0]) as rast
,classified2d_id
from (
select * from
classified2d c2d1
where c2d1.rect && p_aoi
and c2d1.value<>40
and c2d1.mean_depth>0) as c2d
where
ST_Intersects(c2d.rect ,p_aoi)
and c2d.mean_length>5*c2d.mean_width/1000 and
ST_Area(rect)>0.00010681
order by 2 ) as p
where not ST_IsEmpty(rast)
)
--insert into rasters.r select classified2d_id,rast, 'Brak'::text
from classified;
select ST_Clip(rast,p_aoi) into l_tmp_classified from classified;
...
return l_tmp_classified
....
}}}
--
Ticket URL: <https://trac.osgeo.org/postgis/ticket/3483#comment:6>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.
More information about the postgis-tickets
mailing list