UDF output

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)

Hi Andrea,

could you put a ‘rename_labels’ between the apply and filter_bands?
Like:
hillshaded.rename_labels(“bands”, [“band1”, “band2”])

I believe that after the UDF, openEO has lost track of band names.

best regards,
Jeroen

Hi Jeroen,

unfortunately, It does not work.

hillshaded = merged_cube.apply(process=process)
hillshaded_band = hillshaded.rename_labels("bands", ["B02", "B03", "B04", "B08", "sunAzimuthAngles", "sunZenithAngles","DEM"])
hillshaded_mask = hillshaded_band.filter_bands(bands=["B02"])

hillshaded_save = hillshaded_mask.save_result(format='netCDF')
my_job  = hillshaded_save.send_job(title="hillshaded_mask")
results = my_job.start_and_wait().get_results()
results.download_files("hillshaded_mask")

Id: vito-j-6cbaf6fced7048dbbf096aac540932c5

Error message:
**OpenEoApiError** : [500] Internal: Server error: OpenEoApiError('[500] Internal: Server error: ValueError("Invalid band name/index \'B02\'. Valid names: [\'DEM\']") (ref: r-562c70f0ee3a47109bd49381954c3eda)') (ref: r-cb7ea7ad237d4f08b5a940c1e81c59c4)

FYI - I downloaded the ouput of udf, I can see that there are 7 bands included.
image

Hi Andrea,
I found the bug that causes this, and just committed a likely fix. It will take some time for automated testing before it becomes available on openeo-dev.

This is the issue:

One thing is that you will want to use apply_neighborhood as opposed to apply for the UDF, we discussed this via email with your colleague! (He has the working example now.)

best regards,
Jeroen

Thanks for looking into that. That was quick :slight_smile:

Do you think it is working now? @jeroen.dries