Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue with grid.transform #600

Open
mkyoungs opened this issue Apr 25, 2023 · 8 comments
Open

Issue with grid.transform #600

mkyoungs opened this issue Apr 25, 2023 · 8 comments
Labels

Comments

@mkyoungs
Copy link

mkyoungs commented Apr 25, 2023

Hi, I am trying to calculate thickness weighted u velocity in MITgcm and transform it back to depth space using the time averaged density structure. See the code below.
zexpand = coords.Z.expand_dims(dim={"i": 1080,'j':310})
ztrans = grid.transform(zexpand, 'Z', coords.layer_1RHO_center.values, target_data=(ds.dens-1000))
ztrans = ztrans.rename(dens='l1_c')
ures = grid.interp(ds.LaUH1RHO.mean(dim='time')/ds.LaHw1RHO.mean(dim='time'),'X')
ures_z = grid.transform(ures,'l1_c',coords.Z,target_data=ztrans)

Everything seems to work right until the last step, when I get the error:
2607 def transform(self, da, axis, target, **kwargs):
2608 """Convert an array of data to new 1D-coordinates along 'axis'.
2609 The method takes a multidimensional array of data 'da' and
2610 transforms it onto another data_array 'target_data' in the
(...)
2688
2689 """
-> 2691 ax = self.axes[axis]
2692 return ax.transform(da, target, **kwargs)

KeyError: 'l1_c'

I am not sure what this means or how to fix it. Does anyone have any advice?

@rcaneill
Copy link
Contributor

rcaneill commented Apr 26, 2023

Hello Madeleine,

I think that I know what is happening: the xgcm grid contains yours axes (usually X, Y, Z, T), and these are the axes that you refer to when using the grid methods.

Here when you use grid.transform(ures,'l1_c',coords.Z,target_data=ztrans), you tell your grid to use the axis l1_c, which you have not defined in your grid.
So I guess that you will need to redefine your grid before doing the last transformation, e.g.

ds['ztrans'] = ztrans
ds['ures'] = ures
grid_with_density_axis = xgcm.Grid(ds, coords={'S':{'center':'l1_c'}}, metrics={'S':['nae_of_density_metric']})
# Use S for density axis, S as in sigma

And then use grid_with_density_axis.transform(ures,'S',coords.Z,target_data=ztrans)

I haven't tested my code snippet, so you'll probably need to adapt slightly

@mkyoungs
Copy link
Author

mkyoungs commented Apr 26, 2023

Thanks for the help! I think the issue may stem from the fact that although the 1RHO coordinate (l1_c) and the adjacent outer points are categorized as grid coordinates in ds, they are not when I create the grid object. How would I just add the 1RHO coordinate to the grid object manually to deal with the lack of automatic reading in? The same issue seems to be showing up in the documentation in the MITgcm example in the section 'Creating the grid object. ' It claims to create the 1RHO axis in the grid object, but the output in the lines above don't include it.

@rcaneill
Copy link
Contributor

Could you show the attributes you have for l1_c ?

Otherwise, you can specify by hand the axes, by giving the coords argument to xgcm.Grid: https://xgcm.readthedocs.io/en/latest/api.html#grid

@mkyoungs
Copy link
Author

I can specify it manually, but then I would have to specify every axis manually. I am confused about how to do that with all of the metrics involved. Shouldn't this be updating automatically? See the dataset output. It seems to recognize the density coordinate system.
Screenshot 2023-04-27 at 12 52 37 PM

@mkyoungs
Copy link
Author

I realized I may not have answered your question. Here are the attributes.
Screenshot 2023-04-27 at 12 59 49 PM

@rcaneill
Copy link
Contributor

So the attributes don't have the mandatory axis attribute. You can add it by hand before calling the grid:

ds.l1_c.attrs['axis'] = 'S'
ds.l1_i.attrs['axis'] = 'S'
ds.l1_b.attrs['axis'] = 'S'

You will also need to add the c_grid_axis_shift attribute to -0.5 for left points, and 0.5 for right points (or whatever for inner/outer point), which I guess looking at your 1st snapshot would be:

ds.l1_i.attrs['c_grid_axis_shift'] = 0.5
ds.l1_b.attrs['c_grid_axis_shift'] = 0.5

Here is the link to the doc:
https://xgcm.readthedocs.io/en/latest/grids.html#comodo-data

@rcaneill
Copy link
Contributor

Regarding metrics, I don't think that xgcm parses them automatically, you have to provide them, see Cell [4] in https://xgcm.readthedocs.io/en/latest/grid_metrics.html#Using-metrics-with-xgcm

metrics = {
    ('X',): ['dxC', 'dxG'], # X distances
    ('Y',): ['dyC', 'dyG'], # Y distances
    ('Z',): ['drW', 'drS', 'drC'], # Z distances
    ('X', 'Y'): ['rA', 'rAz', 'rAs', 'rAw'] # Areas
}
grid = Grid(ds, metrics=metrics)

And you will need to also add the metric for the density coordinate

@mkyoungs
Copy link
Author

Yes, that worked! Thank you.

@mkyoungs mkyoungs reopened this May 12, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants