Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 12 additions & 9 deletions openeo_driver/save_result.py
Original file line number Diff line number Diff line change
Expand Up @@ -314,8 +314,7 @@ def __init__(
"""
super().__init__(data=timeseries)
if not isinstance(regions, (GeometryCollection, DriverVectorCube)):
# TODO: raise exception instead of warning?
_log.warning(
raise ValueError(

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@soxofaan I did not find this warning in kibana. Safe to convert it to Exception?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes makes sense

f"AggregatePolygonResult: GeometryCollection or DriverVectorCube expected but got {type(regions)}"
)
self._regions = regions
Expand Down Expand Up @@ -445,17 +444,21 @@ def _create_point_timeseries_xarray(self, feature_ids, timestamps, lats, lons, a
return the_array

def to_netcdf(self, destination: Optional[str] = None) -> str:
def features_ids_from_index(geometries):
return ["feature_%d" % i for i in range(len(geometries.geoms))]

if isinstance(self._regions, GeometryCollection):
points = [r.representative_point() for r in self._regions.geoms]
feature_ids = features_ids_from_index(self._regions)
else:
feature_ids = ["feature_%d" % i for i in range(len(self._regions.geoms))]
elif isinstance(self._regions, DriverVectorCube):
points = [r.representative_point() for r in self._regions.get_geometries()]
feature_ids = self._regions.get_ids()
feature_ids = (list(feature_ids) if feature_ids is not None
else features_ids_from_index(self._regions.get_geometries()))
feature_ids = (
list(feature_ids)
if feature_ids is not None
else ["feature_%d" % i for i in range(self._regions.geometry_count())]
)
else:
raise ValueError(
f"AggregatePolygonResult: GeometryCollection or DriverVectorCube expected but got {type(self._regions)}"
)

lats = [p.y for p in points]
lons = [p.x for p in points]
Expand Down
26 changes: 26 additions & 0 deletions tests/test_save_result_netcdf.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,11 @@
from pathlib import Path

import geopandas as gpd
import numpy as np
import pytest
import xarray as xr
from numpy.testing import assert_array_equal, assert_allclose
from openeo_driver.datacube import DriverVectorCube
from openeo_driver.save_result import AggregatePolygonResult, AggregatePolygonResultCSV
from shapely.geometry import GeometryCollection, Polygon

Expand Down Expand Up @@ -150,3 +152,27 @@ def test_aggregate_polygon_result_CSV(tmp_path):
timeseries_ds.red.sel( t='2017-09-05')

assert_allclose( timeseries_ds.red.sel(feature=2).sel( t='2017-09-06').data,4645.719597)


def test_aggregate_polygon_result_vector_cube_no_ids_to_netcdf(tmp_path):
"""Test for https://github.com/Open-EO/openeo-python-driver/issues/519"""
timeseries = {
"2020-06-01T00:00:00Z": [[77.2], [77.15]],
"2020-06-02T00:00:00Z": [[77.21], [77.16]],
}
# GeoDataFrame without an 'id' column so that get_ids() returns None
gdf = gpd.GeoDataFrame(
geometry=[
Polygon([(77.11, 28.56), (77.29, 28.56), (77.29, 28.69), (77.11, 28.69)]),
Polygon([(77.11, 28.56), (77.20, 28.56), (77.20, 28.62), (77.11, 28.62)]),
]
)
regions = DriverVectorCube(gdf)

result = AggregatePolygonResult(timeseries, regions=regions)
filename = result.to_netcdf(tmp_path / "timeseries_vector_cube_no_ids.nc")

timeseries_ds = xr.open_dataset(filename)
# Feature IDs should be auto-generated as "feature_0", "feature_1", ...
assert list(timeseries_ds.coords["feature_names"].data) == ["feature_0", "feature_1"]
assert_array_equal(timeseries_ds.band_0.sel(feature=0).data, [77.2, 77.21])