Adds support for parametric vertical coordinate #528
Conversation
cf_xarray/parametric.py
Outdated
| return getattr(m, stdname) | ||
|
|
||
|
|
||
| def derive_dimension_order(output_order, **dim_map): |
There was a problem hiding this comment.
This is quite clever, but I'd prefer if the user just did the transpose on their own after this decoding.
There was a problem hiding this comment.
Is it fine that this would break away from the convention since it provides the output ordering in the definition?
There was a problem hiding this comment.
Yes with Xarray we care much less about dimension order.
There was a problem hiding this comment.
@jasonb5 / @dcherian – thanks for all of your work on this. I wanted to weigh-in as a user. It would be unexpected to get a matrix of [time, lat, lon, *plev*] instead of [time, *plev*, lat, lon]. I think users would prefer the former (matches most climate data, CF conventions, and almost certainly the order of the original data, which is typically [time, *lev*, lat lon]).
If this doesn't match cf-xarray's style, xcdat could apply this functionality after calling cf-xarray's transform.
There was a problem hiding this comment.
Yes xcdat could do .cf.transpose('T', 'Z', 'Y', 'X')
almost certainly the order of the original data, which is typically [time, lev, lat lon]).
Is this ordering not preserved here?
There was a problem hiding this comment.
I just did a quick check: it looks like it goes from (ds.ta.dims):
('time', 'lev', 'lat', 'lon')
to (p = ds.a * ds.p0 + ds.b * ds.ps; p.dims):
('lev', 'time', 'lat', 'lon')
In other words, the dataarray of interest is [time, lev, lat, lon], but the transform will put lev at the front [lev, time, lat, lon]
|
Thanks for taking this on! I've left a few high-level comments. Mostly I think it would be a lot clearer if we define a dataclass per transform that took care of the all the specialized logic: #528 (comment) |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #528 +/- ##
==========================================
+ Coverage 85.78% 93.47% +7.68%
==========================================
Files 13 13
Lines 2364 2190 -174
Branches 183 0 -183
==========================================
+ Hits 2028 2047 +19
+ Misses 303 143 -160
+ Partials 33 0 -33
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. |
cf_xarray/parametric.py
Outdated
| z_shape = list(self.eta.shape) | ||
|
|
||
| z_shape.insert(1, self.sigma.shape[0]) | ||
|
|
||
| z_dims = list(self.eta.dims) | ||
|
|
||
| z_dims.insert(1, self.sigma.dims[0]) |
There was a problem hiding this comment.
| z_shape = list(self.eta.shape) | |
| z_shape.insert(1, self.sigma.shape[0]) | |
| z_dims = list(self.eta.dims) | |
| z_dims.insert(1, self.sigma.dims[0]) | |
| z_shape = (self.sigma.shape[0], *self.eta.shape) | |
| z_dims = (self.sigma.dims[0], *self.eta.dims) |
cf_xarray/parametric.py
Outdated
| z = xr.DataArray(np.empty(z_shape), dims=z_dims) | ||
|
|
||
| z_sigma = self.eta + self.sigma * ( | ||
| np.minimum(self.depth_c, self.depth) + self.eta | ||
| ) | ||
|
|
||
| z = xr.where(~np.isnan(self.sigma), z_sigma, z) | ||
|
|
||
| z = xr.where(np.isnan(self.sigma), self.zlev, z) |
There was a problem hiding this comment.
| z = xr.DataArray(np.empty(z_shape), dims=z_dims) | |
| z_sigma = self.eta + self.sigma * ( | |
| np.minimum(self.depth_c, self.depth) + self.eta | |
| ) | |
| z = xr.where(~np.isnan(self.sigma), z_sigma, z) | |
| z = xr.where(np.isnan(self.sigma), self.zlev, z) | |
| z_sigma = self.eta + self.sigma * ( | |
| np.minimum(self.depth_c, self.depth) + self.eta | |
| ) | |
| z = xr.where(np.isnan(self.sigma), self.zlev, z_sigma) |
Does this work?
|
ignore the 3.9 failures, I'll drop it soon. |
|
Thanks! This is getting close. Just some style nits. I haven't looked at the tests but I'm assuming you have it covered with some real life usage too? :) |
|
EDIT: I added a checklist to the top. |
| f + (self.sigma - 1) * (self.depth - f), | ||
| ) | ||
|
|
||
| return z.squeeze().assign_attrs(standard_name=self.computed_standard_name) |
There was a problem hiding this comment.
| return z.squeeze().assign_attrs(standard_name=self.computed_standard_name) | |
| return z.assign_attrs(standard_name=self.computed_standard_name) |
Can we delete all the squeeze too please
There was a problem hiding this comment.
@dcherian Just want to check, if I remove the squeeze we'll end up with placeholder dimensions e.g. (lev: 3, dim_0: 1, lat: 2, lon: 2). If that's fine then I'll go ahead and remove the call.
There was a problem hiding this comment.
where is the dim_0 being created?
Co-authored-by: Deepak Cherian <[email protected]>
for more information, see https://pre-commit.ci
|
|
||
|
|
||
| def test_ocean_double_sigma_coordinate(): | ||
| k_c = xr.DataArray( |
There was a problem hiding this comment.
shouldn't these be scalars? How are they stored in the netCDF files you're looking at?
Adds support for parametric vertical coordinates from CFConventions.
squeeze