Error with ST_Value, ST_NearestValue: out of range raster coordinates

Question

I have a UK 5m DTM coverage loaded into Postgis with:

raster2pgsql -s 27700 -I -C *.tif -t 100x100 ....

which looks more or less like the below for a full 1000 x 1000 pixel tile, before tiling:

This works as expected, however when I run a query for ST_Value, such as:

SELECT ST_Value(rast, geom, true) 
  FROM terrain_5, addresses 
 WHERE ST_Intersects(rast, geom) 
   AND address_id = xxxx;

I get the error:

NOTICE: Attempting to get pixel value with out of range raster coordinates: (2, 100)

If I run the same query, but using ST_NearestValue instead I get two rows back, eg,

-3.40282346638529e+38
 70.7099990844727

I have checked for the points where this happens and they are always round numbers divisible by 500.

If instead I request the metadata, ie, the same query but using (ST_MetaData(rast)).* I get two rows,

upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands

166500 | 27500 | 100 | 100 | 5 | -5 | 0 | 0 | 27700 | 1

166500 | 27000 | 100 | 100 | 5 | -5 | 0 | 0 | 27700 | 1

It is pretty obvious that the problem is due to 100x100 tiles of pixel size 5m, resulting in edge conditions every time a coordinate being tested for intersection is divisible by 500.

I have solved the problem by rewriting the above query using max(ST_NearestValue(rast, geom)), ie,

SELECT max(ST_NearestValue(rast, geom, true) 
  FROM terrain_5, addresses 
 WHERE ST_Intersects(rast, geom) 
   AND address_id = xxxx;

the final parameter makes no difference, as there are no no data values in the coverage, as can be seen from the input raster.

If I run the histogram, excluding, no data, I get and add up the values:

SELECT sum(count) 
  FROM 
      (SELECT 
            (ST_Histogram(rast, 1, true, 10)).count 
         FROM terrain_5, addresses
        WHERE ST_Intersects(rast, geom) 
         AND address_id = xxxx
      ) foo;

I get 20000, ie, the sum of the pixels in both the intersected tiles that intersect the point divisible by 500.

I have seen this question, but I don't think it really answers my case.

I have a solution, using max(ST_NearestValue(rast, geom)), but this does strike me as unexpected behaviour, or I am missing something really obvious. Sadly, when I ran this on the national dataset, the value that was chosen from the two results was the low one, so I have a lot of houses whose height above sea level is -3.40282346638529e+38, when even by British standards of building in flood plans, strikes me as excesssive. I only noticed when looking at summary statistics after the fact, as my original update query swallowed the errors.

Answer 1

This isn't really an answer as such but I think I have experienced this bug too and have reported it here.

It could also be something in one of the underlying libraries.

I will see what else I can do to try and debug.