Wednesday, 5 September 2018

PostGIS raster query "out of range raster coordinates"


I have a large raster (national scale) that I have loaded into PostGIS using raster2pgsql and a 100x100 tile size: raster2pgsql -I -C -s 27700 my.tif -t 100x100. I am now trying to query a large number of points to find the raster values at each location (values are held in a separate decode table keyed by the pixel value), my query is as follows:


INSERT INTO raster_values SELECT p.id, d.* FROM my_point_table p JOIN my_raster r ON ST_Intersects(p.the_geom, r.rast) JOIN my_raster_decode d ON (ST_Value(r.rast, p.the_geom) = d.value);


However when I run the query every so often I get a warning from PostGIS:



NOTICE: Attempting to get pixel value with out of range raster coordinates: (28, 100) NOTICE: Attempting to get pixel value with out of range raster coordinates: (64, 100) NOTICE: Attempting to get pixel value with out of range raster coordinates: (100, 1) NOTICE: Attempting to get pixel value with out of range raster coordinates: (100, 17)


The presence of a 100 in every error would indicate there is an issue with the bounds I am using but I am struggling to work out what I have done wrong. I have checked the SRID of both the point and raster table and confirmed that these are consistent (27700):


"enforce_srid_rast" CHECK (st_srid(rast) = 27700) the_geom | geometry(Geometry,27700)



Answer



Try the following:


INSERT INTO raster_values
SELECT pt.id, d.*
FROM
(
SELECT p.id id, ST_Intersection(rast,1,ST_AsRaster(p.the_geom, rast),1) as intersectx,

ST_Value(r.rast , p.the_geom) val
FROM my_raster r, my_point_table p
WHERE ST_Intersects(the_geom, rast)
) as pt,
my_raster_decode d
WHERE pt.intersectx IS NOT NULL AND pt.val = d.value;

ST_Intersects returns a true or false, it is not really joining, more filtering for where the geometries intersect.


EDIT:


I've edited the query adding a subquery to convert point to raster (for st_intersection) and also to calculate the ST_Value. Only Where this subquery is valid will the result be tested against the 'my_raster_decode' table.



If there is no intersection between the points and the raster then the query should return nothing, not just throw an error.


No comments:

Post a Comment