3

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)

  • Sounds to me like the raster tiles are not actually in the right place... Do they appear at the correct co-ord location if you load them into a GIS? – Matt Dec 19 '17 at 11:10
1

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.

  • Thanks @cushen, this certainly makes for a neater query however unfortunately I am still seeing the same warnings. – Mark V Dec 20 '17 at 9:55
  • @MarkV too bad that didn't work; I've edited the answer, hopefully it should now work. I've tried out and example with points that don't intersect a raster here locally and all seemed ok. – Cushen Dec 20 '17 at 23:56
  • so sorry I completely missed your answer. It arrived as I stopped work for the Christmas holidays. I will be running the data load again this week so will try your suggested approach. Thanks so much for taking the time to help – Mark V Feb 6 '18 at 6:33
  • Thanks so much @cushen for the detailed answer and sorry it has taken me so long to get back to this. I was waiting on a new dataset which just arrived last week. The solution works perfectly thank you – Mark V Feb 12 '18 at 15:29
  • No worries @MarkV. Glad it worked out for you! – Cushen Feb 13 '18 at 6:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.