2

My query in (Dutch SRID, in meter):

 SELECT ST_NearestValue(ahn.rast, 1, ST_SetSRID(ST_Point(135742.6965165309,453864.2795337387), 28992)
   FROM ahn2.ahn
  WHERE ST_Intersects(ahn.rast, ST_SetSRID(ST_Point(135742.6965165309,453864.2795337387), 28992)::geometry, 1))

returns nothing in pgAdmin4:

enter image description here

I am sure the point is located in the middle of a raster. Slightly changing the coordinates results in:

  • 135742.6965165309, 453864.2795337387 (original)
  • 135740.6965165309, 453864.2795337387 (returns value)
  • 135744.6965165309, 453864.2795337387 (returns value)
  • 135742.6965165309, 453865.2795337387 (returns value)
  • 135742.6965165309, 453863.2795337387 (returns value)
  • 135741.6965165309, 453864.2795337387 (also nothing)

Why is my query returning nothing and how can it be fixed?

Background:

I'm working on a script that needs to calculate whether or not a certain locations has a shadow. To do that a line of sight (LOS) between the ground and the sun is calculated. This line is in fact a list of coordinates (x, y, and z). To determine whether or not there is shadow the height of the terrain needs to be determined for each x and y coordinate. If it is below the z (from the line) than the LOS is not blocked. The above query serves to extract the terrain height.

Version info

POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

  • what does just (WHERE ST_Intersects) return? – Moh May 7 '18 at 14:10
  • I don't know if it will make much difference, but you can certainly make you query easier to read by dropping some decimals and moving the search point into a CTE makes the spatial index behave better. WITH point(geom) AS (SELECT ST_SetSRID(ST_Point(135742.6965,453864.2795), 28992)) SELECT ST_NearestValue(rast, geom) FROM ahn2.ahn, point WHERE ST_Intersects(rast, geom); I have also noticed some odd behaviour with ST_NearestValue, but in my case it was to do with tiled rasters and the search point geometries being on tiled boundaries. – John Powell May 8 '18 at 10:28
  • @User2009, this returns a true/false for each row in the raster table (so for each tile, right?). As expected there is one 'true' and several thousand 'false'. Thanks. – LMB May 8 '18 at 18:07
  • @JohnPowellakaBarça, funny enough the proposed improvement of legibility also 'fixed' the problem for this particular set of coordinates. Unfortunately this query is run many times with many different coordinates. The problem still persist for other coordinates. – LMB May 8 '18 at 18:08
  • 1
    Yes, do you have an id field of some description, so that you can run distinct or max/group by on, as I suspect you will in fact have 100 rather than 106 distinct ids. I suspect it might have something to do with pathological behaviour on tile boundaries, though my understanding of C is not good enough to check this myself. See this possibly related question. – John Powell May 9 '18 at 21:10

Your Answer

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

Browse other questions tagged or ask your own question.