diff --git a/xee/ext.py b/xee/ext.py index c54b1cc..9ca80f5 100644 --- a/xee/ext.py +++ b/xee/ext.py @@ -240,7 +240,7 @@ def __init__( default_scale = self.SCALE_UNITS.get(self.scale_units, 1) if scale is None: scale = default_scale - default_transform = affine.Affine.scale(scale, -1 * scale) + default_transform = affine.Affine.scale(scale, scale) transform = affine.Affine(*proj.get('transform', default_transform)[:6]) self.scale_x, self.scale_y = transform.a, transform.e @@ -420,13 +420,25 @@ def project(self, bbox: types.BBox) -> types.Grid: appropriate region of data to return according to the Array's configured projection and scale. """ - # The origin of the image is in the top left corner. X is the minimum value - # and Y is the maximum value. - x_origin, _, _, y_origin = self.bounds # x_min, x_max, y_min, y_max + x_min, y_min, x_max, y_max = self.bounds x_start, y_start, x_end, y_end = bbox width = x_end - x_start height = y_end - y_start + # Find the actual coordinates of the first or last point of the bounding box + # (bbox) based on the positive and negative scale in the actual Earth Engine + # (EE) image. Since EE bounding boxes can be flipped (negative scale), we + # cannot determine the origin (transform translation) simply by identifying + # the min and max extents. Instead, we calculate the translation by + # considering the direction of scaling (positive or negative) along both + # the x and y axes. + translate_x = self.scale_x * x_start + ( + x_min if self.scale_x > 0 else x_max + ) + translate_y = self.scale_y * y_start + ( + y_min if self.scale_y > 0 else y_max + ) + return { # The size of the bounding box. The affine transform and project will be # applied, so we can think in terms of pixels. @@ -435,11 +447,8 @@ def project(self, bbox: types.BBox) -> types.Grid: 'height': height, }, 'affineTransform': { - # Since the origin is in the top left corner, we want to translate - # the start of the grid to the positive direction for X and the - # negative direction for Y. - 'translateX': x_origin + self.scale_x * x_start, - 'translateY': y_origin + self.scale_y * y_start, + 'translateX': translate_x, + 'translateY': translate_y, # Define the scale for each pixel (e.g. the number of meters between # each value). 'scaleX': self.scale_x, diff --git a/xee/ext_integration_test.py b/xee/ext_integration_test.py index 4b6a35f..e2d7466 100644 --- a/xee/ext_integration_test.py +++ b/xee/ext_integration_test.py @@ -338,8 +338,8 @@ def test_open_dataset__sanity_check(self): ds = self.entry.open_dataset( pathlib.Path('LANDSAT') / 'LC08' / 'C01' / 'T1', drop_variables=tuple(f'B{i}' for i in range(3, 12)), - scale=25.0, # in degrees n_images=3, + projection=ee.Projection('EPSG:4326', [25, 0, 0, 0, -25, 0]), ) self.assertEqual(dict(ds.dims), {'time': 3, 'lon': 14, 'lat': 7}) self.assertNotEmpty(dict(ds.coords)) @@ -352,6 +352,24 @@ def test_open_dataset__sanity_check(self): self.assertFalse(v.isnull().all(), 'All values are null!') self.assertEqual(v.shape, (3, 14, 7)) + def test_open_dataset__sanity_check_with_negative_scale(self): + ds = self.entry.open_dataset( + pathlib.Path('LANDSAT') / 'LC08' / 'C01' / 'T1', + drop_variables=tuple(f'B{i}' for i in range(3, 12)), + scale=-25.0, # in degrees + n_images=3, + ) + self.assertEqual(dict(ds.dims), {'time': 3, 'lon': 14, 'lat': 7}) + self.assertNotEmpty(dict(ds.coords)) + self.assertEqual( + list(ds.data_vars.keys()), + [f'B{i}' for i in range(1, 3)] + ['BQA'], + ) + for v in ds.values(): + self.assertIsNotNone(v.data) + self.assertTrue(v.isnull().all(), 'All values must be null!') + self.assertEqual(v.shape, (3, 14, 7)) + def test_open_dataset__n_images(self): ds = self.entry.open_dataset( pathlib.Path('LANDSAT') / 'LC08' / 'C01' / 'T1', @@ -516,9 +534,9 @@ def test_write_projected_dataset_to_raster(self): ds = xr.open_dataset( col, engine=xee.EarthEngineBackendEntrypoint, - scale=10, crs=crs, geometry=geom, + projection=ee.Projection('EPSG:4326', [10, 0, 0, 0, -10, 0]), ) ds = ds.isel(time=0).transpose('Y', 'X') @@ -553,8 +571,8 @@ def test_write_dataset_to_raster(self): ds = xr.open_dataset( col, engine=xee.EarthEngineBackendEntrypoint, - scale=0.0025, geometry=geom, + projection=ee.Projection('EPSG:4326', [0.0025, 0, 0, 0, -0.0025, 0]), ) ds = ds.isel(time=0).transpose('lat', 'lon')