Discrepancy between xarray-spatial hillshade and GDAL hillshade

[thuydotm]

I'm testing xarray-spatial against QGIS and seeing that when running hillshade on the same input data, the results by xarray-spatial and GDAL/QGIS are very different. Source code for QGIS hillshade can be found at: https://github.com/qgis/QGIS/blob/master/python/plugins/processing/algs/gdal/hillshade.py. This needs more research to see how the algorithm was implemented in QGIS to clearly understand the difference between the 2 libraries.

Input data:

array([[      nan,       nan,       nan,       nan,       nan,       nan],
       [704.237  , 242.24084, 429.3324 , 779.8816 , 193.29506, 984.6926 ],
       [226.56795, 815.7483 , 290.6041 ,  76.49687, 820.89716,  32.27882],
       [344.8238 , 256.34998, 806.8326 , 602.0442 , 721.1633 , 496.95636],
       [185.43515, 834.10425, 387.0871 , 716.0262 ,  49.61273, 752.95483],
       [302.4271 , 151.49211, 442.32797, 358.4702 , 659.8187 , 447.1241 ],
       [148.04834, 819.2133 , 468.97913, 977.11694, 597.69666, 999.14185],
       [268.1575 , 625.96466, 840.26483, 448.28333, 859.2699 , 528.04095]],
      dtype=float32)

xarray-spatial hillshade:

array([[       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan, 0.75030494, 0.06941041, 0.90643436, 0.15474272,            nan],
       [       nan, 0.80836594, 0.72366774, 0.14052185, 0.774778  ,            nan],
       [       nan, 0.93396175, 0.7071851 , 0.42872226, 0.9455124 ,            nan],
       [       nan, 0.85551083, 0.6819584 , 0.46013114, 0.23561102,            nan],
       [       nan, 0.41484872, 0.3213355 , 0.5821109 , 0.21879822,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan]],
 dtype=float32)

QGIS/GDAL hillshade

array([[  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [107.84468 ,  70.09885 ,   0.      ,  17.661407,   0.      ,          0.      ],
       [ 80.06987 ,  71.644684,   0.      ,   0.      ,   0.      ,          0.      ],
       [ 85.574615, 106.36669 ,  96.23605 ,  28.27108 ,  90.29079 ,         85.07072 ],
       [ 81.44522 ,  77.092354,   8.479876,   0.      ,   0.      ,          0.      ],
       [ 62.541145,   2.647696,   0.      ,   0.      ,   0.      ,          6.515689],
       [ 74.07955 ,  78.71434 ,   0.      ,  84.590744,  34.814816,         44.81609 ]],
 dtype=float32)

[brendancol]

@thuydotm same azimuth and angle correct? Can you run with GDAL from the CLI and then post the command?

[thuydotm]

@brendancol sure, here is the command:

gdaldem hillshade /private/var/folders/pf/9s4nnc0d0z7ghsq4f972gd9h0000gn/T/processing_AitUZB/665b1e0551a047d181859f56c34e8ea3/OUTPUT.tif /private/var/folders/pf/9s4nnc0d0z7ghsq4f972gd9h0000gn/T/processing_AitUZB/4009f64aa4fc4ed38fcb7f9590a0dea5/OUTPUT.tif -of GTiff -b 1 -z 1.0 -s 1.0 -az 225.0 -alt 25.0

[brendancol]

@thuydotm I think the will require looking at the GDAL slope and aspect functions as well to see how those non-trivial calcs are made:

• compare / contrast with GDAL slope implementation
• compare / contrast with GDAL aspect implementation

[remi-braun]

Here is some results I have, hope this helps :)

TL;DR:

• error comes from the gradient computation: the resolution is missing
• output is scaled differently
• xarray-spatial's numpy function can be simplified to avoid unnecessary trigonometric steps but seems coherent with GDAL.

---

Coherence between the core functions

Looking at (only) the numpy algorithm from xarray-spatial:

def _run_numpy(data, azimuth=225, angle_altitude=25):
    data = data.astype(np.float32)

    azimuth = 360.0 - azimuth
    x, y = np.gradient(data)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-x, y)
    azimuthrad = azimuth*np.pi/180.
    altituderad = angle_altitude*np.pi/180.
    shaded = np.sin(altituderad) * np.sin(slope) +  np.cos(altituderad) * np.cos(slope) *  np.cos((azimuthrad - np.pi/2.) - aspect)
    result = (shaded + 1) / 2
    result[(0, -1), :] = np.nan
    result[:, (0, -1)] = np.nan
    return result

Aspect

xarray-spatial has the same aspect function than GDAL's even if they are not written exactly the same way

Core function (shaded)

0. azimuth = 360.0 - azimuth is useless since azimuth is only used in a cosinus and cos(azimuth) = cos(2pi - azimuth)
1. If we ingest the slope and aspect in the shaded function:

shaded =  np.sin(alt_rad) * np.sin(np.pi/2. - np.arctan(np.sqrt(x2_y2))) + np.cos(alt_rad) * np.cos(np.pi/2. - np.arctan(np.sqrt(x2_y2))) * np.cos((az_rad - np.pi/2.) - aspect)

2. This can be simplified, knowing that:

• np.sin(np.pi/2. - np.arctan(np.sqrt(x2_y2))) = np.cos(np.arctan(np.sqrt(x2_y2)))
• np.cos(np.pi/2. - np.arctan(np.sqrt(x2_y2))) = np.sin(np.arctan(np.sqrt(x2_y2)))
• np.cos((az_rad - np.pi/2.) - aspect) = - np.sin((aspect - az)

into:

shaded =  np.sin(alt_rad) * np.cos(np.arctan(np.sqrt(x2_y2))) - np.cos(alt_rad) * np.sin(np.arctan(np.sqrt(x2_y2))) * np.sin(aspect - az_rad)

3. Moreover, since cos(atan(x)) = 1 / sqrt(1+x^2) and sin(atan(x)) = x / sqrt(1+x^2)

shaded =  np.sin(alt_rad) * 1 / sqrt(1+x2_y2) - np.cos(alt_rad) * sqrt(x2_y2) / sqrt(1+x2_y2) * np.sin(aspect - az_rad)

which can be factorised into:

shaded=(np.sin(alt_rad) - np.cos(alt_rad) * sqrt(x2_y2) * np.sin(aspect - az_rad)/sqrt(1+x2_y2)

4. Knowing that here aspect = atan2(-x, y) and in gdal aspect = atan2(x, y) which is the opposite, we can derive this fct.

shaded=(np.sin(alt_rad) - np.cos(alt_rad) * sqrt(x2_y2) * np.sin(aspect - az_rad) / sqrt(1 + x2_y2)

This is what we can see in GDAL 1.7.2, here

shaded=(np.sin(alt_rad) - np.cos(alt_rad) * np.sqrt(x2_y2) * np.sin(aspect - az_rad)) / np.sqrt(1 + x2_y2)

5. Newer versions of GDAL are going further into simplification, simplifying also the aspect sinus.
Here you have the whole simplification process https://github.com/OSGeo/gdal/blob/e4d1f7ff474ceeffd69b1ceef3fa635428788e8c/apps/gdaldem_lib.cpp#L770

sin(aspect - az) = sin(aspect)*cos(az) - cos(aspect)*sin(az))

and as sin(aspect)=sin(atan2(y,x)) = y / sqrt(xx_plus_yy)
and cos(aspect)=cos(atan2(y,x)) = x / sqrt(xx_plus_yy)

sin(aspect - az) = (y * cos(az) - x * sin(az)) / sqrt(xx_plus_yy)

---

Gradient issue

I replicated the GDAL function here:

    ds = rasterio.open("dem.tif")
    array = ds.read(masked=True)

    # Compute angles
    az_rad = azimuth * DEG_2_RAD
    alt_rad = (90 - zenith) * DEG_2_RAD

    # Compute slope and aspect
    dx, dy = np.gradient(np.where(array.mask, 0.0, array.data), *ds.res)
    x2_y2 = dx**2 + dy**2
    aspect = np.arctan2(dx, dy)

    # Compute hillshade (GDAL algo)
    hshade = (
        np.sin(alt_rad) + np.cos(alt_rad) * np.sqrt(x2_y2) * np.sin(aspect - az_rad)
    ) / np.sqrt(1 + x2_y2)
    hshade = np.where(hshade <= 0, 1.0, 254.0 * hshade + 1)

I can validate it by displaying in QGIS the input DEM with hillshade and mine (they are similar):

You can see the dataset's resolution inside the gradient function.

If I remove it and scale the two hillshades (mine and xarray-spatial's), I get a similar output:

Here you can find the used dataset:
hillshade.zip

[remi-braun]

Here is a proposition for the updated code, but I'm using rioxarray and it has a GDAL dependency to retrieve the resolution, so I won't make any PR because I don't know how to retrieve it otherwise.

from functools import partial
import rioxarray as rxr

xds = rxr.open_rasterio("dem.tif")
hillshade(xds, azimuth=315, zenith=45)

def hillshade(
    xds: xr.DataArray, azimuth: float, angle_altitude: float, **kwargs
) :
    def _run_numpy(data, azimuth, angle_altitude, res):
        # Compute angles
        az_rad = azimuth * DEG_2_RAD
        alt_rad = angle_altitude * DEG_2_RAD
    
        # Compute slope and aspect
        dx, dy = np.gradient(data, *res)
        x2_y2 = dx ** 2 + dy ** 2
        aspect = np.arctan2(dx, dy)
    
        # Compute hillshade (GDAL algo)
        hshade = (np.sin(alt_rad) + np.cos(alt_rad) * np.sqrt(x2_y2) * np.sin(aspect - az_rad)) / np.sqrt(1 + x2_y2)
        hshade = np.where(hshade <= 0, 1.0, 254.0 * hshade + 1)
    
        hshade[(0, -1), :] = np.nan
        hshade[:, (0, -1)] = np.nan
    
        return hshade
    
    _func = partial(_run_numpy, azimuth=azimuth, angle_altitude=angle_altitude, res=np.abs(xds.rio.resolution()))
    out = xds.data.map_overlap(
    out = xds.data.map_overlap(
        _func,
        depth=(1, 1),
        boundary=np.nan,
        meta=np.array(())
    )

    xds = xds.copy(data=out).rename(kwargs.get("name", "hillshade"))

    return xds