[raster] Instability with st_intersection/st_intersects

[kib]

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"

[kib]

Please also not that the ordering/sorting does not work

[dustymugs]

Try switching around the arguments in ST_Intersects() from raster, geometry to geometry, raster. It could be that there is a different answer resulting from the rasterization of the geometry. Have you tried PostGIS 2.1 from svn trunk?

Ordering/sorting not working? You'd need to elaborate on that…

[dustymugs]

Also check the output of the subqueries of bar and foo…

		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)

[pracine]

Could also be that nodata values are distributed differently in the two bands.

[kib]

Replying to pracine:

Could also be that nodata values are distributed differently in the two bands.

In this case NODATA is at the same positions in both grids - they are equal wrt. nodata, position and tile size.

[kib]

Replying to dustymugs:

Try switching around the arguments in ST_Intersects() from raster, geometry to geometry, raster. It could be that there is a different answer resulting from the rasterization of the geometry.

There is no difference when switching around geometry and raster

Have you tried PostGIS 2.1 from svn trunk?

No - I have only tried with 2.0.x

Ordering/sorting not working? You'd need to elaborate on that…

The SQL I gave had an order by lon,lat but looking at the output you see that the results is not in that ordering I can sort myself, so that is of minor import, but may be important in figuring out what goes on…

[kib]

Replying to dustymugs:

Also check the output of the subqueries of bar and foo…

Sorry but I do not read WKT, but if some here this it is relevant, say so and I will supply the data!

[dustymugs]

Replying to kib:

Replying to dustymugs:

Also check the output of the subqueries of bar and foo…

Sorry but I do not read WKT, but if some here this it is relevant, say so and I will supply the data!

Just see how many records are returned in the subquery.

[dustymugs]

Ordering/sorting not working? You'd need to elaborate on that…

The SQL I gave had an order by lon,lat but looking at the output you see that the results is not in that ordering I can sort myself, so that is of minor import, but may be important in figuring out what goes on…

Be more specific instead of just "order by lon, lat" since there are two of each. ORDER BY bar.lon, bar.lat.

[dustymugs]

Any way to provide a concise test case?

[kib]

Replying to dustymugs:

Replying to kib:

Replying to dustymugs:

Also check the output of the subqueries of bar and foo…

Sorry but I do not read WKT, but if some here this it is relevant, say so and I will supply the data!

Just see how many records are returned in the subquery.

OK in that case there are as many 15 and 12 - as in the with-queries.

[pracine]

This happens because ST_Intersection() return less polygons for one raster because it returns neigbour pixel sharing a same value as only one polygon. Use ST_PixelAsPolygons(ST_Clip(d52_2012112701_4.rast, pos.geom)) instead of ST_Intersection(d52_2012112701_4.rast, pos.geom)…

[kib]

Replying to dustymugs:

Ordering/sorting not working? You'd need to elaborate on that…

The SQL I gave had an order by lon,lat but looking at the output you see that the results is not in that ordering I can sort myself, so that is of minor import, but may be important in figuring out what goes on…

Be more specific instead of just "order by lon, lat" since there are two of each. ORDER BY bar.lon, bar.lat.

Makes no difference - just tried:-(

[dustymugs]

Be more specific instead of just "order by lon, lat" since there are two of each. ORDER BY bar.lon, bar.lat.

Makes no difference - just tried:-(

I'm betting this has something to do with floating point issues. Try something like

ORDER BY round(bar.lon::numeric(16,12), bar.lat::numeric(16,12))

[dustymugs]

Replying to pracine:

This happens because ST_Intersection() return less polygons for one raster because it returns neigbour pixel sharing a same value as only one polygon. Use ST_PixelAsPolygons(ST_Clip(d52_2012112701_4.rast, pos.geom)) instead of ST_Intersection(d52_2012112701_4.rast, pos.geom)…

Good call pracine. I forgot how ST_Intersection works. This would explain the behavior.

[kib]

Replying to pracine:

This happens because ST_Intersection() return less polygons for one raster because it returns neigbour pixel sharing a same value as only one polygon. Use ST_PixelAsPolygons(ST_Clip(d52_2012112701_4.rast, pos.geom)) instead of ST_Intersection(d52_2012112701_4.rast, pos.geom)…

Cool that works - and does sort fine too!

lon	lat	current_dir	current_speed
14.7918498789346	54.825	null	null
14.7918498789346	54.875	null	null
14.7918498789346	54.925	null	null
14.7918498789346	54.975	null	null
14.8751840193705	54.825	57.5229305513203	0.0245266498532146
14.8751840193705	54.875	31.0229305513203	0.0245266498532146
14.8751840193705	54.925	2.64793055132031	0.0167141498532146
14.8751840193705	54.975	330.39793055132	0.0108547748532146
14.9585181598063	54.825	52.8979305513203	0.0186672748532146
14.9585181598063	54.875	22.5229305513203	0.0225735248532146
14.9585181598063	54.925	351.39793055132	0.0206203998532146
14.9585181598063	54.975	346.27293055132	0.0128078998532146
(12 rows)

A little strange though with the null values:

select ST_X(g.geom) as lon, 
ST_Y(g.geom) as lat, 
ST_Value(d1142_2013031301_4.rast,1,g.geom) as "current-dir", 
ST_Value(d1142_2013031301_5.rast,1,g.geom) as "current-speed" 
from d1142_2013031301_4 
join d1142_2013031301_5 using (rid)  
CROSS JOIN (VALUES (ST_GeomFromText('POINT(14.79 54.825)',4326))) As g(geom) 
where ST_Intersects(d1142_2013031301_4.rast, g.geom) 
order by lon, lat;

  lon  │  lat   │   current-dir    │   current-speed    
───────┼────────┼──────────────────┼────────────────────
 14.79 │ 54.825 │ 79.7729305513203 │ 0.0284328998532146

[dustymugs]

Nothing to fix here.

[pracine]

A little strange though with the null values:

The NULL come out because I didn't set a band number in ST_PixelAsPolygons or ST_Clip. Just set 1 as band number in one or the other and NULL won't come out.

And next time write to the list before filling a ticket

[pracine]

Replying to pracine:

A little strange though with the null values:

The NULL come out because I didn't set a band number in ST_PixelAsPolygons or ST_Clip. Just set 1 as band number in one or the other and NULL won't come out.

Wrong. ST_PixelAsPolygon returns all pixel NULL or not. So just add a WHERE clause.

[kib]

Replying to pracine:

And next time write to the list before filling a ticket

Thanks for helping! Sorry about the noise.

[kib]

Replying to pracine:

Replying to pracine:

A little strange though with the null values:

The NULL come out because I didn't set a band number in ST_PixelAsPolygons or ST_Clip. Just set 1 as band number in one or the other and NULL won't come out.

Wrong. ST_PixelAsPolygon returns all pixel NULL or not. So just add a WHERE clause.

OK - I will just have to clip in SQL too :-/ (I have NODATA so nulls are a valid result)

where lon between 14.8 and 15 and lat between 54.8 and 55.0

Thanks for helping out!

[pracine]

Wrong. ST_PixelAsPolygon returns all pixel NULL or not. So just add a WHERE clause.

OK - I will just have to clip in SQL too :-/ (I have NODATA so nulls are a valid result)

where lon between 14.8 and 15 and lat between 54.8 and 55.0

I meant "WHERE NOT (d52_2012112701_4).val IS NULL".