Calculate hillshade from DEM

hillshade.py

import cv2
import numpy as np


def get_hillshade(dem: np.ndarray, azimuth: float, zenith: float, resolution: int):
    """Calculates hillshade from DEM."""

    # calculate gradient using open-cv
    grad_x = cv2.Sobel(dem, cv2.CV_32F, 1, 0, ksize=3) / resolution / 8
    grad_y = cv2.Sobel(dem, cv2.CV_32F, 0, 1, ksize=3) / resolution / 8
    slope = np.arctan((grad_x**2 + grad_y**2) ** 0.5)  # slope in radians from 0 to pi

    # calculate aspect in radians from 0 to 2pi
    aspect = np.zeros_like(slope)
    aspect[grad_x != 0] = np.arctan2(grad_y[grad_x != 0], -grad_x[grad_x != 0])  # returns in range [-pi, pi]
    aspect[(grad_x != 0) & (aspect < 0)] += 2 * np.pi  # shift values for negative aspect
    aspect[(grad_x == 0) & (grad_y > 0)] = np.pi / 2  # set aspect for vertical gradient when grad_y > 0
    aspect[(grad_x == 0) & (grad_y < 0)] = 2 * np.pi - np.pi / 2  # set aspect for vertical gradient when grad_y < 0

    # calculate and return hillshade
    return (np.cos(zenith) * np.cos(slope)) + (np.sin(zenith) * np.sin(slope) * np.cos(azimuth - aspect))