es.hillshade()

[kstueb]

The function hillshade = es.hillshade(arr)...
...from earthpy.spatial yields a correct hillshade raster only if the arr represents a DTM with a spatial resolution of 1m.

This is easily fixed if the line x, y = np.gradient(arr) in the function hillshade() is replaced with x, y = np.gradient(arr, resolution), where resolution is the spatial resolution of the DTM passed to the function as additional parameter.

https://earthpy.readthedocs.io/en/latest/_modules/earthpy/spatial.html#hillshade

Oversteepened hillshade generated with es.hillshade() from a 25m resolution DTM:

Correct hillshade from the same DTM calculated with the proposed correction to es.hillshade():

Both hillshades plotted with ep.plot_bands()

[kstueb]

Proposed corrected function:

def hillshade(arr, azimuth=30, altitude=30, resolution=1):
    """Create hillshade from a numpy array containing elevation data.

    Parameters
    ----------
    arr : numpy array of shape (rows, columns)
        Numpy array with elevation values to be used to created hillshade.
    azimuth : float (default=30)
        The desired azimuth for the hillshade.
    altitude : float (default=30)
        The desired sun angle altitude for the hillshade.
    resolution : float (default=1)
        The spatial resolution of arr in the same units as the elevation values. 

    Returns
    -------
    numpy array
        A numpy array containing hillshade values.

    Example
    -------
    .. plot::

        >>> import matplotlib.pyplot as plt
        >>> import rasterio as rio
        >>> import earthpy.spatial as es
        >>> from earthpy.io import path_to_example
        >>> with rio.open(path_to_example('rmnp-dem.tif')) as src:
        ...     dem = src.read()
        >>> print(dem.shape)
        (1, 187, 152)
        >>> squeezed_dem = dem.squeeze() # remove first dimension
        >>> print(squeezed_dem.shape)
        (187, 152)
        >>> shade = es.hillshade(squeezed_dem)
        >>> plt.imshow(shade, cmap="Greys")
        <matplotlib.image.AxesImage object at 0x...>
    """
    try:
        # Possibly allowing different resolutions in x and y:
        # Note that resolution would have to be passed as a tuple (-y, x)! 
        #if isinstance(resolution, (list, tuple)):
        #    x, y = np.gradient(arr, *resolution) 
        #else:
        x, y = np.gradient(arr, resolution)
    except ValueError:
        raise ValueError("Input array should be two-dimensional")

    if azimuth <= 360.0:
        azimuth = 360.0 - azimuth
        azimuthrad = azimuth * np.pi / 180.0
    else:
        raise ValueError(
            "Azimuth value should be less than or equal to 360 degrees"
        )

    if altitude <= 90.0:
        altituderad = altitude * np.pi / 180.0
    else:
        raise ValueError(
            "Altitude value should be less than or equal to 90 degrees"
        )

    slope = np.pi / 2.0 - np.arctan(np.sqrt(x * x + y * y))
    aspect = np.arctan2(-x, y)

    shaded = np.sin(altituderad) * np.sin(slope) + np.cos(
        altituderad
    ) * np.cos(slope) * np.cos((azimuthrad - np.pi / 2.0) - aspect)

    return 255 * (shaded + 1) / 2