From ea7860484c8bf40042ab221fa4779aa44f883aa6 Mon Sep 17 00:00:00 2001 From: konstntokas Date: Tue, 22 Oct 2024 11:21:23 +0200 Subject: [PATCH] one test not passing --- test/core/gridmapping/test_base.py | 2 +- test/core/gridmapping/test_coords.py | 2 +- test/core/resampling/test_rectify.py | 1 - xcube/core/gridmapping/base.py | 3 +- xcube/core/gridmapping/coords.py | 121 +++++++++++++++------------ xcube/core/gridmapping/regular.py | 2 + xcube/core/gridmapping/transform.py | 6 ++ xcube/core/resampling/spatial.py | 4 +- 8 files changed, 80 insertions(+), 61 deletions(-) diff --git a/test/core/gridmapping/test_base.py b/test/core/gridmapping/test_base.py index 2d2efe0d2..49c2d8ec7 100644 --- a/test/core/gridmapping/test_base.py +++ b/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) diff --git a/test/core/gridmapping/test_coords.py b/test/core/gridmapping/test_coords.py index c468f9d42..3b9f60cb2 100644 --- a/test/core/gridmapping/test_coords.py +++ b/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) diff --git a/test/core/resampling/test_rectify.py b/test/core/resampling/test_rectify.py index 05de2e7de..99d39fee1 100644 --- a/test/core/resampling/test_rectify.py +++ b/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 diff --git a/xcube/core/gridmapping/base.py b/xcube/core/gridmapping/base.py index 68b9aa2ca..5419718c3 100644 --- a/xcube/core/gridmapping/base.py +++ b/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, diff --git a/xcube/core/gridmapping/coords.py b/xcube/core/gridmapping/coords.py index cf8e0e81a..75e519496 100644 --- a/xcube/core/gridmapping/coords.py +++ b/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 diff --git a/xcube/core/gridmapping/regular.py b/xcube/core/gridmapping/regular.py index 21796b2e7..24a03a181 100644 --- a/xcube/core/gridmapping/regular.py +++ b/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) diff --git a/xcube/core/gridmapping/transform.py b/xcube/core/gridmapping/transform.py index 32fd72006..11eea404a 100644 --- a/xcube/core/gridmapping/transform.py +++ b/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, ) diff --git a/xcube/core/resampling/spatial.py b/xcube/core/resampling/spatial.py index 83dd0fecc..c95a64298 100644 --- a/xcube/core/resampling/spatial.py +++ b/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,