one test not passing

test/core/gridmapping/test_base.py

@@ -336,7 +336,7 @@ def test_to_regular(self):
             pyproj.CRS.from_string("EPSG:32633"), transformed_gm_regular.crs
         )
         self.assertEqual((827, 1163), transformed_gm_regular.size)
-        self.assertEqual((827, 1163), transformed_gm_regular.tile_size)
+        self.assertEqual((1000, 1000), transformed_gm_regular.tile_size)
         self.assertEqual(False, transformed_gm_regular.is_j_axis_up)
         self.assertEqual(False, transformed_gm_regular.is_lon_360)
 

test/core/gridmapping/test_coords.py

@@ -300,7 +300,7 @@ def test_to_regular(self):
         gm_irr = GridMapping.from_coords(lon, lat, GEO_CRS)
         gm_reg_actual = gm_irr.to_regular()
         gm_reg_expected = GridMapping.regular(
-            size=(4, 4), xy_min=(-2, 48), xy_res=4.0, crs=GEO_CRS
+            size=(4, 4), tile_size=(2, 2), xy_min=(-2, 48), xy_res=4.0, crs=GEO_CRS
         )
         self.assertEqual(gm_reg_expected.size, gm_reg_actual.size)
         self.assertEqual(gm_reg_expected.tile_size, gm_reg_actual.tile_size)

test/core/resampling/test_rectify.py

@@ -3,7 +3,6 @@
 # https://opensource.org/licenses/MIT.
 
 import unittest
-from typing import Tuple
 
 import numpy as np
 import pytest

xcube/core/gridmapping/base.py

@@ -10,7 +10,6 @@
 from typing import Callable
 from collections.abc import Mapping
 from typing import Optional
-from typing import Tuple
 from typing import Union
 
 import numpy as np
@@ -696,6 +695,7 @@ def transform(
         self,
         crs: Union[str, pyproj.crs.CRS],
         *,
+        xy_res: Union[Number, tuple[Number, Number]] = None,
         tile_size: Union[int, tuple[int, int]] = None,
         xy_var_names: tuple[str, str] = None,
         tolerance: float = DEFAULT_TOLERANCE,
@@ -719,6 +719,7 @@ def transform(
         return transform_grid_mapping(
             self,
             crs,
+            xy_res=xy_res,
             tile_size=tile_size,
             xy_var_names=xy_var_names,
             tolerance=tolerance,

xcube/core/gridmapping/coords.py

@@ -19,6 +19,8 @@
 from .helpers import _default_xy_var_names
 from .helpers import _normalize_crs
 from .helpers import _normalize_int_pair
+from .helpers import _normalize_number_pair
+from .helpers import Number
 from .helpers import _to_int_or_float
 from .helpers import from_lon_360
 from .helpers import round_to_fraction
@@ -82,6 +84,7 @@ def new_grid_mapping_from_coords(
     y_coords: xr.DataArray,
     crs: Union[str, pyproj.crs.CRS],
     *,
+    xy_res: Union[Number, tuple[Number, Number]] = None,
     tile_size: Union[int, tuple[int, int]] = None,
     tolerance: float = DEFAULT_TOLERANCE,
 ) -> GridMapping:
@@ -169,62 +172,70 @@ def new_grid_mapping_from_coords(
         height, width = x_coords.shape
         size = width, height
 
-        x = da.asarray(x_coords)
-        y = da.asarray(y_coords)
-
-        x_x_diff = _abs_no_nan(da.diff(x[0, :]))
-        x_y_diff = _abs_no_nan(da.diff(x[:, 0]))
-        y_x_diff = _abs_no_nan(da.diff(y[0, :]))
-        y_y_diff = _abs_no_nan(da.diff(y[:, 0]))
-
-        if not is_lon_360 and crs.is_geographic:
-            is_anti_meridian_crossed = da.any(da.max(x_x_diff) > 180) or da.any(
-                da.max(x_y_diff) > 180
+        if xy_res is not None:
+            x_res, y_res = _normalize_number_pair(xy_res)
+        else:
+            x = da.asarray(x_coords)
+            y = da.asarray(y_coords)
+
+            x_x_diff = _abs_no_nan(da.diff(x, axis=1))
+            x_y_diff = _abs_no_nan(da.diff(x, axis=0))
+            y_x_diff = _abs_no_nan(da.diff(y, axis=1))
+            y_y_diff = _abs_no_nan(da.diff(y, axis=0))
+
+            if not is_lon_360 and crs.is_geographic:
+                is_anti_meridian_crossed = da.any(da.max(x_x_diff) > 180) or da.any(
+                    da.max(x_y_diff) > 180
+                )
+                if is_anti_meridian_crossed:
+                    x_coords = to_lon_360(x_coords)
+                    x = da.asarray(x_coords)
+                    x_x_diff = _abs_no_nan(da.diff(x, axis=1))
+                    x_y_diff = _abs_no_nan(da.diff(x, axis=0))
+                    is_lon_360 = True
+
+            x_res = x_x_diff[0, 0]
+            y_res = y_y_diff[0, 0]
+            is_regular = (
+                da.allclose(x_x_diff[0, :], x_res, atol=tolerance)
+                and da.allclose(x_x_diff[-1, :], x_res, atol=tolerance)
+                and da.allclose(y_y_diff[:, 0], y_res, atol=tolerance)
+                and da.allclose(y_y_diff[:, -1], y_res, atol=tolerance)
             )
-            if is_anti_meridian_crossed:
-                x_coords = to_lon_360(x_coords)
-                x = da.asarray(x_coords)
-                x_x_diff = _abs_no_nan(da.diff(x[0, :]))
-                x_y_diff = _abs_no_nan(da.diff(x[:, 0]))
-                is_lon_360 = True
-
-        x_res = x_x_diff[0]
-        y_res = y_y_diff[0]
-        is_regular = (
-            da.allclose(x_x_diff, x_res, atol=tolerance)
-            and da.allclose(x_x_diff, x_res, atol=tolerance)
-            and da.allclose(y_y_diff, y_res, atol=tolerance)
-            and da.allclose(y_y_diff, y_res, atol=tolerance)
-        )
-
-        if not is_regular:
-            # Find resolution via area
-            x_abs_diff = da.sqrt(da.square(x_x_diff) + da.square(x_y_diff))
-            y_abs_diff = da.sqrt(da.square(y_x_diff) + da.square(y_y_diff))
-            if crs.is_geographic:
-                # Convert degrees into meters
-                x_abs_diff_r = da.radians(x_abs_diff)
-                y_abs_diff_r = da.radians(y_abs_diff)
-                x_abs_diff = _ER * da.cos(x_abs_diff_r) * y_abs_diff_r
-                y_abs_diff = _ER * y_abs_diff_r
-            xy_areas = (x_abs_diff * y_abs_diff).flatten()
-            xy_areas = da.where(xy_areas > 0, xy_areas, np.nan)
-            # Get indices of min and max area
-            xy_area_index_min = da.nanargmin(xy_areas)
-            xy_area_index_max = da.nanargmax(xy_areas)
-            # Convert area to edge length
-            xy_res_min = math.sqrt(xy_areas[xy_area_index_min])
-            xy_res_max = math.sqrt(xy_areas[xy_area_index_max])
-            # Empirically weight min more than max
-            xy_res = 0.7 * xy_res_min + 0.3 * xy_res_max
-            if crs.is_geographic:
-                # Convert meters back into degrees
-                # print(f'xy_res in meters: {xy_res}')
-                xy_res = math.degrees(xy_res / _ER)
-                # print(f'xy_res in degrees: {xy_res}')
-            # Because this is an estimation, we can round to a nice number
-            xy_res = round_to_fraction(xy_res, digits=1, resolution=0.5)
-            x_res, y_res = float(xy_res), float(xy_res)
+            if not is_regular:
+                # Let diff arrays have same shape as original by
+                # doubling last rows and columns.
+                x_x_diff_c = da.concatenate([x_x_diff, x_x_diff[:, -1:]], axis=1)
+                y_x_diff_c = da.concatenate([y_x_diff, y_x_diff[:, -1:]], axis=1)
+                x_y_diff_c = da.concatenate([x_y_diff, x_y_diff[-1:, :]], axis=0)
+                y_y_diff_c = da.concatenate([y_y_diff, y_y_diff[-1:, :]], axis=0)
+                # Find resolution via area
+                x_abs_diff = da.sqrt(da.square(x_x_diff_c) + da.square(x_y_diff_c))
+                y_abs_diff = da.sqrt(da.square(y_x_diff_c) + da.square(y_y_diff_c))
+                if crs.is_geographic:
+                    # Convert degrees into meters
+                    x_abs_diff_r = da.radians(x_abs_diff)
+                    y_abs_diff_r = da.radians(y_abs_diff)
+                    x_abs_diff = _ER * da.cos(x_abs_diff_r) * y_abs_diff_r
+                    y_abs_diff = _ER * y_abs_diff_r
+                xy_areas = (x_abs_diff * y_abs_diff).flatten()
+                xy_areas = da.where(xy_areas > 0, xy_areas, np.nan)
+                # Get indices of min and max area
+                xy_area_index_min = da.nanargmin(xy_areas)
+                xy_area_index_max = da.nanargmax(xy_areas)
+                # Convert area to edge length
+                xy_res_min = math.sqrt(xy_areas[xy_area_index_min])
+                xy_res_max = math.sqrt(xy_areas[xy_area_index_max])
+                # Empirically weight min more than max
+                xy_res = 0.7 * xy_res_min + 0.3 * xy_res_max
+                if crs.is_geographic:
+                    # Convert meters back into degrees
+                    # print(f'xy_res in meters: {xy_res}')
+                    xy_res = math.degrees(xy_res / _ER)
+                    # print(f'xy_res in degrees: {xy_res}')
+                # Because this is an estimation, we can round to a nice number
+                xy_res = round_to_fraction(xy_res, digits=1, resolution=0.5)
+                x_res, y_res = float(xy_res), float(xy_res)
 
         if tile_size is None and x_coords.chunks is not None:
             j_chunks, i_chunks = x_coords.chunks

xcube/core/gridmapping/regular.py

@@ -114,6 +114,8 @@ def to_regular_grid_mapping(
     tile_size: Union[int, tuple[int, int]] = None,
     is_j_axis_up: bool = False,
 ) -> GridMapping:
+    if tile_size is None:
+        tile_size = grid_mapping.tile_size
     if grid_mapping.is_regular:
         if tile_size is not None or is_j_axis_up != grid_mapping.is_j_axis_up:
             return grid_mapping.derive(tile_size=tile_size, is_j_axis_up=is_j_axis_up)

xcube/core/gridmapping/transform.py

@@ -13,6 +13,7 @@
 from .base import GridMapping
 from .coords import new_grid_mapping_from_coords
 from .helpers import _assert_valid_xy_names
+from .helpers import Number
 from .helpers import _normalize_crs
 
 
@@ -40,6 +41,7 @@ def transform_grid_mapping(
     grid_mapping: GridMapping,
     crs: Union[str, pyproj.crs.CRS],
     *,
+    xy_res: Union[Number, tuple[Number, Number]] = None,
     tile_size: Union[int, tuple[int, int]] = None,
     xy_var_names: tuple[str, str] = None,
     tolerance: float = DEFAULT_TOLERANCE,
@@ -80,10 +82,14 @@ def _transform(block: np.ndarray) -> np.ndarray:
     #   on x and y arrays, this will take twice as long as if operation
     #   would be performed on the xy_coords dask array directly.
 
+    if tile_size is None:
+        tile_size = grid_mapping.tile_size
+
     return new_grid_mapping_from_coords(
         x_coords=xy_coords[0].rename(xy_var_names[0]),
         y_coords=xy_coords[1].rename(xy_var_names[1]),
         crs=target_crs,
+        xy_res=xy_res,
         tile_size=tile_size,
         tolerance=tolerance,
     )

xcube/core/resampling/spatial.py

@@ -196,7 +196,7 @@ def resample_in_space(
             )
             downscaled_dataset = resample_dataset(
                 source_ds,
-                ((x_scale, 1, 0), (1, y_scale, 0)),
+                ((1 / x_scale, 1, 0), (1, 1 / y_scale, 0)),
                 size=downscaled_size,
                 tile_size=source_gm.tile_size,
                 xy_dim_names=source_gm.xy_dim_names,
@@ -205,7 +205,7 @@ def resample_in_space(
             downscaled_gm = GridMapping.from_dataset(
                 downscaled_dataset,
                 tile_size=source_gm.tile_size,
-                prefer_crs=source_gm.crs,
+                crs=source_gm.crs,
             )
         return rectify_dataset(
             downscaled_dataset,