GDALIsLineOfSightVisible(): points exactly on the DEM surface are never visible

[ctoney]

What is the bug?

A point exactly on the surface of a DEM is never visible with GDALIsLineOfSightVisible(), even if viewed from a point directly above. Is this the intended behavior?

Small code example is below. I've also done checks on a real DEM using Z value equal to the pixel value at the point XY. A point exactly on the surface is never visible.

This is not intuitive and not well documented. The behavior is also different from that of gdal_viewshed which treats the surface itself as visible, and uses 0 as the default height above the DEM surface for the target (-tz).

The documentation for GDALIsLineOfSightVisible() references a paper for the algorithm: https://www.researchgate.net/publication/2411280_Efficient_Line-of-Sight_Algorithms_for_Real_Terrain_Data.

The paper itself is contradictory on whether a point exactly on the surface intercepts the line-of-sight. In section 4. Height Map, it first says:

Thus, given the coordinates (x,y) of two combat elements, the z coordinate is obtained for the visibility calculation from the height map. After that, a 3D Bresenham algorithm is used to “traverse” the straight-line segment joining the two coordinates. At each step, the height of the current coordinate of the line segment is verified to check if it is higher (larger) than the corresponding height in the map.

And just below that:

However, in order to assert that the combat elements are not visible, it is enough that one of those heights is the same as or inferior to the height of the terrain.

This seems to be how the algorithm is implemented in GDALIsLineOfSightVisible().

But further down in 6. Conclusions, rather than "same as or inferior to the height of the terrain" it refers instead only to "a height smaller than the height of the image (height map)":

The visibility is calculated through an algorithm of discrete traversal along the straight-line segment between the two combat elements, such as, for example, Bresenham algorithm. In case it meets along the traversal a height smaller than the height of the image (height map), the combat elements will not be visible (see Figure 5); otherwise, they are visible.

That is, the first interpretation gives z > terrainHeight while the second would imply z >= terrainHeight.

The current implementation appears to require all line segments to be above the terrain, including the start and end points. So a small fudge value must always be added to the DEM elevation for Z in order to check the visibility of a point on the surface using GDALIsLineOfSightVisible(). This will be true even if the other point is directly overhead, or even above any point on the DEM surface (e.g., the position of an aerial vehicle).

Maybe this is due to the use case of aerial drone path planning where the surface itself would not be used directly and there would always be some height above the DEM surface that must be avoided? Even so, using z >= terrainHeight would not impact that use case if there is always some addition to the DEM values anyway. Otherwise, it seems awkward to use. Is this the expected/intended behavior?

Steps to reproduce the issue

>>> from osgeo import gdal

>>> mem_ds = gdal.GetDriverByName("MEM").Create("", 2, 1)

>>> res = gdal.IsLineOfSightVisible(mem_ds.GetRasterBand(1), 0, 0, 1, 0, 0, 0)
>>> res.is_visible
False

>>> res = gdal.IsLineOfSightVisible(mem_ds.GetRasterBand(1), 0, 0, 1, 0, 0, 0.01)
>>> res.is_visible
True

Versions and provenance

Linux Ubuntu 24.04 (ubuntugis-unstable package)

$ gdalinfo --version
GDAL 3.10.3, released 2025/04/01

Additional context

No response

[rouault]

CC @Ryanf55

[Ryanf55]

Hello thank you for the detailed report. This seems like a correct change. I'll fix it, add docs, and a test.

[ctoney]

Maybe start/end points that fall exactly on the surface should be visible only when the traversed points in between are all above terrain? (i.e., potentially still use z > terrainHeight for all checks in between, and z >= terrainHeight only for start/end points.)

I'm not sure.

Tradeoffs seem to be that it would make the the code slightly more complicated, but which method makes the most sense logically (should a tangential point in between intercept visibility?), and is floating point comparison a concern in this case with z >= terrainHeight?