6

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:

enter image description here

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.

  • just a thougtht; does it make a difference if you use a true CROSS JOIN (it shouldn´t, though, but the binding is strong in you, young CROSS JOIN....)? – ThingumaBob Apr 19 '18 at 11:23
  • @ThingamaBob. It doesn't, no. I would be quite concerned if it did, as I have always understood that the , is rewritten as CROSS JOIN by the query planner. – John Powell Apr 19 '18 at 12:37
  • yeah, I know...it only makes a difference for referencing tables in multiple joins by condition. why does ST_NearestValue return two rows? can you get a ST_Value if you use columnx, rowy instead? – ThingumaBob Apr 19 '18 at 12:49
  • @ThingumaBob. I have no idea why it returns two rows (this is the point of my question :-)), and more to the point, why one of them is -3.40282346638529e+38. It seems like an out by one error based on tiling. – John Powell Apr 19 '18 at 13:57
2

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.

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.