From f1b63cfb5e4094590424100193eac05a7dfd3596 Mon Sep 17 00:00:00 2001 From: Rylan Boothman Date: Thu, 9 Nov 2023 17:27:51 -0500 Subject: [PATCH 1/7] transform dataset coordinates into specified crs --- setup.py | 2 + xee/ext.py | 12 +++++- xee/ext_integration_test.py | 74 +++++++++++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index ab6ce73..8af797d 100644 --- a/setup.py +++ b/setup.py @@ -28,6 +28,8 @@ "absl-py", "pytest", "pyink", + "rioxarray", + "rasterio", ] examples_require = [ diff --git a/xee/ext.py b/xee/ext.py index 99a826f..db5ff66 100644 --- a/xee/ext.py +++ b/xee/ext.py @@ -253,6 +253,12 @@ def __init__( else: self.mask_value = mask_value + self.coordinate_transformer = pyproj.Transformer.from_crs( + "EPSG:4326", + self.crs, + always_xy=True + ) + @functools.cached_property def get_info(self) -> dict[str, Any]: """Make all getInfo() calls to EE at once.""" @@ -555,11 +561,13 @@ def get_variables(self) -> utils.Frozen[str, xarray.Variable]: lon_grid = self.project((0, 0, v0.shape[1], 1)) lat_grid = self.project((0, 0, 1, v0.shape[2])) lon = self.image_to_array( - lnglat_img, grid=lon_grid, dtype=np.float32, bandIds=['longitude'] + lnglat_img, grid=lon_grid, dtype=np.float32, bandIds=['longitude', 'latitude'] ) + lon = self.coordinate_transformer.transform(lon[0], lon[1])[0] lat = self.image_to_array( - lnglat_img, grid=lat_grid, dtype=np.float32, bandIds=['latitude'] + lnglat_img, grid=lat_grid, dtype=np.float32, bandIds=['longitude', 'latitude'] ) + lat = self.coordinate_transformer.transform(lat[0], lat[1])[1] width_coord = np.squeeze(lon) height_coord = np.squeeze(lat) diff --git a/xee/ext_integration_test.py b/xee/ext_integration_test.py index e340f52..7b4b95c 100644 --- a/xee/ext_integration_test.py +++ b/xee/ext_integration_test.py @@ -22,6 +22,10 @@ import numpy as np import xarray as xr from xarray.core import indexing +import os +import rioxarray +import rasterio +import tempfile import xee import ee @@ -397,6 +401,76 @@ def test_validate_band_attrs(self): for _, value in variable.attrs.items(): self.assertIsInstance(value, valid_types) + def test_write_projected_dataset_to_raster(self): + # ensure that a projected dataset written to a raster intersects with the + # point used to create the initial image collection + with tempfile.TemporaryDirectory() as temp_dir: + temp_file = os.path.join(temp_dir, "test.tif") + + crs = "epsg:32610" + proj = ee.Projection(crs) + point = ee.Geometry.Point([-122.44, 37.78]) + geom = point.buffer(1024).bounds() + + col = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") + col = col.filterBounds(point) + col = col.filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", 5)) + col = col.limit(10) + + ds = xr.open_dataset( + col, + engine=xee.EarthEngineBackendEntrypoint, + scale=10, + crs=crs, + geometry=geom, + ) + + ds = ds.isel(time=0).transpose("Y", "X") + ds.rio.set_spatial_dims(x_dim="X", y_dim="Y", inplace=True) + ds.rio.write_crs(crs, inplace=True) + ds.rio.reproject(crs, inplace=True) + ds.rio.to_raster(temp_file) + + with rasterio.open(temp_file) as raster: + # see https://gis.stackexchange.com/a/407755 for evenOdd explanation + bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) + intersects = bbox.intersects(point, 1, proj=proj) + self.assertTrue(intersects.getInfo()) + + def test_write_dataset_to_raster(self): + # ensure that a dataset written to a raster intersects with the point used + # to create the initial image collection + with tempfile.TemporaryDirectory() as temp_dir: + temp_file = os.path.join(temp_dir, "test.tif") + + crs = "EPSG:4326" + proj = ee.Projection(crs) + point = ee.Geometry.Point([-122.44, 37.78]) + geom = point.buffer(1024).bounds() + + col = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") + col = col.filterBounds(point) + col = col.filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", 5)) + col = col.limit(10) + + ds = xr.open_dataset( + col, + engine=xee.EarthEngineBackendEntrypoint, + scale=0.0025, + geometry=geom, + ) + + ds = ds.isel(time=0).transpose("lat", "lon") + ds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True) + ds.rio.write_crs(crs, inplace=True) + ds.rio.reproject(crs, inplace=True) + ds.rio.to_raster(temp_file) + + with rasterio.open(temp_file) as raster: + # see https://gis.stackexchange.com/a/407755 for evenOdd explanation + bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) + intersects = bbox.intersects(point, 1, proj=proj) + self.assertTrue(intersects.getInfo()) if __name__ == '__main__': absltest.main() From d698f43cc1c6e7f2d0175bab1626efa2beab6ae9 Mon Sep 17 00:00:00 2001 From: Rylan Boothman Date: Mon, 13 Nov 2023 08:24:51 -0500 Subject: [PATCH 2/7] pyink linting --- xee/ext.py | 14 +++-- xee/ext_integration_test.py | 121 ++++++++++++++++++------------------ 2 files changed, 70 insertions(+), 65 deletions(-) diff --git a/xee/ext.py b/xee/ext.py index db5ff66..08133e8 100644 --- a/xee/ext.py +++ b/xee/ext.py @@ -254,9 +254,7 @@ def __init__( self.mask_value = mask_value self.coordinate_transformer = pyproj.Transformer.from_crs( - "EPSG:4326", - self.crs, - always_xy=True + 'EPSG:4326', self.crs, always_xy=True ) @functools.cached_property @@ -561,11 +559,17 @@ def get_variables(self) -> utils.Frozen[str, xarray.Variable]: lon_grid = self.project((0, 0, v0.shape[1], 1)) lat_grid = self.project((0, 0, 1, v0.shape[2])) lon = self.image_to_array( - lnglat_img, grid=lon_grid, dtype=np.float32, bandIds=['longitude', 'latitude'] + lnglat_img, + grid=lon_grid, + dtype=np.float32, + bandIds=['longitude', 'latitude'], ) lon = self.coordinate_transformer.transform(lon[0], lon[1])[0] lat = self.image_to_array( - lnglat_img, grid=lat_grid, dtype=np.float32, bandIds=['longitude', 'latitude'] + lnglat_img, + grid=lat_grid, + dtype=np.float32, + bandIds=['longitude', 'latitude'], ) lat = self.coordinate_transformer.transform(lat[0], lat[1])[1] width_coord = np.squeeze(lon) diff --git a/xee/ext_integration_test.py b/xee/ext_integration_test.py index 7b4b95c..a5995c9 100644 --- a/xee/ext_integration_test.py +++ b/xee/ext_integration_test.py @@ -405,72 +405,73 @@ def test_write_projected_dataset_to_raster(self): # ensure that a projected dataset written to a raster intersects with the # point used to create the initial image collection with tempfile.TemporaryDirectory() as temp_dir: - temp_file = os.path.join(temp_dir, "test.tif") - - crs = "epsg:32610" - proj = ee.Projection(crs) - point = ee.Geometry.Point([-122.44, 37.78]) - geom = point.buffer(1024).bounds() - - col = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") - col = col.filterBounds(point) - col = col.filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", 5)) - col = col.limit(10) - - ds = xr.open_dataset( - col, - engine=xee.EarthEngineBackendEntrypoint, - scale=10, - crs=crs, - geometry=geom, - ) - - ds = ds.isel(time=0).transpose("Y", "X") - ds.rio.set_spatial_dims(x_dim="X", y_dim="Y", inplace=True) - ds.rio.write_crs(crs, inplace=True) - ds.rio.reproject(crs, inplace=True) - ds.rio.to_raster(temp_file) - - with rasterio.open(temp_file) as raster: - # see https://gis.stackexchange.com/a/407755 for evenOdd explanation - bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) - intersects = bbox.intersects(point, 1, proj=proj) - self.assertTrue(intersects.getInfo()) + temp_file = os.path.join(temp_dir, 'test.tif') + + crs = 'epsg:32610' + proj = ee.Projection(crs) + point = ee.Geometry.Point([-122.44, 37.78]) + geom = point.buffer(1024).bounds() + + col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') + col = col.filterBounds(point) + col = col.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 5)) + col = col.limit(10) + + ds = xr.open_dataset( + col, + engine=xee.EarthEngineBackendEntrypoint, + scale=10, + crs=crs, + geometry=geom, + ) + + ds = ds.isel(time=0).transpose('Y', 'X') + ds.rio.set_spatial_dims(x_dim='X', y_dim='Y', inplace=True) + ds.rio.write_crs(crs, inplace=True) + ds.rio.reproject(crs, inplace=True) + ds.rio.to_raster(temp_file) + + with rasterio.open(temp_file) as raster: + # see https://gis.stackexchange.com/a/407755 for evenOdd explanation + bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) + intersects = bbox.intersects(point, 1, proj=proj) + self.assertTrue(intersects.getInfo()) def test_write_dataset_to_raster(self): # ensure that a dataset written to a raster intersects with the point used # to create the initial image collection with tempfile.TemporaryDirectory() as temp_dir: - temp_file = os.path.join(temp_dir, "test.tif") - - crs = "EPSG:4326" - proj = ee.Projection(crs) - point = ee.Geometry.Point([-122.44, 37.78]) - geom = point.buffer(1024).bounds() - - col = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") - col = col.filterBounds(point) - col = col.filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", 5)) - col = col.limit(10) - - ds = xr.open_dataset( - col, - engine=xee.EarthEngineBackendEntrypoint, - scale=0.0025, - geometry=geom, - ) + temp_file = os.path.join(temp_dir, 'test.tif') + + crs = 'EPSG:4326' + proj = ee.Projection(crs) + point = ee.Geometry.Point([-122.44, 37.78]) + geom = point.buffer(1024).bounds() + + col = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') + col = col.filterBounds(point) + col = col.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 5)) + col = col.limit(10) + + ds = xr.open_dataset( + col, + engine=xee.EarthEngineBackendEntrypoint, + scale=0.0025, + geometry=geom, + ) + + ds = ds.isel(time=0).transpose('lat', 'lon') + ds.rio.set_spatial_dims(x_dim='lon', y_dim='lat', inplace=True) + ds.rio.write_crs(crs, inplace=True) + ds.rio.reproject(crs, inplace=True) + ds.rio.to_raster(temp_file) + + with rasterio.open(temp_file) as raster: + # see https://gis.stackexchange.com/a/407755 for evenOdd explanation + bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) + intersects = bbox.intersects(point, 1, proj=proj) + self.assertTrue(intersects.getInfo()) - ds = ds.isel(time=0).transpose("lat", "lon") - ds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True) - ds.rio.write_crs(crs, inplace=True) - ds.rio.reproject(crs, inplace=True) - ds.rio.to_raster(temp_file) - - with rasterio.open(temp_file) as raster: - # see https://gis.stackexchange.com/a/407755 for evenOdd explanation - bbox = ee.Geometry.Rectangle(raster.bounds, proj=proj, evenOdd=False) - intersects = bbox.intersects(point, 1, proj=proj) - self.assertTrue(intersects.getInfo()) if __name__ == '__main__': absltest.main() From e2606594754e8ba20081a9c99927d37645c66f49 Mon Sep 17 00:00:00 2001 From: Rylan Boothman Date: Fri, 17 Nov 2023 17:12:09 -0500 Subject: [PATCH 3/7] add rasterio and rioxarray to test dependencies --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 942cb56..81a5fcb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,8 @@ tests = [ "absl-py", "pytest", "pyink", + "rasterio", + "rioxarray", ] examples = [ "apache_beam[gcp]", From c6d574d675d864ba051ef7effc51715a88c96d31 Mon Sep 17 00:00:00 2001 From: Rylan Boothman Date: Fri, 17 Nov 2023 17:12:35 -0500 Subject: [PATCH 4/7] fix imports --- xee/ext_integration_test.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/xee/ext_integration_test.py b/xee/ext_integration_test.py index a5995c9..b622a78 100644 --- a/xee/ext_integration_test.py +++ b/xee/ext_integration_test.py @@ -16,16 +16,15 @@ import json import os import pathlib +import tempfile from absl.testing import absltest from google.auth import identity_pool import numpy as np import xarray as xr from xarray.core import indexing -import os import rioxarray import rasterio -import tempfile import xee import ee From 00ba44476bb5141a15d54f074f7f565bce677a01 Mon Sep 17 00:00:00 2001 From: Rylan Boothman Date: Mon, 11 Dec 2023 13:51:32 -0500 Subject: [PATCH 5/7] use pixelCoordinates instead of pixelLonLat --- xee/ext.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/xee/ext.py b/xee/ext.py index b98e394..2afa104 100644 --- a/xee/ext.py +++ b/xee/ext.py @@ -253,10 +253,6 @@ def __init__( else: self.mask_value = mask_value - self.coordinate_transformer = pyproj.Transformer.from_crs( - 'EPSG:4326', self.crs, always_xy=True - ) - @functools.cached_property def get_info(self) -> Dict[str, Any]: """Make all getInfo() calls to EE at once.""" @@ -555,23 +551,15 @@ def get_variables(self) -> utils.Frozen[str, xarray.Variable]: f'ImageCollection due to: {e}.' ) - lnglat_img = ee.Image.pixelLonLat() + lnglat_img = ee.Image.pixelCoordinates(ee.Projection(self.crs_arg)) lon_grid = self.project((0, 0, v0.shape[1], 1)) lat_grid = self.project((0, 0, 1, v0.shape[2])) lon = self.image_to_array( - lnglat_img, - grid=lon_grid, - dtype=np.float32, - bandIds=['longitude', 'latitude'], + lnglat_img, grid=lon_grid, dtype=np.float32, bandIds=['x'] ) - lon = self.coordinate_transformer.transform(lon[0], lon[1])[0] lat = self.image_to_array( - lnglat_img, - grid=lat_grid, - dtype=np.float32, - bandIds=['longitude', 'latitude'], + lnglat_img, grid=lat_grid, dtype=np.float32, bandIds=['y'] ) - lat = self.coordinate_transformer.transform(lat[0], lat[1])[1] width_coord = np.squeeze(lon) height_coord = np.squeeze(lat) From 81766cf344a94a2f85df328adc728d5487c5532e Mon Sep 17 00:00:00 2001 From: Rahul Mahrsee Date: Wed, 20 Dec 2023 10:43:31 +0000 Subject: [PATCH 6/7] lint fixes. --- xee/ext.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/xee/ext.py b/xee/ext.py index 32c13a2..1618191 100644 --- a/xee/ext.py +++ b/xee/ext.py @@ -573,9 +573,7 @@ def _process_coordinate_data( self._get_tile_from_ee, list(zip(data, itertools.cycle([coordinate_type]))), ): - tiles[i] = ( - arr.tolist() if coordinate_type == 'x' else arr.tolist()[0] - ) + tiles[i] = arr.tolist() if coordinate_type == 'x' else arr.tolist()[0] return np.concatenate(tiles) def get_variables(self) -> utils.Frozen[str, xarray.Variable]: @@ -609,7 +607,6 @@ def get_variables(self) -> utils.Frozen[str, xarray.Variable]: lat_total_tile = math.ceil(v0.shape[2] / height_chunk) lat = self._process_coordinate_data( lat_total_tile, height_chunk, v0.shape[2], 'y' - ) width_coord = np.squeeze(lon) From 24ff5825558f01fb6bc1fc2395192bcecb554e32 Mon Sep 17 00:00:00 2001 From: Nathaniel Schmitz Date: Thu, 21 Dec 2023 11:18:29 -0500 Subject: [PATCH 7/7] Conditionally run rasterio unit tests. --- xee/ext_integration_test.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/xee/ext_integration_test.py b/xee/ext_integration_test.py index 35475d2..d1eafc0 100644 --- a/xee/ext_integration_test.py +++ b/xee/ext_integration_test.py @@ -23,12 +23,17 @@ import numpy as np import xarray as xr from xarray.core import indexing -import rioxarray -import rasterio import xee import ee +_SKIP_RASTERIO_TESTS = False +try: + import rasterio # pylint: disable=g-import-not-at-top + import rioxarray # pylint: disable=g-import-not-at-top,unused-import +except ImportError: + _SKIP_RASTERIO_TESTS = True + _CREDENTIALS_PATH_KEY = 'GOOGLE_APPLICATION_CREDENTIALS' _SCOPES = [ 'https://www.googleapis.com/auth/cloud-platform', @@ -400,6 +405,7 @@ def test_validate_band_attrs(self): for _, value in variable.attrs.items(): self.assertIsInstance(value, valid_types) + @absltest.skipIf(_SKIP_RASTERIO_TESTS, 'rioxarray module not loaded') def test_write_projected_dataset_to_raster(self): # ensure that a projected dataset written to a raster intersects with the # point used to create the initial image collection @@ -436,6 +442,7 @@ def test_write_projected_dataset_to_raster(self): intersects = bbox.intersects(point, 1, proj=proj) self.assertTrue(intersects.getInfo()) + @absltest.skipIf(_SKIP_RASTERIO_TESTS, 'rioxarray module not loaded') def test_write_dataset_to_raster(self): # ensure that a dataset written to a raster intersects with the point used # to create the initial image collection