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