Problem with ST_Union

[krzysiek]

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"

[robe]

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

[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)?

[krzysiek]

Is dump of the table ok ?

Schema : rasters Table r filesize: 42MB:

www.imago3d.pl/rasters_r.zip

[krzysiek]

Did you managed to get the data ?

[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

 ....

[robe]

FWIW the test file doesn't seem to be present anymore so can't debug this. I'm tempted to mark this as a won't fix.

[krzysiek]

The problem was the accuracy of the georeference of the images (as far as I remember). The imaging system was used to visualise a surface of the road. It was build of two cameras. During the callibration of the system we got numbers with 14 significand digits. After fifth or sixth digit there were zeros only until 12 or latter. Since we were following the road the sysyem tend to cumulate error and it got bigger than one pixel after a while and this was detected as a problem described in this ticket. The occurence of the probem seemd to be random and I've noticed the source after a few sleepless nights. Looks like you can close the ticket.

Thank you and regards Krzysiek

[robe]

Closing because can't replicate