-
Notifications
You must be signed in to change notification settings - Fork 2
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
Cice grid #5
Cice grid #5
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #5 +/- ##
===========================================
- Coverage 100.00% 96.98% -3.02%
===========================================
Files 3 5 +2
Lines 190 232 +42
===========================================
+ Hits 190 225 +35
- Misses 0 7 +7 ☔ View full report in Codecov by Sentry. |
@aekiss - Sorry - third time lucky for this PR! Can you review please? |
I did only very short comparisons using a 1 degree config in OM2. The sea ice area and volume look practically the same: Although not identical (anomaly in aice of new grid - old grid after 1 year): I ran the 0.25 om2 for 10 years with the new grid. I didn't do a comparison run, but did compare to some save runs. Again nothing particularly of note in sea ice area and volume. There are more differences in concentration, although I think these will be due to the change in land mask between runs: Link for reference: https://github.com/anton-seaice/sandbox/blob/c745853d63766ca0bb0d7ac477c52d78c3ad69cf/new_cice_grid/new-grid-intake-cat.ipynb |
Thanks @anton-seaice this looks awesome. Thanks especially for the detailed documentation and provenance. For clarity, should this bit in your initial post
be changed to
to match the link targets? |
tests/test_cice_grid.py
Outdated
test_grid["ulat"] = deg2rad(ds.y.isel(nxp=slice(2, None, 2), nyp=slice(2, None, 2))) | ||
test_grid["ulon"] = deg2rad(ds.x.isel(nxp=slice(2, None, 2), nyp=slice(2, None, 2))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Am I right in guessing that isel slices start from 2 rather than 0 because 0 corresponds to the NE corner of cells that are not part of the actual grid?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If so, perhaps add a comment to that effect as it is a bit confusing otherwise
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Am I right in guessing that isel slices start from 2 rather than 0 because 0 corresponds to the NE corner of cells that are not part of the actual grid?
Yes that's right, in the old OM2 grid ulat
starts from zero led to mismatch in longitudes
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did a bit of a refactor to hopefully clarify. Its hard to write this stuff concisely in comments.
tests/test_cice_grid.py
Outdated
area_wrapped = mom_grid.ds.area | ||
area_wrapped = xr.concat([ds.area.isel(nx=slice(1, None)), ds.area.isel(nx=0)], dim="nx") | ||
|
||
top_row = xr.concat([ds.area.isel(ny=-1, nx=slice(-2, 0, -1)), ds.area.isel(ny=-1, nx=[-1, 0])], dim="nx") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
top_row = xr.concat([ds.area.isel(ny=-1, nx=slice(-2, 0, -1)), ds.area.isel(ny=-1, nx=[-1, 0])], dim="nx") | |
top_row = area_wrapped.isel(ny=-1, nx=slice(-1, 0, -1)) |
Do we just want a reversed copy of the wrapped top row here, or have I misunderstood?
(Or perhaps the folded top row should be concatenated before wrapping? The order of these operations affects the index of the tripole singularities - see fig 9.9 of Griffies - note that the singularities fall on u points 0 and ni/2)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I tried this suggestion and the resulting areas were offset by 1 cell ... I will have to investigate. I am confident the areas from this method are right, but maybe the method is wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah ...
xarray is confusing sometimes.
ds.area.isel(ny=-1, nx=[-1, 0])
and ds.area.isel(ny=-1, nx=[0,-1])
are equal :)
Anyway the slightly better version of this line would be:
xr.concat([ds.area.isel(ny=-1, nx=slice(-2, None, -1)), ds.area.isel(ny=-1, nx=-1)], dim="nx")
But the even better version, is to fold first, then wrap, and then its generally clearer in total!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Changed to fold first, then wrap
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for doing this @anton-seaice! An unsolicited review from me with a comment about the organisation of the code across repos, rather than about the specifics of the code itself.
I think it's a great idea to use and extend the esmgrids
package, since it includes a lot of the functionality we need already. Currently however, there's a bit of scope creep across different repos/packages:
esmgrids
generates the CICE grid from the MOM grid and adds a little bit of provenance metadata (git info)om3-utils/cice_grid.py
usesesmgrids
and adds a bit more metadata (input files, hashes, run command, crs)om3-utils/test_cice_grid.py
tests mostly thatesmgrids
is doing what it's supposed to.
I propose moving all of this functionality into the esmgrids
package. It certainly makes sense for the tests to live there, and I don't see why all the metadata couldn't be added there.
As for cmdline scripts, we could either provide these as console-scripts
in the esmtools
package, or we could add them to om3-scripts
.
Of course, we’ll need to set up proper test and release workflows for esmgrid
so that we can use it effectively and reliably across other projects, but this needs to be done regardless.
That makes sense. |
Yeah, I think there's a little bit of work needed across the board (both infrastructure and code). That would be great if you can look at installation and testing. If we're going to use |
Thanks @dougiesquire I mostly agree, I want to be sure that we aren't breaking anything historical (e.g. for OM2) and note that we would need to tidyup the old testing setup etc in esmgrids. There also is some overlap / redundancy with mesh generation (om3-scripts). Hopefully we can be a little bit careful with metadata, especially the CRS, to make sure we can support regional and non-mercator grids in the future. Can we move this to a new issue, so we can use this PR to make a cice grid for OM2 release and to get us set-up to do testing with esmgrids ? |
Agreed. I thought I opened an issue about this, but I can't find it now... Anyway, I've just opened this: COSIMA/om3-scripts#11
So you would like to add this to |
Not ideal, but also not a big deal.
Sorry, I got distracted with the profiling stuff. I'll get to it ASAP. Should indeed be straightforward. |
FYI: There are some historical tests in the esmgrids repo, but they rely on some test data which I don't think we have anymore. |
Assuming the tests were passing last time there were any changes in the repo, one could simply take some random grid and regenerate the test data, no? |
Looks like it for the MOM to cice grid: I've replaced this with using the ocean_grid_generator in this PR I don't know what this input is though: (p.s. we can definitely just leave it failing for now.) |
Sorry, I wasn't meaning to suggest things should be happening faster. I was trying (unclearly) to ask whether getting testing running on |
Mostly i'm worried it might send us down a rabbit hole for a while :) |
FYI @anton-seaice, Aidan joined our team meeting yesterday. Their plan is to release the 1deg and 0.1deg configs this week and then follow up with the 025deg configs next week once this is all bedded down. So there's a little bit of time. |
Co-authored-by: Andrew Kiss <[email protected]>
@aekiss ... I responded to the comments. Looks like well move this to esmgrids but let me know if the new changes make sense |
This is being implemented via COSIMA/esmgrids#6 |
This change updates the workflow to generate a cice grid from a mom-supergrid file containing a tripolar grid. - fixes angle T - adds cf compliant variable names - add versioning to output - adds an end-end test COSIMA/om3-utils#5 describes full scope of changes --------- Co-authored-by: dougiesquire <[email protected]> Co-authored-by: Dougie Squire <[email protected]>
Following on from COSIMA/om3-scripts#8 but now with tests!
Created a cice grid generation script, using the esmgrids package.
The script creates a cice grid from the mom supergrid and the ocean mask file. The scripts adds the git commit of this script and the input path and md5sum of the mom files to the netcdf output. I added as close as we can get to cf-compliance and a CRS (See COSIMA/om3-scripts#7).
The differences between versions are in ... notebook and explored below in more detail. By the look of things, 1 deg and 0.25 degree files are inherited / historical and 0.1 deg was made for OM2 using the esmgrids package.
Differences ( I filtered differences less than 2e-6):
1deg
new vars not in old?
{'crs'}
missing vars in new?
{'latt_bonds', 'lonu_bonds', 'hue', 'hun', 'latu_bonds', 'lont_bonds'}
tlon anom min: nan, anom max: nan
angleT anom min: nan, anom max: nan
tlat anom min: nan, anom max: nan
uarea anom min: -1186349698.6215067, anom max: 33617654.79803729
hte anom min: nan, anom max: nan
ulon anom min: nan, anom max: nan
tarea anom min: -19085417.985637665, anom max: 9877771.11514306
angle anom min: nan, anom max: nan
htn anom min: nan, anom max: nan
ulat anom min: nan, anom max: nan
025deg
new vars not in old?
{'crs'}
missing vars in new?
set()
tlon anom min: nan, anom max: nan
angleT anom min: nan, anom max: nan
tlat anom min: nan, anom max: nan
uarea anom min: -69552757.44506374, anom max: 7366673.800787419
hte anom min: nan, anom max: nan
ulon anom min: 3.620558794636963e-05, anom max: 3.141592653589793
tarea anom min: nan, anom max: nan
angle anom min: -1.5707963267948966, anom max: 1.5707963267948966
htn anom min: nan, anom max: nan
ulat anom min: -0.0018422347075643941, anom max: 0.0018422347075643941
01deg
new vars not in old?
{'crs'}
missing vars in new?
{'clat_u', 'clon_t', 'clat_t', 'clon_u'}
tlon anom min: nan, anom max: nan
angleT anom min: -3.139931603248902, anom max: 1.1133153100325133
tlat anom min: nan, anom max: nan
uarea anom min: nan, anom max: nan
hte anom min: nan, anom max: nan
ulon anom min: nan, anom max: nan
tarea anom min: nan, anom max: nan
angle anom min: nan, anom max: nan
htn anom min: nan, anom max: nan
ulat anom min: nan, anom max: nan
Even ignoring the join, its not clear why 0.25 deg uarea has changed, but it looks like it was set equal to tarea in the old grid
Known issues:
See companion PR for esmgrids repo
To-do: