From e11d8781d497c9f8ca7e8232947e7973b2c09a98 Mon Sep 17 00:00:00 2001 From: Stephen Po-Chedley Date: Mon, 5 Jun 2023 15:10:57 -0700 Subject: [PATCH 1/3] fix for #494 --- tests/test_spatial.py | 21 +++++---------- xcdat/spatial.py | 61 ++++++++++++++++--------------------------- 2 files changed, 29 insertions(+), 53 deletions(-) diff --git a/tests/test_spatial.py b/tests/test_spatial.py index 1a486c5b..4e13c9c7 100644 --- a/tests/test_spatial.py +++ b/tests/test_spatial.py @@ -269,26 +269,17 @@ def setup(self): decode_times=True, cf_compliant=False, has_bounds=True ) - def test_bounds_reordered_when_upper_indexed_first(self): + def test_value_error_thrown_for_multiple_out_of_order_lon_bounds(self): domain_bounds = xr.DataArray( name="lon_bnds", - data=np.array( - [[-89.375, -90], [0.0, -89.375], [0.0, 89.375], [89.375, 90]] - ), + data=np.array([[3, 1], [5, 3], [5, 7], [7, 9]]), coords={"lat": self.ds.lat}, dims=["lat", "bnds"], ) - result = self.ds.spatial._force_domain_order_low_to_high(domain_bounds) - - expected_domain_bounds = xr.DataArray( - name="lon_bnds", - data=np.array( - [[-90, -89.375], [-89.375, 0.0], [0.0, 89.375], [89.375, 90]] - ), - coords={"lat": self.ds.lat}, - dims=["lat", "bnds"], - ) - assert result.identical(expected_domain_bounds) + # Check _get_longitude_weights raises error when there are + # > 1 out-of-order bounds for the dataset. + with pytest.raises(ValueError): + self.ds.spatial._get_longitude_weights(domain_bounds, region_bounds=None) def test_raises_error_if_dataset_has_multiple_bounds_variables_for_an_axis(self): ds = self.ds.copy() diff --git a/xcdat/spatial.py b/xcdat/spatial.py index c51af518..186a4088 100644 --- a/xcdat/spatial.py +++ b/xcdat/spatial.py @@ -281,12 +281,6 @@ def get_weights( "to reference a specific data variable's axis bounds." ) - # The logic for generating longitude weights depends on the - # bounds being ordered such that d_bounds[:, 0] < d_bounds[:, 1]. - # They are re-ordered (if need be) for the purpose of creating - # weights. - d_bounds = self._force_domain_order_low_to_high(d_bounds) - r_bounds = axis_bounds[key]["region"] weights = axis_bounds[key]["weights_method"](d_bounds, r_bounds) @@ -324,33 +318,6 @@ def _validate_axis_arg(self, axis: List[SpatialAxis]): # Check the axis coordinate variable exists in the Dataset. get_dim_coords(self._dataset, key) - def _force_domain_order_low_to_high(self, domain_bounds: xr.DataArray): - """Reorders the ``domain_bounds`` low-to-high. - - This method ensures all lower bound values are less than the upper bound - values (``domain_bounds[:, 1] < domain_bounds[:, 1]``). - - Parameters - ---------- - domain_bounds: xr.DataArray - The bounds of an axis. - - Returns - ------ - xr.DataArray - The bounds of an axis (re-ordered if applicable). - """ - index_bad_cells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0] - - if len(index_bad_cells) > 0: - new_domain_bounds = domain_bounds.copy() - new_domain_bounds[index_bad_cells, 0] = domain_bounds[index_bad_cells, 1] - new_domain_bounds[index_bad_cells, 1] = domain_bounds[index_bad_cells, 0] - - return new_domain_bounds - - return domain_bounds - def _validate_region_bounds(self, axis: SpatialAxis, bounds: RegionAxisBounds): """Validates the ``bounds`` arg based on a set of criteria. @@ -441,12 +408,34 @@ def _get_longitude_weights( ------- xr.DataArray The longitude axis weights. + + Raises + ------ + ValueError + If the there are multiple instances in which the + domain_bounds[:, 0] > domain_bounds[:, 1] """ p_meridian_index: Optional[np.ndarray] = None d_bounds = domain_bounds.copy() + # check if there is more than one potential prime meridian cell + # there should only be one for rectilinear data + pmcells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0] + if len(pmcells) > 1: + raise ValueError( + "More than one longitude bound is out of order. Only one bound \n\ + spanning the prime meridian is permitted" + ) + # convert longitude bounds to 360 degree domain + d_bounds: xr.DataArray = self._swap_lon_axis(d_bounds, to=360) # type: ignore + # check for bounds spanning prime meridian + p_meridian_index = _get_prime_meridian_index(d_bounds) + # if bounds span a prime meridian, ensure bounds are organized to span 0-360 + # (without losing weight) + if p_meridian_index is not None: + d_bounds = _align_lon_bounds_to_360(d_bounds, p_meridian_index) + if region_bounds is not None: - d_bounds: xr.DataArray = self._swap_lon_axis(d_bounds, to=360) # type: ignore r_bounds: np.ndarray = self._swap_lon_axis( region_bounds, to=360 ) # type:ignore @@ -455,10 +444,6 @@ def _get_longitude_weights( if is_region_circular: r_bounds = np.array([0.0, 360.0]) - p_meridian_index = _get_prime_meridian_index(d_bounds) - if p_meridian_index is not None: - d_bounds = _align_lon_bounds_to_360(d_bounds, p_meridian_index) - d_bounds = self._scale_domain_to_region(d_bounds, r_bounds) weights = self._calculate_weights(d_bounds) From 398855015eeea2a3110494cf8c484bbf90c75e37 Mon Sep 17 00:00:00 2001 From: Stephen Po-Chedley Date: Wed, 7 Jun 2023 20:50:35 -0700 Subject: [PATCH 2/3] add test for out of order bounds --- tests/test_spatial.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_spatial.py b/tests/test_spatial.py index 4e13c9c7..fe0361cd 100644 --- a/tests/test_spatial.py +++ b/tests/test_spatial.py @@ -180,6 +180,24 @@ def test_spatial_average_for_lat_region(self): assert result.identical(expected) + def test_spatial_average_for_domain_wrapping_p_meridian_non_cf_conventions( + self, + ): + ds = self.ds.copy() + + # get spatial average for original dataset + ref = ds.spatial.average("ts").ts + + # change first bound from -0.9375 to 359.0625 + lon_bnds = ds.lon_bnds.copy() + lon_bnds[0, 0] = 359.0625 + ds["lon_bnds"] = lon_bnds + + # check spatial average with new (bad) bound + result = ds.spatial.average("ts").ts + + assert result.identical(ref) + @requires_dask def test_spatial_average_for_lat_region_and_keep_weights_with_dask(self): ds = self.ds.copy().chunk(2) From 852cf430469811be6ca2304b69c99361f0c77b96 Mon Sep 17 00:00:00 2001 From: Stephen Po-Chedley Date: Tue, 13 Jun 2023 12:29:18 -0700 Subject: [PATCH 3/3] Apply suggestions from code review Co-authored-by: Tom Vo --- xcdat/spatial.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/xcdat/spatial.py b/xcdat/spatial.py index 186a4088..2c50595a 100644 --- a/xcdat/spatial.py +++ b/xcdat/spatial.py @@ -418,20 +418,15 @@ def _get_longitude_weights( p_meridian_index: Optional[np.ndarray] = None d_bounds = domain_bounds.copy() - # check if there is more than one potential prime meridian cell - # there should only be one for rectilinear data - pmcells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0] - if len(pmcells) > 1: + pm_cells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0] + if len(pm_cells) > 1: raise ValueError( - "More than one longitude bound is out of order. Only one bound \n\ - spanning the prime meridian is permitted" + "More than one longitude bound is out of order. Only one bound " + "value spanning the prime meridian is permitted in data on " + "a rectilinear grid." ) - # convert longitude bounds to 360 degree domain d_bounds: xr.DataArray = self._swap_lon_axis(d_bounds, to=360) # type: ignore - # check for bounds spanning prime meridian p_meridian_index = _get_prime_meridian_index(d_bounds) - # if bounds span a prime meridian, ensure bounds are organized to span 0-360 - # (without losing weight) if p_meridian_index is not None: d_bounds = _align_lon_bounds_to_360(d_bounds, p_meridian_index)