Hi all,
I am working on a UDF that reproduces the gdal function to create the hillshade.
First, I get the DEM
eoconn = openeo.connect('https://openeo.dataspace.copernicus.eu/', auto_validate=False)
eoconn.authenticate_oidc()
dem_datacube = eoconn.load_collection(
"COPERNICUS_30",
spatial_extent={'west':11.241593,
'east':11.444267,
'south':46.914255,
'north':46.992429,
'crs':4326},
bands=["DEM"]
)
Then, I apply the udf
udf_path="./udf-hillshade.py"
hillshade_func = openeo.UDF.from_file(udf_path, context={"from_parameter": "context"})
hillshade = dem_datacube.apply(process=hillshade_func)
where the udf is defined in this way
import openeo
from openeo.udf.debug import inspect
from openeo.udf import XarrayDataCube
import math
import numpy as np
import xarray as xr
from skimage.restoration import inpaint
def apply_datacube(cube: XarrayDataCube,
context: dict) -> XarrayDataCube:
"""
UDF to get the hillshade, i.e., a shaded relief map.
It takes the illumination angles as input, i.e.,
altitude of the sun (default 315 degrees: the terrain is lit as if the
light source is coming from the upper left)
azimuth of the sun (default 45 degrees: the light is pitched 45 degrees
between the horizon and directly overhead)
First, the slope and aspect are calculated by following Horn (1981).
Input and output datacube have the same dimensions.
"""
altitude = 45
azimuth = 315
nodataval = 255
# conversion from degrees to radians
deg2rad = 180/np.pi
# zenith angle of the sun calculated from the altitude
zenith_rad = (90 - altitude)/deg2rad
# zenith angle of the sun calculated from the altitude
azimuth_rad = azimuth/deg2rad
# load the datacube
array = cube.get_array()
data_array = array.values[0]
# get the coodinates of the datacube
xmin = array.coords['x'].min()
xmax = array.coords['x'].max()
ymin = array.coords['y'].min()
ymax = array.coords['y'].max()
coord_x = np.linspace(start = xmin,
stop = xmax,
num = array.shape[-2])
coord_y = np.linspace(start = ymin,
stop = ymax,
num = array.shape[-1])
# get the datacube spatial resolution
resolution = array.coords['x'][-1] - array.coords['x'][-2]
# Calculate gradients
dz_dx = ((data_array[:-2, 2:] + 2 * data_array[1:-1, 2:] + data_array[2:, 2:]) -
(data_array[:-2, :-2] + 2 * data_array[1:-1, :-2] + data_array[2:, :-2]))/(8 * resolution.values)
dz_dy = ((data_array[2:, :-2] + 2 * data_array[2:, 1:-1] + data_array[2:, 2:]) -
(data_array[:-2, :-2] + 2 * data_array[:-2, 1:-1] + data_array[:-2, 2:])) / (8 * resolution.values)
# Compute slope and aspect
slope = np.arctan(np.sqrt(dz_dx**2 + dz_dy**2))
aspect = np.arctan2(dz_dy, -dz_dx)
aspect = np.where(aspect < np.pi/2, np.pi/2 - aspect, 5*np.pi/2 - aspect)
# Compute hillshade
hillshade = (np.cos(zenith_rad) * np.cos(slope) +
np.sin(zenith_rad) * np.sin(slope) * np.cos(azimuth_rad - aspect))
hillshade = np.clip(hillshade, 0, 1)
# Scale to 0-255
hillshade = nodataval * hillshade
hillshade[hillshade<0] = 0
hillshade[slope==0] = 0
# round to integers
hillshade = np.round(hillshade)
# Create output array and insert hillshade values
output = np.zeros_like(data_array) + nodataval
output[1:-1, 1:-1] = hillshade
# Create a mask for the 255 values
mask = (output == 255).astype(np.uint8)
# Replace the 255 values using biharmonic inpainting
interpolated_array = inpaint.inpaint_biharmonic(output, mask)
interpolated_array = interpolated_array.astype(np.uint8)
# predicted datacube: same dimensions as the input datacube
predicted_cube = xr.DataArray(interpolated_array,
dims=['x', 'y'],
coords=dict(x=coord_x, y=coord_y))
return XarrayDataCube(predicted_cube)
I have tested the UDF and it is working. However, I needed to add something to interpolate missing values on the edges. When I add this last step about the inpainting - last lines of the code-, it returns a partial output
You can see that sume chunks seem not to be processed. Can you help and tell me why?
Thanks,
Valentina
