cos-lat averaging #125

RondeauG · 2022-12-06T14:44:14Z

Pull Request Checklist:

pre-commit hooks are installed/active in my local clone ($ pre-commit install)
This PR addresses an already opened issue (for bug fixes / features)
- This PR fixes #xyz
(If applicable) Documentation has been added / updated (for bug fixes / features)
HISTORY.rst has been updated (with summary of main changes)
- Link to issue (:issue:number) and pull request (:pull:number) has been added

What kind of change does this PR introduce?

Added cos-lat method to aggregate.spatial_mean()
Deprecated mean and added warnings to xesmf
Reworked the warning messages in _subset_file_coverage to be more explicit
Fixed a bug in compute_indicators when less than 3 timesteps were produced.

Does this PR introduce a breaking change?

No, but deprecates mean and changes the name of interp_coord to interp_centroid (still supports interp_coord with a deprecation warning)

Other information:

RondeauG · 2022-12-06T14:45:34Z

This is a partial fix to #94. Fixing xesmf.SpatialAverager is harder than expected and might require a PR in xESMF itself, but I didn't want to wait that long to implement cos-lat averaging.

aulemahal

A few comments. I don't think the "cos-lat" think is compatible with curvilinear grids, but I'm not sure why it shouldn't be.

xscen/aggregate.py

aulemahal · 2022-12-09T20:39:09Z

xscen/indicators.py

+                if len(out.time) < 3:
+                    freq = (
+                        ind.injected_parameters["freq"]
+                        if "freq" in ind.injected_parameters
+                        else ind.parameters["freq"]["default"]
+                        if "freq" in ind.parameters
+                        else ind.src_freq
+                    )
+                else:
+                    freq = xr.infer_freq(out.time)
+            else:
+                freq = "fx"


Do you think ii would be better to wrap this in a private function instead of having it twice?

mccrayc

Two key points for the cos-lat averaging:

No need to normalize the weights by their sum, I'm pretty sure xarray does this already.
In agreement with what Pascal said, I believe cos-lat weighting should work fine with 2D curvilinear grids with a couple minor modifications (see comments).

xscen/aggregate.py

aulemahal · 2022-12-12T15:52:21Z

Je me reprends, cos-lat as implemented is indeed not usable with curvilinear grids. The thing is, cos(lat) is a weight that varies in latitude and thus assumes all cells have the same area in °² (degrees squared). Which is not true for rotated pole for example.

I found this blog : https://towardsdatascience.com/the-correct-way-to-average-the-globe-92ceecd172b7
It shows a method approximating the cell area with a lat-dependent earth radius. It however uses a cell area approximation that only works for non-curvilinear grids. You could use cf_xarray 0.7.6 to compute 2D bounds and then cell areas (in °).

@RondeauG your code was indeed correct. But I think curvilinear grid support is essential as one of the principal use case is with MRCC5 data.

aulemahal · 2022-12-12T18:49:49Z

@RondeauG I found the challenge interesting and spent a bit too much time on this, but here are 3 ways of computing cell areas for the weighted average:
https://gist.github.com/aulemahal/13c907cf2f7ccd3a84be34a7a009516c

We might not have pygeos already, but it's an optional dependency of clisops, so it might be worth adding it explicitly, it is much faster!

However, there's another issue with cf-xarray... Computing the 2D bounds yields strange results when the cell crosses the first lon meridian. In the MRCC5 case, the domain corner crosses the +180 / -180 line and cell areas are super large because the corners are wrong.

See:

# Conflicts: # HISTORY.rst # xscen/aggregate.py

RondeauG · 2022-12-20T22:19:36Z

As a reference to compare to, I implemented the cos-lat method for 2D grids using ds.cf['latitude'] instead of ds.cf['X']. Here's an example for RegCM4 / NAM-22:

RondeauG · 2022-12-20T22:20:22Z

I'll look into the links provided, though I feel like this might be another method entirely? As in, we'd have cos-lat or SpatialAverager or pygeos/ESMF/whatever?

RondeauG · 2022-12-20T22:22:36Z

Also, with the changes, get_spatial_dims(ds) is pretty useless. I could restore it within method=="mean" as it was before, but I'm tempted to just get rid of mean entirely since it does the wrong thing to begin with.

aulemahal · 2022-12-22T15:38:39Z

this might be another method entirely?

My "shapely" and "pygeos" functions are both based on the cos-lat approximation, so I would argue it's just a curvilinear-compatible version of cos-lat.

The ESMF version is indeed another technique because it computes the cell area using other approximations.

Note for shapely: the newest version (2.0) now merges the functionality of pygeos so there might be a way to code this idea without adding a new dependency.

aulemahal · 2022-12-22T16:08:12Z

I updated my gist with a function that works with shapely 2.0. It is indeed the exact same as with pygeos and it's as fast.

# Conflicts: # HISTORY.rst

RondeauG · 2023-03-07T16:57:25Z

@aulemahal @mccrayc
This might be a lack of comprehension from my part on rotated grids, but let's say that I have this region:

The first part of Pascal's code (using the pygeos version for example) produces this:
xr.DataArray(pygeos.area(pygeos.polygons(pygeos.linearrings(ds.lon_bounds, ds.lat_bounds))), dims=ds.lon.dims, coords=ds.lon.coords).plot()

np.cos(np.deg2rad(ds.lat)) creates more or less the opposite:

Meaning that when you multiply both, you get:

...which makes no intuitive sense to me? I would have expected weights to be curved, not linear like that?

aulemahal · 2023-03-07T21:00:17Z

No I think this is ok! You are in a Oblique Mercator projection, which projects the globe unto a cylinder. The projection being centered at y = 0, in the center of your plot, it makes sense to me that the cell area decreases with abs(y).

https://en.wikipedia.org/wiki/Oblique_Mercator_projection

Also, I found a way to compute areas (in meters) with pyproj from lat and lon bounds. This was added to the gist, and gives the following:

# Conflicts: # HISTORY.rst # environment-dev.yml # environment.yml

aulemahal

This looks good!

Semi-related:
As we currently have no bulletproof method, the choice makes sense. However, it seems to me that there's no reason to have such an array of methods for such a simple task. I guess that when we have a simple, light area-weighted mean method, we could remove everything else? Is there a use for a collection of imperfect algos ?

This would mean :

fix bounds generation when crossing meridians in cf_xarray
generate areacella in m² from those bounds (using pyproj and shapely)

xscen/aggregate.py

Co-authored-by: Pascal Bourgault <bourgault.pascal@ouranos.ca>

for more information, see https://pre-commit.ci

RondeauG · 2023-03-09T14:21:04Z

Semi-related: As we currently have no bulletproof method, the choice makes sense. However, it seems to me that there's no reason to have such an array of methods for such a simple task. I guess that when we have a simple, light area-weighted mean method, we could remove everything else? Is there a use for a collection of imperfect algos ?

That's why I deprecate mean, since cos-lat is a superior version of the same thing. I could see an argument to drop interp_centroid in favor of always using xesmf too.

However, I think that cos-lat and xesmf (SpatialAverager) fundamentaly do different things. cos-lat keeps the grid cells as is, while SpatialAverager will break pixels into smaller parts to closely fit to a given polygon.

aulemahal · 2023-03-09T18:32:01Z

@mccrayc When you have time, a second review is needed since you "requested changes".

mccrayc

Looks good to me!

xscen/aggregate.py

xscen/indicators.py

xscen/aggregate.py

Co-authored-by: Pascal Bourgault <bourgault.pascal@ouranos.ca>

added cos-lat averaging, other bugfixes

c967e91

RondeauG requested review from aulemahal and mccrayc December 6, 2022 14:46

upd History

ecb08e3

aulemahal requested changes Dec 9, 2022

View reviewed changes

mccrayc requested changes Dec 11, 2022

View reviewed changes

xscen/aggregate.py Outdated Show resolved Hide resolved

xscen/aggregate.py Outdated Show resolved Hide resolved

RondeauG added 2 commits December 20, 2022 15:42

Merge branch 'main' into fix-94

dd316d5

# Conflicts: # HISTORY.rst # xscen/aggregate.py

support for 2D cos-lat

08404ab

RondeauG and others added 4 commits January 24, 2023 10:46

Merge branch 'main' into fix-94

21825d7

# Conflicts: # HISTORY.rst

Merge branch 'main' into fix-94

decd909

Merge branch 'main' into fix-94

2e2e952

# Conflicts: # HISTORY.rst

Merge remote-tracking branch 'origin/fix-94' into fix-94

5cbf42b

# Conflicts: # HISTORY.rst

updated requirements, suggested changes from code review

4f5de5e

RondeauG added 3 commits March 7, 2023 16:24

upd notebook

519113a

Merge branch 'main' into fix-94

a8caee9

# Conflicts: # HISTORY.rst # environment-dev.yml # environment.yml

fix xesmf not setting 'axis'

37c06e5

RondeauG requested review from aulemahal and mccrayc March 7, 2023 22:38

aulemahal approved these changes Mar 8, 2023

View reviewed changes

xscen/aggregate.py Outdated Show resolved Hide resolved

xscen/aggregate.py Outdated Show resolved Hide resolved

Update xscen/aggregate.py

268c2c7

Co-authored-by: Pascal Bourgault <bourgault.pascal@ouranos.ca>

[pre-commit.ci] auto fixes from pre-commit.com hooks

dc2da01

for more information, see https://pre-commit.ci

RondeauG added 2 commits March 9, 2023 09:44

fix bug in flox attributes, change wording in history

966910d

Merge remote-tracking branch 'origin/fix-94' into fix-94

096051e

mccrayc approved these changes Mar 9, 2023

View reviewed changes

RondeauG commented Mar 10, 2023

View reviewed changes

xscen/aggregate.py Outdated Show resolved Hide resolved

Update xscen/aggregate.py

1789aa7

RondeauG commented Mar 10, 2023

View reviewed changes

xscen/indicators.py Outdated Show resolved Hide resolved

Update xscen/indicators.py

6d49582

aulemahal reviewed Mar 10, 2023

View reviewed changes

xscen/aggregate.py Outdated Show resolved Hide resolved

Update xscen/aggregate.py

21f1243

Co-authored-by: Pascal Bourgault <bourgault.pascal@ouranos.ca>

RondeauG merged commit eaff2df into main Mar 10, 2023

RondeauG deleted the fix-94 branch March 10, 2023 15:27

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

cos-lat averaging #125

cos-lat averaging #125

RondeauG commented Dec 6, 2022 •

edited

Loading

RondeauG commented Dec 6, 2022

aulemahal left a comment

aulemahal Dec 9, 2022

mccrayc left a comment •

edited

Loading

aulemahal commented Dec 12, 2022

aulemahal commented Dec 12, 2022

RondeauG commented Dec 20, 2022

RondeauG commented Dec 20, 2022

RondeauG commented Dec 20, 2022

aulemahal commented Dec 22, 2022 •

edited

Loading

aulemahal commented Dec 22, 2022

RondeauG commented Mar 7, 2023

aulemahal commented Mar 7, 2023

aulemahal left a comment

RondeauG commented Mar 9, 2023

aulemahal commented Mar 9, 2023

mccrayc left a comment

cos-lat averaging #125

cos-lat averaging #125

Conversation

RondeauG commented Dec 6, 2022 • edited Loading

Pull Request Checklist:

What kind of change does this PR introduce?

Does this PR introduce a breaking change?

Other information:

RondeauG commented Dec 6, 2022

aulemahal left a comment

Choose a reason for hiding this comment

aulemahal Dec 9, 2022

Choose a reason for hiding this comment

mccrayc left a comment • edited Loading

Choose a reason for hiding this comment

aulemahal commented Dec 12, 2022

aulemahal commented Dec 12, 2022

RondeauG commented Dec 20, 2022

RondeauG commented Dec 20, 2022

RondeauG commented Dec 20, 2022

aulemahal commented Dec 22, 2022 • edited Loading

aulemahal commented Dec 22, 2022

RondeauG commented Mar 7, 2023

aulemahal commented Mar 7, 2023

aulemahal left a comment

Choose a reason for hiding this comment

RondeauG commented Mar 9, 2023

aulemahal commented Mar 9, 2023

mccrayc left a comment

Choose a reason for hiding this comment

RondeauG commented Dec 6, 2022 •

edited

Loading

mccrayc left a comment •

edited

Loading

aulemahal commented Dec 22, 2022 •

edited

Loading