[postgis-tickets] [PostGIS] #2229: raster: Instability with st_intersection/st_intersects

[PostGIS]

#2229: raster: Instability with st_intersection/st_intersects
---------------------+------------------------------------------------------
 Reporter:  kib      |       Owner:  pramsey
     Type:  defect   |      Status:  new    
 Priority:  high     |   Milestone:         
Component:  postgis  |     Version:  2.0.x  
 Keywords:           |  
---------------------+------------------------------------------------------
 I have two rasters which are geographically equal - jost two bands from
 same grib-file.

 When I select an area from those the positions are different. I use this
 general SQL:
 {{{
 with bar as (select st_x(st_centroid((d52_2012112701_4).geom)) as lon,
        st_y(st_centroid((d52_2012112701_4).geom)) as lat,
        (d52_2012112701_4).val as "current_dir"
 from (select st_intersection(d52_2012112701_4.rast,1,pos.geom) as
 d52_2012112701_4
       from (select ST_GeomFromText('MULTIPOLYGON(((15.000 54.800,15.000
 55.000,14.800 55.000,14.800 54.800,15.000 54.800)))',4326) as geom) as pos
 join
            d52_2012112701_4 on
 st_intersects(d52_2012112701_4.rast,pos.geom) )bar),
 foo as (select st_x(st_centroid((d52_2012112701_5).geom)) as lon,
        st_y(st_centroid((d52_2012112701_5).geom)) as lat,
        (d52_2012112701_5).val as "current_speed"
 from (select st_intersection(d52_2012112701_5.rast,1,pos.geom) as
 d52_2012112701_5
       from (select ST_GeomFromText('MULTIPOLYGON(((15.000 54.800,15.000
 55.000,14.800 55.000,14.800 54.800,15.000 54.800)))',4326) as geom) as pos
 join
            d52_2012112701_5 on
 st_intersects(d52_2012112701_5.rast,pos.geom) ) foo)
 select bar.lon, bar.lat, bar."current_dir", foo."current_speed"
 from bar join foo using (lat,lon)
 order by lon,lat;
 }}}
 If I select from only one raster I get these points:
 || lon   ||  lat||current_dir          ||
 || 14.8167584745763 || 54.875 ||146.137054443359     ||
 || 14.8167584745763 || 54.975 || 120.262062072754    ||
 ||14.8167584745763 ||   54.8  ||145.637054443359         ||
 || 14.8167584745763 || 54.825 || 154.512054443359     ||
 || 14.8167584745763 || 54.925 ||145.262054443359    ||
 || 14.8751840193705 ||   54.8  ||157.012054443359    ||
 || 14.8751840193705 || 54.875  ||146.012054443359        ||
 || 14.8751840193705 || 54.825 || 153.262054443359        ||
 ||14.8751840193705 || 54.925  ||136.262054443359      ||
 || 14.8751840193705 || 54.975  ||112.887062072754         ||
 ||14.9584255447942  ||  54.8  ||167.012054443359        ||
 || 14.9584255447942 || 54.825  ||162.012054443359        ||
 || 14.9584255447942 || 54.925  ||143.887054443359      ||
 || 14.9584255447942 || 54.975  ||115.262062072754      ||
 || 14.9584255447942  ||54.875  ||154.637054443359       ||
 ||(15 rows)    ||

 While picking the other raster gives other points:
 ||       lon        ||  lat   ||   current_speed   ||
 || 14.8167584745763 || 54.875 || 0.119599260389805||
 || 14.8167584745763 || 54.975 || 0.115693010389805||
 || 14.8167584745763 || 54.825 || 0.154755517840385||
 || 14.8167584745763 || 54.925 || 0.113739885389805||
 || 14.8751840193705 || 54.825 || 0.117646135389805||
 || 14.8751840193705 || 54.925 || 0.107880510389805||
 || 14.8751840193705 || 54.975 || 0.154755517840385||
 ||             14.9 ||   54.8 || 0.131318017840385||
 || 14.9167584745763 || 54.875 || 0.105927385389805||
 || 14.9584255447942 || 54.825 || 0.115693010389805||
 || 14.9584255447942 || 54.925 || 0.113739885389805||
 || 14.9584255447942 || 54.975 || 0.172333642840385||
 ||(12 rows)||


 This instability does of cause gives problem when joining on positions
 where I would expect to have the exact same positions from both of these
 rasters :-( Here I only get ''' 10 rows''' where I would expect '''15'''


 Versions:
 "PostgreSQL 9.1.8 on x86_64-unknown-linux-gnu, compiled by gcc
 (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS="2.0.3 r11128"
 GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2,
 released 2012/10/08" LIBXML="2.8.0" LIBJSON="UNKNOWN" RASTER"

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/2229>
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.