Opened 16 months ago

Last modified 16 months ago

#3483 new defect

Problem with ST_Union

Reported by: krzysiek Owned by: dustymugs
Priority: blocker Milestone: PostGIS 2.1.9
Component: raster Version: 2.1.x
Keywords: ST_Union, raster, alignment Cc: kmatuk@…

Description

I have loaded some rasters into Postgis. I did a lot of ST_XXXX (scale, snap to grid,resample) magic to make them aligned. After executing:

select r1.rid,r2.rid ,ST_NotSameAlignmentReason(r1.rast,r2.rast) ,ST_SameAlignment(r1.rast,r2.rast) from

rasters.r r1, rasters.r r2

where not ST_SameAlignment(r1.rast,r2.rast)

I got nothing, if I remove:

where not ST_SameAlignment(r1.rast,r2.rast)

all I got is

"The rasters are aligned"

Which seems ok, but when I execute:

SELECT ST_Union(rast,'FIRST') from rasters.r

I got:

rt_raster_from_two_rasters: The two rasters provided do not have the same alignment

Why ?

Postgres 9.4

"POSTGIS="2.1.8 r13780" GEOS="3.5.0-CAPI-1.9.0 r4090" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER"

Change History (7)

comment:1 Changed 16 months ago by robe

Component: postgisraster
Owner: changed from pramsey to dustymugs

comment:2 Changed 16 months ago by robe

Do you have sample data you can share we can test with?

comment:3 Changed 16 months ago by dustymugs

Can you provide some additional information about the rasters?

How many rasters?

Can you provide the metadata of a tile and possibly of all the tiles (assuming there is a reasonable number)?

comment:4 Changed 16 months ago by krzysiek

Is dump of the table ok ?

Schema : rasters Table r filesize: 42MB:

www.imago3d.pl/rasters_r.zip

comment:5 Changed 16 months ago by krzysiek

Did you managed to get the data ?

comment:6 Changed 16 months ago 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 workarounds 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

 ....

Last edited 16 months ago by krzysiek (previous) (diff)

comment:7 Changed 16 months ago by krzysiek

Priority: mediumblocker
Note: See TracTickets for help on using tickets.