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

Fixed online regridding code (file_regrid.py) for restart-to-restart file regridding only #266

Merged
merged 21 commits into from
Oct 11, 2023

Conversation

yantosca
Copy link
Contributor

@yantosca yantosca commented Oct 5, 2023

This PR fixes several issues with GCPy's online regridding program using xESMF in GCPy module gcpy/file_regrid.py. Several utility functions have been added to the gcpy/grid.py, gcpy/regrid.py, and gcpy/cstools.py modules.

NOTE: At this point, regridding between restart files has been tested. There are some issues with diagnostic regridding (involving the ordering of vertical levels not being consistent with the lev:positive attribute, see: geoschem/GCHP#21). We do not recommend regridding between GCHP diagnostic files.

Regridding 4x5 GCClassic to 2x2.5 GCClassic

NOTE: The grayscale came in when converting from PDF to PNG. Not sure why.
classic-to-classic sfc

Regridding C24 checkpoint to 4x5 GCClassic

checkpoint-to-classic sfc

Regridding 4x5 GCClassic to C24 checkpoint

classic-to-checkpoint sfc

Regridding C48 checkpoint to C24 checkpoint

checkpoint-to-checkpoint sfc

Regridding C48 checkpoint to C24 stretched grid

@lizziel: Not sure if the regridding is wrong or the plotting is wrong
checkpoint-to-stretched sfc

gcpy/cstools.py
- Restore code that raises an ImportError if pyproj and shapely
  are not found in the user's Python environment.  This had previously
  been disabled so as not to break cloud benchmarking.
- Add new function "get_cubed_sphere_res" to return the number of
  grid cells per cubed-sphere face
- Add new function "is_cubed_sphere_rst_grid" to determine if a
  cubed-sphere grid has checkpoint/restart file dimensions, i.e.
  (time, lev, lat, lon).
- Add new function "is_cubed_sphere_diag_grid" to determine if a
  cubed-sphere grid has MAPL History file dimensions, i.e.
  (time, lev, nf, Xdim, Ydim).
- Function "is_cubed_sphere_grid" now uses the is_cubed_sphere_diag_grid
  and is_cubed_sphere_rst_grid functions
- Implemented code cleanup suggestions from Pylint
- Updated comments and Pydoc headers
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/regrid.py
- Use the optimal import order as determined by Pylint
- Split long lines up into smaller lines where feasible
- Renamed "unravel_checkpoint_lat" to "unravel_checkpoint", and
  add handling for the "lon" dimension.
- Renamed "ravel_checkpoint_lat" to "ravel_checkpoint", and
  add handling for the "lon" dimension.
- Removed "return" after "else", and re-indented code
- Removed unused variables
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/examples/diagnostics/compare_diags.py
- Rename and flip GCHP restart variables before plotting
- Remove "print_to_screen" from the logic that determines if a
  table of differences will be printed.

gcpy/examples/plotting/plot_comparisons.py
- Rename and flip GCHP restart variables before plotting
- Restrict the comparisons to one variable, for simplicity
- Now use the "tkagg" X11 backend so that plots show on the screen
- Renamed "main" to "plot_comparisons".
- Use argparse to parse command-line arguments in a new "main" routine
- Accept ref, dev, varname, level arguments from the command line

gcpy/examples/plotting/plot_single_panel.py
- Rename and flip GCHP restart variables before plotting
- Restrict the comparisons to one variable, for simplicity
- Now use the "tkagg" X11 backend so that plots show on the screen
- Renamed "main" to "plot_single_level".
- Use argparse to parse command-line arguments in a new "main" routine
- Accept infile, varname, level arguments from the command line

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/regrid_restart_file.py
- Add a "to_gchp" argument so that we can change 'lev:positive="up"'
  to 'lev:positive="down"' if we are regridding to a GCHP restart file.
  This will prevent the 'lev:positive="up"' attribute from GEOS-Chem
  Classic restart files from being carried over to the output GCHP
  restart files.
- Trim trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
CHANGELOG.md
- Add note about setting lev:positive="down" in regrid_restart_file.py

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Import is_cubed_sphere_diag_grid from cstools.py
- Now pass weightsdir as a command-line argument, and pass it to
  the functions that create the regridders
- Now compute cs_res_in properly with get_cubed_sphere_res_in
- Added verbose printout
- Removed obsolete variables
- Specify NETCDF4 format for the xarray.Dataset.to_netcdf() function
- Renamed "ds" -> "dset" where expedient
- Added functions to save metadata
- Added function to update GCHP grid and coordinates
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/regrid.py
- Restore unravel_checkpoint_lat and ravel_checkpoint_lat functions
  as they were before recent commits.

gcpy/util.py
- Restore original functionality of reshape_MAPL_CS
- Remove multi_index_lat keyword, as this routine is no longer
  called from file_regrid.py.  It is only needed for plotting and
  will be abstracted to the plotting code in a future commit.

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Simplified the top-of-file Pydoc header
- Remove the import xesmf as xe; it's not used.
- Remove the test for xesmf 0.2.1 and later (we use 0.5.1)
- Remove vert_params_out argument from file_regrid routine and
  ancillary subroutines
- Instead of ds_in and ds_out, many routines a single variable dset
  (this should prevent duplicating memory)
- Remove calls to "drop_lon_and_lat"
- Do not call get_vert_grid to get the vertical dimensions, as we are
  only regridding horizontally.
- Added new routine flip_lev_coord_if_necessary, which flips the
  vertical levels and adjusts the lev:positive attribute accordingly.
  The resetting lev:positive has now been removed from the routine
  save_ll_metadata.
- Update comments, trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/cstools.py
- Function get_cubed_sphere_res was returning the Xdim coordinate
  (which is an array value) rather than the Xdim dimension (which is
  a scalar) as the cubed-sphere resolution (for the case when the
  grid is a checkpoint grid).  This was causing upstream errors
  in boolean tests.  Now fixed.
- In routine is_cubed_sphere_rst_grid, make sure "lon" exists as
  a dimension before referencing it.  This had caused some upstream
  errors.

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Use a single "dset" variable instead of "dset_in" and "dset_out",
  which should save memory.
- Now use the numpy.array_equal() function to determine if the
  sg_params_in and sg_params_out are equal.  The previous code used
  a list comprehension iterating over a zip object, which is overkill.

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/grid.py
- Added new function "get_lev_and_ilev", which returns the GEOS-Chem
  eta values at the vertical grid midpoints and edges.  These are
  needed to specify the GEOS-Chem Classic "lev" and "ilev" netCDF
  coordinate variables.
- Added Pydoc header comment at the top of the module.
- Trimmed trailing whitespace

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/grid.py
- Split get_lev_and_ilev into 2 separate functions: get_ilev_coord
  and get_lev_coord.  This facilitates separate calling.
- Trimmed trailing whitespace

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Import get_ilev_coord, get_lev_coord from gcpy/grid.py
- Now get the original data type of the dataset by using the last
  data variable (should skip the GCHP cubed-sphere variables)
- Add a workaround so that the "lev" dimension won't be overwritten
  after restoring the original data precision
- Drop "lons" and "lats" when regridding the to lat/lon "classic" grid
- Make sure "lev" and/or "ilev" is in the dataset before adjusting
  any attributes
- If regridding from CS/SG -> LL, restore GEOS-Chem Classic eta
  values in the "lev" and/or "ilev" coordinates
- Set variable attributes for "lev" and "ilev" in save_ll_metadata

CHANGELOG.md
- Updated accordingly

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/grid.py
- In routines get_ilev_coord and get_lev_coord:
  - Now return a numpy.ndarray object instead of an xarray.DataArray,
    as it is simpler to flip the dimensions this way
  - Added a "top_down" boolean argument to determine when to
    return the coordinate dimension arranged from the top-of-atm
    to the surface.
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Move test for similar CS/SG grids into regrid_cssg_to_cssg.  This
  allow us to skip regridding when the cubed-sphere grids are similar,
  but still allows us to reformat between checkpoint & diagnostic
  dimension format.
- Add workaround to save the "lon" coord values in CS/SG checkpoint
  format as a np.linspace rather than longitude values.  This matches
  the format of restart files created by MAPL.
- Remove "lons" and "lats" variables from CS/SG checkpoint output
- Remove "lon" and "lat" variables from CS/SG diagnostic output
- Update calls to get_lev_coord and get_ilev_coord from regrid.py
  with the appropriate setting of top_down.
- Converted ' quotes to " quotes in strings
- Trimmed trailing whitespace

gcpy/regrid.py
- Updated Pydoc at top-of-module
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Add fixes so that all combinations of regridding should now work
  properly (pending verification)
- Make sure that the proper coordinate values are saved out for
  each type of file
- Abstracted if block into routine order_dims_time_lev_lat_lon
- Rewrote adjust_cssg_grid_and_coords to work for both checkpoint
  and diagnostic grids
- Updated comments and Pydoc
- Trimmed trailing whitespace

gcpy/grid.py
- Added missing "top_down" argument to function get_ilev_coords

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Add modifications in routine "adjust_cssg_grid_and_coords" so that
  when regridding from checkpoint to diagnostic, the "lons" and "lats"
  coordinates are preserved in the output.
- Also pass dim_format_in to adjust_cssg_grid_and_coords in order
  to decide which variables to copy to "lons" and "lats"
- Remove error trap that is no longer needed
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
@yantosca yantosca added topic: Regridding Issues pertaining to horizontal & vertical regridding category: Bug Fix Fixes a bug that was previously reported labels Oct 5, 2023
@yantosca yantosca added this to the 1.4.0 milestone Oct 5, 2023
@yantosca yantosca requested review from lizziel and msulprizio October 5, 2023 17:35
@yantosca yantosca self-assigned this Oct 5, 2023
gcpy/cstools.py
- Added inquiry function to determine if a GCHP grid is
  index top-down or bottom-up (because the lev:positive
  attribute is currently not reliable).

Signed-off-by: Bob Yantosca <[email protected]>
gcpy/file_regrid.py
- Now import is_gchp_lev_positive_down from cstools.py and use it to
  determine if the lev coord needs to be flipped before saving to
  checkpoint format
- Fixed incorrect logic an an if block

Signed-off-by: Bob Yantosca <[email protected]>
Copy link
Contributor

@lizziel lizziel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a partial review so far. I will get to the rest tonight or tomorrow morning.

gcpy/cstools.py Outdated Show resolved Hide resolved
gcpy/cstools.py Outdated Show resolved Hide resolved
gcpy/cstools.py Show resolved Hide resolved
gcpy/cstools.py Outdated Show resolved Hide resolved
gcpy/cstools.py Outdated Show resolved Hide resolved
@lizziel
Copy link
Contributor

lizziel commented Oct 10, 2023

There is definitely something weird going on with CS comparison plots, probably related to this issue: #167. Regridding both C24 and C48 to lat-lon would probably be better for validation. Regardless, that stretched grid plot looks strange, but maybe a separate issue.

gcpy/cstools.py
- Updated the algorithm in function is_gchp_lev_positive_down so that
  it will only identify emissions variables as variables starting
  with "Emis" (but not having "Emis" anywhere else in the variable name).
- Updated Pydoc headers following suggestions by @lizziel
- Trimmed trailing whitespace

Signed-off-by: Bob Yantosca <[email protected]>
Copy link
Contributor

@lizziel lizziel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for all this work @yantosca! I added comments and requested a few minor changes, mostly in comments. I created a feature request to expand on this work here: #267.

gcpy/cstools.py Show resolved Hide resolved
gcpy/regrid_restart_file.py Show resolved Hide resolved
gcpy/regrid.py Outdated Show resolved Hide resolved
gcpy/file_regrid.py Show resolved Hide resolved
gcpy/file_regrid.py Show resolved Hide resolved
gcpy/file_regrid.py Outdated Show resolved Hide resolved
gcpy/file_regrid.py Outdated Show resolved Hide resolved
gcpy/file_regrid.py Show resolved Hide resolved
gcpy/file_regrid.py Outdated Show resolved Hide resolved
gcpy/file_regrid.py Show resolved Hide resolved
@yantosca
Copy link
Contributor Author

Thanks @lizziel. I will address your comments. In the meantime I made the comparison plots of c48 -> 2x25 and c24-> 2x25. Regridding looks OK there.
c24 -> 2x2.5
Screenshot 2023-10-11 at 11 32 46 AM

c48 -> 2x2.5 (note: Ref is mislabeled but the resolution is OK)
Screenshot 2023-10-11 at 11 34 10 AM

gcpy/cstools.py
- Added a TODO statement to remind us to rethink the logic
  in function is_cubed_sphere_rst_grid if changes are made to
  GCHP restart files/internal state variables

gcpy/file_regrid.py
- Added TODO comments in several routines
- Updated Pydoc headers to be more correct
- In flip_lev_coord_if_necessary, add a call to is_gchp_lev_positive_down
  to flip "lev" and "ilev" coords when "dim_format_out" is "classic"
- Updated comments, trimmed trailing whitespace

gcpy/regrid.py
- Remove commented-out code

gcpy/regrid_restart_file.py
- Added a TODO comment to make the flipping of levels more robust

Signed-off-by: Bob Yantosca <[email protected]>
@yantosca yantosca requested a review from lizziel October 11, 2023 19:07
@yantosca yantosca merged commit 7949550 into dev Oct 11, 2023
@yantosca yantosca deleted the bugfix/lev-positive-attribute-in-restart-regridding branch October 11, 2023 19:54
yantosca added a commit that referenced this pull request Oct 11, 2023
This merge brings the feature/plot-subdir up-to-date with the dev
branch (which now has the fixes for regridding that were added
by PR #266).

Signed-off-by: Bob Yantosca <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
category: Bug Fix Fixes a bug that was previously reported topic: Regridding Issues pertaining to horizontal & vertical regridding
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[QUESTION]Selection of interpolation method in ESMF_RegridWeightGen
3 participants