Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: support bbox in in geometry #138

Merged
merged 17 commits into from
Mar 21, 2024
Merged
Show file tree
Hide file tree
Changes from 16 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
14 changes: 14 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,20 @@ i = ee.ImageCollection(ee.Image("LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318"))
ds = xarray.open_dataset(i, engine='ee')
```

Open any Earth Engine ImageCollection to match an existing transform:

```python
raster = rioxarray.open_rasterio(...) # assume crs + transform is set
ds = xr.open_dataset(
'ee://ECMWF/ERA5_LAND/HOURLY',
engine='ee',
geometry=tuple(raster.rio.bounds()), # must be in EPSG:4326
projection=ee.Projection(
crs=str(raster.rio.crs), transform=raster.rio.transform()[:6]
),
)
```

See [examples](https://github.com/google/Xee/tree/main/examples) or [docs](https://github.com/google/Xee/tree/main/docs) for more uses and integrations.

## How to run integration tests
Expand Down
47 changes: 33 additions & 14 deletions xee/ext.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
import math
import os
import sys
from typing import Any, Dict, Iterable, List, Literal, Optional, Tuple, Union
from typing import Any, Dict, Iterable, List, Literal, Optional, Sequence, Tuple, Union
from urllib import parse
import warnings

Expand Down Expand Up @@ -139,7 +139,9 @@ def open(
crs: Optional[str] = None,
scale: Optional[float] = None,
projection: Optional[ee.Projection] = None,
geometry: Optional[ee.Geometry] = None,
geometry: Union[
ee.Geometry, Tuple[float, float, float, float], None
] = None,
primary_dim_name: Optional[str] = None,
primary_dim_property: Optional[str] = None,
mask_value: Optional[float] = None,
Expand Down Expand Up @@ -176,7 +178,9 @@ def __init__(
crs: Optional[str] = None,
scale: Union[float, int, None] = None,
projection: Optional[ee.Projection] = None,
geometry: Optional[ee.Geometry] = None,
geometry: Optional[
Union[ee.Geometry, Tuple[float, float, float, float]]
] = None,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Could you update the type annotation to the following. Geometry | Tuple | None is preferred over Optional[Geometry | Tuple].

Union[ee.Geometry, Tuple[float, float, float, float], None]

Here and other places in the PR.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

implemented here b640d95

primary_dim_name: Optional[str] = None,
primary_dim_property: Optional[str] = None,
mask_value: Optional[float] = None,
Expand Down Expand Up @@ -231,19 +235,31 @@ def __init__(
# Parse the dataset bounds from the native projection (either from the CRS
Copy link
Collaborator

Choose a reason for hiding this comment

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

Consider breaking this out into a separate function. This logic is much more complex now.

Looks like the only input is geometry. Return tuple should be a 4-length tuple.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

added in b640d95

# or the image geometry) and translate it to the representation that will be
# used for all internal `computePixels()` calls.
try:
if isinstance(geometry, ee.Geometry):
if geometry is None:
try:
x_min_0, y_min_0, x_max_0, y_max_0 = self.crs.area_of_use.bounds
except AttributeError:
# `area_of_use` is probable `None`. Parse the geometry from the first
Copy link
Collaborator

Choose a reason for hiding this comment

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

Nit: probably

Copy link
Contributor Author

Choose a reason for hiding this comment

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

implemented here b640d95

# image instead (calculated in self.get_info())
x_min_0, y_min_0, x_max_0, y_max_0 = _ee_bounds_to_bounds(
self.get_info['bounds']
)
else:
x_min_0, y_min_0, x_max_0, y_max_0 = self.crs.area_of_use.bounds
except AttributeError:
# `area_of_use` is probable `None`. Parse the geometry from the first
# image instead (calculated in self.get_info())
elif isinstance(geometry, ee.Geometry):
x_min_0, y_min_0, x_max_0, y_max_0 = _ee_bounds_to_bounds(
self.get_info['bounds']
)
elif isinstance(geometry, Union[List, Tuple]):
Copy link
Collaborator

Choose a reason for hiding this comment

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

You should be able to check against Sequence here instead.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is fine, but a string is also sequence and this would catch and raise in that case. Updated anyways.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

implemented here b640d95

if len(geometry) != 4:
raise ValueError(
'geometry must be a tuple or list of length 4, or a ee.Geometry, '
f'but got {geometry!r}'
)
x_min_0, y_min_0, x_max_0, y_max_0 = geometry
else:
raise ValueError(
'geometry must be a tuple or list of length 4, a ee.Geometry, or'
f' None but got {type(geometry)}'
)

x_min, y_min = self.transform(x_min_0, y_min_0)
x_max, y_max = self.transform(x_max_0, y_max_0)
Expand Down Expand Up @@ -586,7 +602,7 @@ def _get_tile_from_ee(
)
target_image = ee.Image.pixelCoordinates(ee.Projection(self.crs_arg))
return tile_index, self.image_to_array(
target_image, grid=bbox, dtype=np.float32, bandIds=[band_id]
target_image, grid=bbox, dtype=np.float64, bandIds=[band_id]
)

def _process_coordinate_data(
Expand Down Expand Up @@ -693,7 +709,7 @@ def _parse_dtype(data_type: types.DataType):


def _ee_bounds_to_bounds(bounds: ee.Bounds) -> types.Bounds:
coords = np.array(bounds['coordinates'], dtype=np.float32)[0]
coords = np.array(bounds['coordinates'], dtype=np.float64)[0]
x_min, y_min, x_max, y_max = (
min(coords[:, 0]),
min(coords[:, 1]),
Expand Down Expand Up @@ -948,7 +964,9 @@ def open_dataset(
crs: Optional[str] = None,
scale: Union[float, int, None] = None,
projection: Optional[ee.Projection] = None,
geometry: Optional[ee.Geometry] = None,
geometry: Optional[
Union[ee.Geometry, Tuple[float, float, float, float]]
] = None,
primary_dim_name: Optional[str] = None,
primary_dim_property: Optional[str] = None,
ee_mask_value: Optional[float] = None,
Expand Down Expand Up @@ -1007,7 +1025,8 @@ def open_dataset(
coalesce all variables upon opening. By default, the scale and reference
system is set by the the `crs` and `scale` arguments.
geometry (optional): Specify an `ee.Geometry` to define the regional
bounds when opening the data. When not set, the bounds are defined by
bounds when opening the data or a bbox specifying [x_min, y_min, x_max,
y_max] in EPSG:4326. When not set, the bounds are defined by
the CRS's 'area_of_use` boundaries. If those aren't present, the bounds
are derived from the geometry of the first image of the collection.
primary_dim_name (optional): Override the name of the primary dimension of
Expand Down
39 changes: 39 additions & 0 deletions xee/ext_integration_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -357,6 +357,45 @@ def test_honors_projection(self):
self.assertEqual(ds.dims, {'time': 4248, 'lon': 3600, 'lat': 1800})
self.assertNotEqual(ds.dims, standard_ds.dims)

@absltest.skipIf(_SKIP_RASTERIO_TESTS, 'rioxarray module not loaded')
def test_expected_precise_transform(self):
data = np.empty((162, 121), dtype=np.float32)
bbox = (
-53.94158617595226,
-12.078281822698678,
-53.67209159071253,
-11.714464132625046,
)
x_res = (bbox[2] - bbox[0]) / data.shape[1]
y_res = (bbox[3] - bbox[1]) / data.shape[0]
raster = xr.DataArray(
data,
coords={
'y': np.linspace(bbox[3], bbox[1] + x_res, data.shape[0]),
'x': np.linspace(bbox[0], bbox[2] - y_res, data.shape[1]),
},
dims=('y', 'x'),
)
raster.rio.write_crs('EPSG:4326', inplace=True)
ic = (
ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
.filterDate(ee.DateRange('2014-01-01', '2014-01-02'))
.select('precipitation')
)
xee_dataset = xr.open_dataset(
ee.ImageCollection(ic),
engine='ee',
geometry=tuple(raster.rio.bounds()),
Copy link
Collaborator

Choose a reason for hiding this comment

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

Nice!

projection=ee.Projection(
crs=str(raster.rio.crs), transform=raster.rio.transform()[:6]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Also nice!

),
).rename({'lon': 'x', 'lat': 'y'})
self.assertNotEqual(abs(x_res), abs(y_res))
np.testing.assert_equal(
np.array(xee_dataset.rio.transform()),
np.array(raster.rio.transform()),
)

def test_parses_ee_url(self):
ds = self.entry.open_dataset(
'ee://LANDSAT/LC08/C01/T1',
Expand Down
Loading