Minor adjustments and bug fixing of resampling_in_space

CHANGES.md

@@ -55,6 +55,10 @@
 
 ### Fixes
 
+* The function `xcube.core.resample.resample_in_space()` now supports the parameter
+  `source_ds_subset=True` when calling `rectify_dataset`. This feature enables
+  performing the reprojection exclusively on the spatially congruent subset of 
+  the dataset.
 * The function `xcube.core.resample.resample_in_space()` now always operates
   lazily and therefore supports chunk-wise, parallel processing. (#1082)
 * Bux fix in the `has_data` method of the `"https"` data store

test/core/resampling/test_spatial.py

@@ -85,3 +85,84 @@ def test_rectify_and_upscale_dataset(self):
                 dtype=target_ds.rad.dtype,
             ),
         )
+
+    def test_reproject_utm(self):
+        source_ds = self.new_5x5_dataset_regular_utm()
+
+        # test projected CRS similar resolution
+        target_gm = GridMapping.regular(
+            size=(5, 5), xy_min=(4320080, 3382480), xy_res=80, crs="epsg:3035"
+        )
+        target_ds = resample_in_space(source_ds, target_gm=target_gm)
+        np.testing.assert_almost_equal(
+            target_ds.band_1.values,
+            np.array(
+                [
+                    [1, 1, 2, 3, 4],
+                    [6, 6, 7, 8, 9],
+                    [11, 12, 12, 13, 14],
+                    [16, 17, 17, 18, 19],
+                    [21, 17, 17, 18, 19],
+                ],
+                dtype=target_ds.band_1.dtype,
+            ),
+        )
+
+        # test projected CRS finer resolution
+        # test if subset calculation works as expected
+        target_gm = GridMapping.regular(
+            size=(5, 5), xy_min=(4320080, 3382480), xy_res=20, crs="epsg:3035"
+        )
+        target_ds = resample_in_space(source_ds, target_gm=target_gm)
+        np.testing.assert_almost_equal(
+            target_ds.band_1.values,
+            np.array(
+                [
+                    [15, 16, 16, 16, 16],
+                    [15, 16, 16, 16, 16],
+                    [15, 16, 16, 16, 16],
+                    [20, 21, 21, 21, 21],
+                    [20, 21, 21, 21, 21],
+                ],
+                dtype=target_ds.band_1.dtype,
+            ),
+        )
+
+        # test geographic CRS with similar resolution
+        target_gm = GridMapping.regular(
+            size=(5, 5), xy_min=(9.9886, 53.5499), xy_res=0.0006, crs=CRS_WGS84
+        )
+        target_ds = resample_in_space(source_ds, target_gm=target_gm)
+        np.testing.assert_almost_equal(
+            target_ds.band_1.values,
+            np.array(
+                [
+                    [7, 8, 8, 8, 9],
+                    [12, 13, 13, 13, 14],
+                    [12, 13, 13, 13, 14],
+                    [17, 18, 18, 18, 19],
+                    [22, 23, 23, 23, 24],
+                ],
+                dtype=target_ds.band_1.dtype,
+            ),
+        )
+
+        # test geographic CRS with 1/2 resolution
+        # test if subset calculation works as expected
+        target_gm = GridMapping.regular(
+            size=(5, 5), xy_min=(9.9886, 53.5499), xy_res=0.0003, crs=CRS_WGS84
+        )
+        target_ds = resample_in_space(source_ds, target_gm=target_gm)
+        np.testing.assert_almost_equal(
+            target_ds.band_1.values,
+            np.array(
+                [
+                    [12, 12, 12, 13, 13],
+                    [17, 17, 17, 18, 18],
+                    [17, 17, 17, 18, 18],
+                    [22, 17, 17, 18, 18],
+                    [22, 22, 22, 23, 23],
+                ],
+                dtype=target_ds.band_1.dtype,
+            ),
+        )

test/sampledata.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import pandas as pd
+import pyproj
 import xarray as xr
 
 
@@ -368,6 +369,23 @@ def new_2x2_dataset_with_irregular_coords(cls):
             )
         )
 
+    @classmethod
+    def new_5x5_dataset_regular_utm(cls):
+        x = np.arange(565300.0, 565800.0, 100.0)
+        y = np.arange(5934300.0, 5933800.0, -100.0)
+        spatial_ref = np.array(0)
+        band_1 = np.arange(25).reshape((5, 5))
+        ds = xr.Dataset(
+            dict(
+                band_1=xr.DataArray(
+                    band_1, dims=("y", "x"), attrs=dict(grid_mapping="spatial_ref")
+                )
+            ),
+            coords=dict(x=x, y=y, spatial_ref=spatial_ref),
+        )
+        ds.spatial_ref.attrs = pyproj.CRS.from_epsg("32632").to_cf()
+        return ds
+
     @classmethod
     def new_2x2_dataset_with_irregular_coords_antimeridian(cls):
         lon = np.array([[+179.0, -176.0], [+178.0, +180.0]])

xcube/core/gridmapping/cfconv.py

@@ -104,7 +104,6 @@ def get_dataset_grid_mapping_proxies(
 
     # Find coordinate variables.
     #
-
     latitude_longitude_coords = GridCoords()
     rotated_latitude_longitude_coords = GridCoords()
     projected_coords = GridCoords()

xcube/core/resampling/spatial.py

@@ -208,11 +208,34 @@
     # If CRSes are not both geographic and their CRSes are different
     # transform the source_gm so its CRS matches the target CRS:
     transformed_source_gm = source_gm.transform(target_gm.crs, xy_res=target_gm.xy_res)
+    source_ds = source_ds.drop_vars(source_gm.xy_dim_names)
+    for var in source_ds.data_vars:
+        if "grid_mapping" in source_ds[var].attrs:
+            attrs = source_ds[var].attrs
+            source_ds = source_ds.drop_vars(attrs["grid_mapping"])
+            del attrs["grid_mapping"]
+            source_ds[var] = source_ds[var].assign_attrs(attrs)
+    if "crs" in source_ds:
+        source_ds = source_ds.drop_vars("crs")
+    if "spatial_ref" in source_ds:
+        source_ds = source_ds.drop_vars("spatial_ref")
+    source_ds = source_ds.copy()
     transformed_x, transformed_y = transformed_source_gm.xy_coords
+    attrs = dict(grid_mapping="spatial_ref")
+    transformed_x.attrs = attrs
+    transformed_y.attrs = attrs
+    source_ds = source_ds.assign_coords(
+        spatial_ref=xr.DataArray(0, attrs=transformed_source_gm.crs.to_cf()),
+        transformed_x=transformed_x,
+        transformed_y=transformed_y,
+    )
     return resample_in_space(
-        source_ds.assign(transformed_x=transformed_x, transformed_y=transformed_y),
+        source_ds,
         source_gm=transformed_source_gm,
         ref_ds=ref_ds,
         target_gm=target_gm,
+        var_configs=var_configs,
+        encode_cf=encode_cf,
         gm_name=gm_name,
+        rectify_kwargs=rectify_kwargs,
     )