Fix spatial bugs: handle datasets with domain bounds out of order and zonal averaging #340

pochedls · 2022-09-10T17:44:54Z

Description

This PR is intended to fix #324 and #338. It allows users to perform spatial averaging on datasets in which domain_bounds[:, 0] > domain_bounds[:, 1] by re-ordering the domain bounds for the purposes of calculating spatial weights. It also fixes an issue that prevented spatial averaging on lat-plev grids. A section of the code tried to load longitude bounds even when that axis was not specified for averaging; this section of code is updated to only load bounds information if the axis is specified for averaging.

Closes [Bug]: averaging a lat-plev grid fails #324
Closes [Feature]: Handle "out of order" bounds #338

Checklist

My code follows the style guidelines of this project
I have performed a self-review of my own code
My changes generate no new warnings
Any dependent changes have been merged and published in downstream modules

If applicable:

I have added tests that prove my fix is effective or that my feature works
New and existing unit tests pass with my changes (locally and CI/CD build)
I have commented my code, particularly in hard-to-understand areas
I have made corresponding changes to the documentation
I have noted that this is a breaking change for a major release (fix or feature that would cause existing functionality to not work as expected)

codecov-commenter · 2022-09-10T17:48:07Z

Codecov Report

Base: 100.00% // Head: 100.00% // No change to project coverage 👍

Coverage data is based on head (f4a99b3) compared to base (147b80b).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff            @@
##              main      #340   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           14        14           
  Lines         1201      1200    -1     
=========================================
- Hits          1201      1200    -1

Impacted Files	Coverage Δ
xcdat/spatial.py	`100.00% <100.00%> (ø)`

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

pochedls · 2022-09-10T19:52:15Z

xcdat/spatial.py

            "Bounds",
            {"domain": xr.DataArray, "region": Optional[RegionAxisBounds]},
        )
-        axis_bounds: Dict[SpatialAxis, Bounds] = {
-            "X": {
+        axis_bounds: Dict[SpatialAxis, Bounds] = {}
+        if "X" in axis:
+            axis_bounds["X"] = {
                "domain": self._dataset.bounds.get_bounds("X").copy(),
                "region": lon_bounds,
-            },
-            "Y": {
+            }
+        if "Y" in axis:
+            axis_bounds["Y"] = {
                "domain": self._dataset.bounds.get_bounds("Y").copy(),
                "region": lat_bounds,
-            },
-        }
+            }


Previously, the axis_bounds dict tried to populate X and Y every time, but we do not need to populate the bounds information if we are not performing spatial averaging over the X or Y axis. This fixes #324. I think this was a bug and I didn't see a need to modify the unit tests in this case.

Logic looks good. I did some minor refactoring and will push the code for your review.

pochedls · 2022-09-10T19:53:29Z

tests/test_spatial.py

+    def test_bounds_reordered_when_upper_indexed_first(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]]
+            ),
+            coords={"lat": self.ds.lat},
+            dims=["lat", "bnds"],
+        )
+        result = self.ds.spatial._force_domain_order_low_to_high(domain_bounds)

-        with pytest.raises(ValueError):
-            ds.spatial.get_weights(axis=["X"])
+        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)


Previously, we threw an error if the lower domain bound was larger than the upper domain bound. Now we re-order them. This unit test checks to make sure that is done in _force_domain_order_low_to_high.

pochedls · 2022-09-10T19:55:12Z

xcdat/spatial.py

-    def _validate_domain_bounds(self, domain_bounds: xr.DataArray):
-        """Validates the ``domain_bounds`` arg based on a set of criteria.
-
+    def _force_domain_order_low_to_high(self, domain_bounds: xr.DataArray):
+        """Re-orders the ``domain_bounds`` to be ordered low-to-high such that
+        domain_bounds[:, 0] < domain_bounds[:, 1].
        Parameters
        ----------
        domain_bounds: xr.DataArray
-            The bounds of an axis
-
-        Raises
+            The bounds of an axis.
+        Returns
        ------
-        TypeError
-            If the ``domain_bounds`` of a grid cell are not ordered low-to-high.
+        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:
-            raise ValueError(
-                "The bounds have unexpected ordering. A lower bound has a "
-                "greater value than an upper bound."
-            )
+            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
+        else:
+            return domain_bounds


Previously, we threw an error if the lower domain bound was larger than the upper domain bound (in _validate_domain_bounds) This function is replaced with _force_domain_order_low_to_high which re-orders the bounds such that domain_bound[:, 0] < domain_bounds[:, 1].

Note I had originally replaced these lines with logic to simply determine, which of the bounds is higher versus lower (e.g., domain_lowers = np.min(d_bounds, axis=1)), but the order matters further down (here).

pochedls · 2022-09-10T19:58:32Z

xcdat/spatial.py

-        if len(bounds) <= 0 or len(bounds) > 2:
+        if len(bounds) != 2:


I randomly caught this... I think <=0 doesn't really make any sense, right?

I agree with this change, and <=0 does not make sense (arrays can't have < 0 elements)

pochedls · 2022-09-11T00:55:08Z

I think this is in good shape and can be reviewed. I do plan on re-running it over a bunch of cases before merging.

pochedls · 2022-09-13T18:51:44Z

I re-did comparisons with CDAT - we still agree with CDAT in all cases where both can load the dataset (loop over 115 model files and 21 domains).

tomvothecoder

Looks good overall, with some minor refactoring changes in my upcoming commit.

tomvothecoder · 2022-09-14T16:15:24Z

xcdat/spatial.py

-        if len(bounds) <= 0 or len(bounds) > 2:
+        if len(bounds) != 2:


I agree with this change, and <=0 does not make sense (arrays can't have < 0 elements)

tomvothecoder · 2022-09-14T16:34:47Z

tests/test_spatial.py

+    def test_bounds_reordered_when_upper_indexed_first(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]]
+            ),
+            coords={"lat": self.ds.lat},
+            dims=["lat", "bnds"],
+        )
+        result = self.ds.spatial._force_domain_order_low_to_high(domain_bounds)

-        with pytest.raises(ValueError):
-            ds.spatial.get_weights(axis=["X"])
+        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)


tomvothecoder · 2022-09-14T16:35:14Z

xcdat/spatial.py

            "Bounds",
            {"domain": xr.DataArray, "region": Optional[RegionAxisBounds]},
        )
-        axis_bounds: Dict[SpatialAxis, Bounds] = {
-            "X": {
+        axis_bounds: Dict[SpatialAxis, Bounds] = {}
+        if "X" in axis:
+            axis_bounds["X"] = {
                "domain": self._dataset.bounds.get_bounds("X").copy(),
                "region": lon_bounds,
-            },
-            "Y": {
+            }
+        if "Y" in axis:
+            axis_bounds["Y"] = {
                "domain": self._dataset.bounds.get_bounds("Y").copy(),
                "region": lat_bounds,
-            },
-        }
+            }


Logic looks good. I did some minor refactoring and will push the code for your review.

tomvothecoder · 2022-09-14T16:35:22Z

xcdat/spatial.py

-    def _validate_domain_bounds(self, domain_bounds: xr.DataArray):
-        """Validates the ``domain_bounds`` arg based on a set of criteria.
-
+    def _force_domain_order_low_to_high(self, domain_bounds: xr.DataArray):
+        """Re-orders the ``domain_bounds`` to be ordered low-to-high such that
+        domain_bounds[:, 0] < domain_bounds[:, 1].
        Parameters
        ----------
        domain_bounds: xr.DataArray
-            The bounds of an axis
-
-        Raises
+            The bounds of an axis.
+        Returns
        ------
-        TypeError
-            If the ``domain_bounds`` of a grid cell are not ordered low-to-high.
+        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:
-            raise ValueError(
-                "The bounds have unexpected ordering. A lower bound has a "
-                "greater value than an upper bound."
-            )
+            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
+        else:
+            return domain_bounds


tomvothecoder

Hi @pochedls, please review my refactoring changes. The unit tests still pass.

If you agree with them, merge when ready.

tomvothecoder · 2022-09-14T16:38:43Z

xcdat/spatial.py


+            weights = axis_bounds[key]["weights_method"](d_bounds, r_bounds)
            weights.attrs = d_bounds.attrs
            axis_weights[key] = weights


I refactored a bit here to reduce if else statements and avoid fetching domain_bounds for unspecified keys in axis.

This commented was intended to highlight the entire method, not just this line.

tomvothecoder · 2022-09-14T16:39:41Z

xcdat/spatial.py

    def _force_domain_order_low_to_high(self, domain_bounds: xr.DataArray):
-        """Re-orders the ``domain_bounds`` to be ordered low-to-high such that
-        domain_bounds[:, 0] < domain_bounds[:, 1].
+        """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
        else:
            return domain_bounds


Updated docstring and code spacing for consistency with code style

pochedls · 2022-09-14T16:58:10Z

xcdat/spatial.py

+                "weights_method": self._get_latitude_weights,
+                "region": np.array(lat_bounds, dtype="float")
+                if lat_bounds is not None
+                else None,


You learn something new every day (that if can come after the statement).

Yup, they are called ternary operators (aka inline conditionals).
https://datagy.io/inline-if-in-python/

I don't use them too often for variable assignment, but it can be useful where it makes sense (readable and simple).

xcdat/spatial.py

Co-authored-by: Tom Vo <[email protected]>

initial fixes for #324 and #338

55f38a2

pochedls commented Sep 10, 2022

View reviewed changes

comment to explain bound ordering

7b23208

pochedls requested a review from tomvothecoder September 11, 2022 00:53

tomvothecoder added type: bug Inconsistencies or issues which will cause an issue or problem for users or implementors. Priority: Medium labels Sep 12, 2022

tomvothecoder assigned pochedls Sep 12, 2022

tomvothecoder reviewed Sep 14, 2022

View reviewed changes

Add PR review suggestions

6b56cee

tomvothecoder approved these changes Sep 14, 2022

View reviewed changes

pochedls commented Sep 14, 2022

View reviewed changes

tomvothecoder reviewed Sep 14, 2022

View reviewed changes

xcdat/spatial.py Outdated Show resolved Hide resolved

pochedls and others added 2 commits September 14, 2022 12:35

simplify return statement in spatial.py

34dd7f8

Co-authored-by: Tom Vo <[email protected]>

Merge branch 'main' into bug/fix-spatial-domain-ordering

f4a99b3

pochedls merged commit 2ea2f46 into main Sep 15, 2022

pochedls deleted the bug/fix-spatial-domain-ordering branch September 15, 2022 04:34

tomvothecoder mentioned this pull request Sep 16, 2022

[Feature]: Handle "out of order" bounds #338

Closed

4 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix spatial bugs: handle datasets with domain bounds out of order and zonal averaging #340

Fix spatial bugs: handle datasets with domain bounds out of order and zonal averaging #340

pochedls commented Sep 10, 2022 •

edited by tomvothecoder

Loading

codecov-commenter commented Sep 10, 2022 •

edited

Loading

pochedls Sep 10, 2022

tomvothecoder Sep 14, 2022

pochedls Sep 10, 2022

tomvothecoder Sep 14, 2022

pochedls Sep 10, 2022

tomvothecoder Sep 14, 2022

pochedls Sep 10, 2022

tomvothecoder Sep 14, 2022

pochedls commented Sep 11, 2022

pochedls commented Sep 13, 2022

tomvothecoder left a comment

tomvothecoder Sep 14, 2022

tomvothecoder Sep 14, 2022

tomvothecoder Sep 14, 2022

tomvothecoder Sep 14, 2022

tomvothecoder left a comment

tomvothecoder Sep 14, 2022 •

edited

Loading

tomvothecoder Sep 14, 2022

tomvothecoder Sep 14, 2022

pochedls Sep 14, 2022 •

edited

Loading

tomvothecoder Sep 14, 2022 •

edited

Loading

Fix spatial bugs: handle datasets with domain bounds out of order and zonal averaging #340

Fix spatial bugs: handle datasets with domain bounds out of order and zonal averaging #340

Conversation

pochedls commented Sep 10, 2022 • edited by tomvothecoder Loading

Description

Checklist

codecov-commenter commented Sep 10, 2022 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pochedls commented Sep 11, 2022

pochedls commented Sep 13, 2022

tomvothecoder left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tomvothecoder left a comment

Choose a reason for hiding this comment

tomvothecoder Sep 14, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pochedls Sep 14, 2022 • edited Loading

Choose a reason for hiding this comment

tomvothecoder Sep 14, 2022 • edited Loading

Choose a reason for hiding this comment

pochedls commented Sep 10, 2022 •

edited by tomvothecoder

Loading

codecov-commenter commented Sep 10, 2022 •

edited

Loading

tomvothecoder Sep 14, 2022 •

edited

Loading

pochedls Sep 14, 2022 •

edited

Loading

tomvothecoder Sep 14, 2022 •

edited

Loading