array-to-array reducer support in `aggregate_spatial`

[soxofaan]

From https://discuss.eodc.eu/t/how-to-use-quantiles-on-a-2d-data-cube-x-y/703 I learned that the VITO backend supports array-to-array reducers in aggregate_spatial for example as documented at https://open-eo.github.io/openeo-python-client/basics.html#computing-multiple-statistics

Some issues

• This is not-standard usage outside of the openEO processes spec. So if we want to promote usage of this feature we should also take the initiative to push this in the aggregate_spatial spec
• when using array-to-array reducer with aggregate_spatial on a cube with multiple bands, the aggregations for all bands are flattened in a single band(?) dimension. See example below: 2 input bands with (mean, median, count) aggregation gives a band dimension with 6 "bands". It's undocumented for the user how to handle that. Also, aggregate_spatial has also a target_dimension parameter which could play a role here. This also important to take into account when we want to port aggregate_spatial to proper vector cubes

cube = con.load_collection(
    "SENTINEL2_L2A",
    temporal_extent=["2023-09-05", "2023-09-12"],
    bands=["B02", "B03"]
)
geometry = shapely.geometry.box(3, 51, 3.001, 51.001)

from openeo.processes import array_create
aggregations_multiple = cube.aggregate_spatial(
    geometries=geometry, 
    reducer=lambda data: array_create([data.mean(), data.median(), data.count()]),
)
aggregations_multiple.download("forum703-multiple.json")

results in something like

{
  "2023-09-05T00:00:00Z": [[280.33, 230, 84, 613.44, 565, 84]],
  "2023-09-07T00:00:00Z": [[410, 364, 84, 691.833, 634, 84]],
  "2023-09-10T00:00:00Z": [[621.2857, 584, 84, 863.1071, 833, 84]]
}