We would like to use the output of the udf for masking a cube, however, it can not filter a band from the udf ouput. We aim to download only the final ouput so this middle step (udf-hillshade mask) should be present as a variable in the entire flow.
This is the code with the udf:
start_date = '2021-05-01'
spatial_extent= {'west': -100.1202964782715, 'east': -99.9038314819336, 'south': 19.429686571587766, 'north': 19.494427786331926, 'crs': 4326}
## Get the Sentinel-2 data for a 3 month window.
start_date_dt_object = datetime.strptime(start_date, '%Y-%m-%d')
end_date = (start_date_dt_object + relativedelta(months = +1)).date() ## End date, 1 month later
start_date_exclusion = (start_date_dt_object + relativedelta(months = -1)).date() ## exclusion date, to give a 3 month window.
LOOKUPTABLE = {
"tropical": {
"S1": lambda vv: 1 / (1 + exp(- (-7.17 + (-0.48 * vv)))),
"S2": lambda ndvi, ndwi: 1 / (1 + exp(- (0.845 + (2.14 * ndvi) + (13.5 * ndwi)))),
"S1_S2": lambda vv, ndwi: 1 / (1 + exp(- (-2.64 + (-0.23 * vv) + (8.6 * ndwi)))),
},
"subtropical": {
"S1": lambda vv, vh: 1 / (1 + exp(- (-8.1 + (-0.13 * vv) + (-0.27 * vh)))),
"S2": lambda ndvi, ndwi: 1 / (1 + exp(- (0.845 + (2.14 * ndvi) + (13.5 * ndwi)))),
"S1_S2": lambda vv, ndwi: 1 / (1 + exp(- (-2.64 + (-0.23 * vv) + (8.6 * ndwi)))),
}
}
s2_cube = connection.load_collection(
'SENTINEL2_L2A_SENTINELHUB',
spatial_extent=spatial_extent,
temporal_extent=[start_date_exclusion, end_date],
bands=['B02', 'B03', 'B04', 'B08', 'sunAzimuthAngles', 'sunZenithAngles'])
s2_cube_masking = connection.load_collection(
'SENTINEL2_L2A_SENTINELHUB',
spatial_extent=spatial_extent,
temporal_extent=[start_date_exclusion, end_date],
bands=['CLP', 'SCL'])
scl = s2_cube_masking.band("SCL")
mask_scl = (scl == 3) | (scl == 8) | (scl == 9) | (scl == 10) | (scl == 11)
clp = s2_cube_masking.band("CLP")
mask_clp = mask_scl | (clp / 255) > 0.3
s2_cube = s2_cube.mask(mask_clp.resample_cube_spatial(s2_cube))
dem_cube = connection.load_collection(
"COPERNICUS_30",
spatial_extent = spatial_extent,
temporal_extent=["2010-01-01", "2030-12-31"],)
dem_cube = dem_cube.max_time()
dem_cube = dem_cube.resample_cube_spatial(s2_cube)
merged_cube = s2_cube.merge_cubes(dem_cube)
udf_code = """
from openeo.udf import XarrayDataCube
from openeo.udf.debug import inspect
import numpy as np
from hillshade.hillshade import hillshade
def rasterize(azimuth, resolution=None):
# Convert the azimuth into its components on the XY-plane. Depending on the value of the
# azimuth either the x or the y component of the resulting vector is scaled to 1, so that
# it can be used conveniently to walk a grid.
azimuth = np.deg2rad(azimuth)
xdir, ydir = np.sin(azimuth), np.cos(azimuth)
if resolution is not None:
xdir = xdir * resolution[0]
ydir = ydir * resolution[1]
slope = ydir / xdir
if slope < 1. and slope > -1.:
xdir = 1.
ydir = slope
else:
xdir = 1. / slope
ydir = 1.
return xdir, ydir
def _run_shader(sun_zenith, sun_azimuth, elevation_model, resolution_x, resolution_y):
azimuth = np.nanmean(sun_azimuth.astype(np.float32))
zenith = np.nanmean(sun_zenith.astype(np.float32))
if np.isnan(azimuth):
shadow[np.isnan(sun_azimuth)] = 255
else:
resolution = (float(resolution_x), float(resolution_y))
ray_xdir, ray_ydir = rasterize(azimuth, resolution)
# Assume chunking is already done by Dask
ystart = 0
yend = elevation_model.shape[0]
# Make sure inputs have the right data type
zenith = float(zenith)
ray = (float(ray_xdir), float(ray_ydir))
shadow = hillshade(elevation_model.astype(np.float32),
resolution,
zenith,
ray,
ystart,
yend)
shadow = shadow.reshape(elevation_model.shape)
shadow[np.isnan(sun_azimuth)] = 255
return shadow
def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
in_xarray = cube.get_array()
sun_zenith = in_xarray.sel({"bands": "sunZenithAngles"}).values.astype(np.float32)
sun_azimuth = in_xarray.sel({"bands": "sunAzimuthAngles"}).values.astype(np.float32)
elevation_model = in_xarray.sel({"bands": "DEM"}).values.astype(np.float32)
res_y = in_xarray.coords["y"][int(len(in_xarray.coords["y"])/2)+1] - in_xarray.coords["y"][int(len(in_xarray.coords["y"])/2)]
res_x = in_xarray.coords["x"][int(len(in_xarray.coords["x"])/2)+1] - in_xarray.coords["x"][int(len(in_xarray.coords["x"])/2)]
shadow = _run_shader(sun_zenith, sun_azimuth, elevation_model, res_x, res_y)
cube.get_array().values[0] = shadow
return cube
"""
process = openeo.UDF(code=udf_code, runtime="Python", data={"from_parameter": "x"})
# Shadow mask replaces the first band (B02) in the merged cube
hillshaded = merged_cube.apply(process=process)
hillshaded = hillshaded.filter_bands(bands=['B02'])
hillshaded_save = hillshaded_mask .save_result(format='netCDF')
my_job = hillshaded_save.send_job(title="hillshaded_mask ed")
results = my_job.start_and_wait().get_results()
results.download_files("hillshaded_mask ")
This is the problematic part:
hillshaded = merged_cube.apply(process=process)
hillshaded_mask = hillshaded.filter_bands(bands=['B02'])
Error:
OpenEoApiError: [500] Internal: Server error: OpenEoApiError('[500] Internal: Server error: ValueError("Invalid band name/index \'B02\'. Valid names: [\'DEM\']") (ref: r-23d11294a6574728a3fe951a014d4d88)') (ref: r-ab8cd419b1f24346b3d45ab3e7d9bd3b)
