[raster] ST_Intersect returns false positives

[febus982]

Software and DB Structure: PostgreSQL 9.1.3 GEOS 3.3.2 GDAL 1.9.0

Table 1 (vector): Single column (geom) - About 100.000 points geometries exported with GRASS 6.4.2 (v.out.ogr command).

Table 2 (raster): Column 1 (rast): raster column Column 2 (filename): filename

Raster is geotiff file imported with: /usr/lib/postgresql/9.1/bin/raster2pgsql -s 4326 -a -Y -F -t 20x20

View 1:

SELECT v.geom, r.filename, st_value(r.rast, v.geom) AS value

FROM vector v, raster r

WHERE r.filename = 'file.tif'::text AND st_intersects(r.rast, 1, v.geom);

I got multiple points in this view values (some are null values) They are all aligned on same X and/or same Y.

If i materialize this view in a table using: CREATE TABLE table AS SELECT…… (same command) I got a bunch of notices in console: st_intersects NOTICE: Attempting to get pixel value with out of range raster coordinates:

View 2:

SELECT v.geom, r.filename, st_value(r.rast, v.geom) AS value

FROM geoserver_layers.z1 v, raster.raster r

WHERE r.filename = '2012042611_me.tif'::text AND st_intersects(r.rast, 1, v.geom) AND st_value(r.rast, v.geom) IS NOT NULL;

I Tried filtering out NULL values by adding "st_value(r.rast, v.geom) IS NOT NULL" but I got some points missing instead of multiple points;

This happens in all forms of ST_Intersects st_intersects(v.geom, r.rast, 1) st_intersects(r.rast, 1, v.geom)

I Tried also changing tile size in raster2pgsql

[Bborie Park]

Since you're using points, the point is being converted to a raster for ST_Intersects(raster, band, geometry). Hence the point becomes an area which may intersect with a raster that the point itself wouldn't.

As for ST_Intersects(geometry, raster, band), the raster is converted to polygons before running a geometry, geometry intersects test.

So by the looks of it, the problem is that in your query the ST_Intersects(raster, band, geometry) is operating by pixel (an "area") while ST_Value(raster, geometry) is essentially sticking a pin (your point) through the raster. The ST_Intersects returns true while ST_Value returns NULL.

You may want to use ST_Intersection(raster, band, geometry) instead of ST_Value.

[febus982]

I don't know the reason but i had to reboot and all seems working good now. Didn't have time to test with ST_Intersection.

Maybe some memory corruption :/