Viewshed computation of heightmap varies between DEM and GROUND heightmodes

[sa2329]

What is the bug?

When a viewshed heightmap is computed it should be different representations of the same viewshed, depending on whether DEM or GROUND heightmode is used. That is, the viewshed found with heightMode=DEM should equal heightMode=GROUND + elevation value from the DEM. This does not appear to be the case.

The attached figure shows several line-of-sight segments from a viewshed. The second example ("Observer on a local peak") shows inconsistency between GROUND and DEM - in both directions from the observer. The remainder show inconsistency in a single direction from the observer.

The GROUND values appear correct in all cases, while the DEM values appear to be below ground level at times.

Steps to reproduce the issue

Script to reproduce is included below, also attached as zip in case formatting is mangled.
script.zip

#!/usr/bin/env python
"""
Demonstrate viewshed calculations, focusing on edge cases
"""
from matplotlib import pyplot as plt
import numpy as np
from osgeo import gdal


gdal.UseExceptions()  #   Allow GDAL to throw exceptions

OBS_H = 0  # Observer height


# Test cases to examine
cases = [
    (2, [0, 0, 0, 1, 0, 0, 0, 0], "Simple calculation"),
    (3, [1, 1, 0, 1, 0, 1, 2, 2], "Observer on local peak"),
    (0, [0, 0, 0, 1, 1, 0, 0, 0], "Hill in line-of-sight"),
    (0, [0, 0, 1, 2, 3, 4, 5, 6], "Steady slope"), 
    (0, [0, 0, 1, 1.1, 3, 4, 5, 4], "Staggered slope"),
]

def calc_viewshed(vals, startx, heightMode):
    """
    vals : list of elevation values
    heightMode : gdal terrain mode to use in viewshed computation
    """
    # Create gdal dataset with elevation data
    arr = np.stack([vals, vals])  # Viewshed requires at least 2 rows, so duplicate data.  We'll only use the first row
    
    driver = gdal.GetDriverByName("MEM")
    src = driver.Create(
        "",  # empty filename, to create in memory
        xsize=arr.shape[1],
        ysize=arr.shape[0],
        bands=1,
        eType=gdal.GDT_Float32,
    )
    src.SetGeoTransform((-startx, 1, 0, 0.5, 0, -1))
    src.GetRasterBand(1).WriteArray(arr)

    ds = gdal.ViewshedGenerate(
        srcBand=src.GetRasterBand(1),
        driverName="MEM",
        targetRasterName="",
        creationOptions=[],
        observerX=0,
        observerY=0,
        observerHeight=OBS_H,
        targetHeight=-1,  # argument ignored in heightmap mode, but need to provide a value
        visibleVal=255,
        invisibleVal=0,
        outOfRangeVal=-1,
        noDataVal=-1,
        dfCurvCoeff=0,   # flat earth
        mode=gdal.GVM_Edge,
        maxDistance=1000,
        heightMode=heightMode,
    )
    return ds.ReadAsArray()[0]


fig, axes = plt.subplots(len(cases))
for [startx, vals, desc], ax in zip(cases, axes):
    # Create viewsheds
    vs_height_dem_crs = calc_viewshed(vals, startx, gdal.GVOT_MIN_TARGET_HEIGHT_FROM_DEM)
    vs_height_above_ground = calc_viewshed(vals, startx, gdal.GVOT_MIN_TARGET_HEIGHT_FROM_GROUND)
    vs_height_ground_plus = vals + vs_height_above_ground
     
    # Validate behavior
    if np.allclose(vs_height_dem_crs, vs_height_ground_plus):
        desc += ":\nCalculations using height from DEM and GROUND match"
    else:
        desc += ":\nInconsistency between calculations using height from DEM and GROUND"
        
    # Plot terrain
    ax.set_title(desc)
    x = np.arange(0, len(vals))
    ax.plot(x, vals, 'g-x', linewidth=1, label='Ground level ')

    # Target height visibility
    ax.plot(x, vs_height_dem_crs, ':mo', label='Heightmap: DEM CRS')
    ax.plot(x, vs_height_ground_plus, '--y.', label='Heightmap: GROUND')
    ax.plot(x, vs_height_above_ground, ':b', label='Heightmap:\nGROUND + DEM height')

    # Observer point
    h = OBS_H + vals[int(startx)]
    ax.plot([startx]*2, [vals[int(startx)], h], 'r:')
    ax.plot([startx], [h], 'rd', label="Observer location")

    ax.legend(title='Key', loc='upper left')
    
plt.show()

Versions and provenance

GDAL 3.8.3, released 2024/01/04
Python 3.9.18

Additional context

No response

[rouault]

CC @szekerest