Fix mpas mesh tools

[hollyhan]

Hi Xylar and Matt,

In discussing with @matthewhoffman, we found the need to make some changes to two mesh tool scripts so that the mesh/scripfiles generated from these tools have positive longitude convention [0 2pi], which is required by MPAS mesh spec (https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf)

Your review on these changes would be very appreciated.
Thank you!
Holly

[xylar]

Thanks for the PR! This is a great thing to fix. I have a few suggestions that we can discuss when we meet next, or here if you prefer.

[matthewhoffman]

Thanks for making this PR, @hollyhan . It will be worth it in the long run to have taken the time to fix these issues before we forget about them.

I looked over the changes, and I don't have anything to add beyond what @xylar identified. All of his suggestions are good ones.

[xylar]

Other than the two super nit-picky comments I just put in, this looks great! I just want to run a quick test, then I'll approve.

[xylar]

Hmm, my tests weren't successful. I did this to create a test conda environment:

conda create -y -n test_mesh_tools --file conda_package/dev-spec.txt 
conda activate test_mesh_tools
python -m pip install -e conda_package/

I verified that it worked (that it has the expected entry point):

scrip_from_mpas --help

But when I ran it on a file we use in COMPASS, it didn't work:

Traceback (most recent call last):
  File "/home/xylar/miniconda3/envs/test_mesh_tools/bin/scrip_from_mpas", line 33, in <module>
    sys.exit(load_entry_point('mpas-tools', 'console_scripts', 'scrip_from_mpas')())
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 172, in main
    scrip_from_mpas(options.mpasFile, options.scripFile, options.landiceMasks)
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 51, in scrip_from_mpas
    lonCell = numpy.mod(lonCell, 2.0*np.pi)
NameError: name 'numpy' is not defined

It needs to be np instead of numpy. That was my mistake in my original suggestion, so I'm sorry about that.

With this fix, I was able to generate a script file that was identical to the one generated by the old scrip_from_mpas.

Similarly, when I tried to run set_lat_lon_fields_in_planar_grid.py, I see:

Traceback (most recent call last):
  File "/home/xylar/miniconda3/envs/test_mesh_tools/bin/set_lat_lon_fields_in_planar_grid.py", line 7, in <module>
    exec(compile(f.read(), __file__, 'exec'))
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mesh_tools/planar_grid_transformations/set_lat_lon_fields_in_planar_grid.py", line 88, in <module>
    lonCell = numpy.mod(lonCell, 2.0*np.pi)
NameError: name 'numpy' is not defined

This time, numpy isn't getting imported at all, so you need to change numpy to np again and add import numpy as np at the top of the file.

When I made this change, I was able to add latitude and longitude to the mesh for the compass dome test case:

set_lat_lon_fields_in_planar_grid.py -f mpas_grid.nc -p ais-bedmap2

Things looked as expected, though this wasn't a very rigorous test because longitudes are all between 0 and pi/2...

[xylar]

@hollyhan, I retested tested with the changes you made. I had to go to a little extra trouble to get the C++ tools without actually building the conda package:

conda create -y -n test_mpas_tools python=3.9 mpas_tools=0.12.0
conda activate test_mpas_tools

That creates an environment with the latest release of mpas_tools, including the C++ tools like MpasCellCuller.x.

Then I install over top of it with the changes you made:

python -m pip install -e conda_package/

This is a little risky in general for reasons I'm happy to explain (dependencies may have changed since the last release, for example) but it's good enough for testing.

Then, I created a mesh to test with:

# create a planar mesh
planar_hex --nx 30 --ny 34 --dc 2000.0 --npx --npy -o test.nc
# cull the periodic boundaries
MpasCellCuller.x test.nc culled.nc
# center the mesh at (0,0)
translate_planar_grid -f culled.nc -c
# add lat/lon
set_lat_lon_fields_in_planar_grid.py -f culled.nc --proj ais-bedmap2
# make sure lonCell are all positive:
ncdump -v lonCell culled.nc

This all worked great.

Then, I tested the SCRIP file generation:

$ scrip_from_mpas -m culled.nc -s scrip.nc
== Gathering information.  (Invoke with --help for more details. All arguments are optional)
 -- Landice Masks are disabled

Traceback (most recent call last):
  File "/home/xylar/miniconda3/envs/test_mesh_tools/bin/scrip_from_mpas", line 33, in <module>
    sys.exit(load_entry_point('mpas-tools', 'console_scripts', 'scrip_from_mpas')())
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 173, in main
    scrip_from_mpas(options.mpasFile, options.scripFile, options.landiceMasks)
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 54, in scrip_from_mpas
    sphereRadius =  constants['SHR_CONST_REARTH']
NameError: name 'constants' is not defined

This was a surprise because I thought this worked when I tested yesterday. But it seems like we're missing an import here:

from mpas_tools.cime.constants import constants

I added this myself and things worked as expected.

[hollyhan]

Thanks for going through the testing and explaining the steps in detail, @xylar ! Please let me know if there's anything I should do further to get this PR closed.

[xylar]

@hollyhan, as I said above, from mpas_tools.cime.constants import constants needs to be added to mpas_tools/scrip/from_mpas.py. I did this in my local testing but I didn't modify this branch.

Then, you and I need to work together on rebasing this branch to clean thing up into just 2 or 3 commits.

[xylar]

@hollyhan, would you be willing to change from_mpas.py to check that lonCell and lonVertex are in the desired range (0, 2*pi) rather than modifying them? See #450 (comment)

Then, I'll work with you to rebase the PR into 2 or 3 commits as I mentioned last week.

[hollyhan]

@xylar Yes, I'll work on the check for the range of lonCell and lonVertex before we meet on Wednesday.

[xylar]

@hollyhan, so we still missed a step in from_mpas.py. When I test, I see:

$ scrip_from_mpas -m culled.nc -s scrip.nc
== Gathering information.  (Invoke with --help for more details. All arguments are optional)
 -- Landice Masks are disabled

Traceback (most recent call last):
  File "/home/xylar/miniconda3/envs/test/bin/scrip_from_mpas", line 33, in <module>
    sys.exit(load_entry_point('mpas-tools', 'console_scripts', 'scrip_from_mpas')())
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 177, in main
    scrip_from_mpas(options.mpasFile, options.scripFile, options.landiceMasks)
  File "/home/xylar/code/MPAS-Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 58, in scrip_from_mpas
    sphereRadius =  constants['SHR_CONST_REARTH']
NameError: name 'constants' is not defined

Could you add the missing import:

from mpas_tools.cime.constants import constants

and do the following?

git commit --amend from_mpas.py
git push --force hollyhan/MPAS-Tools fix_MPAS_mesh_tools

[hollyhan]

@xylar I've added the missing import and pushed the branch

[xylar]

Yep, @hollyhan, that did the trick.

I have tested this to my satisfaction. Thanks so much for the hard work. It seemed like a simple task but I managed to make it complicated!

[hollyhan]

@xylar, I tested the new from_mpas.py with the thwaites mesh that uses the longitude convention of [-pi pi].
The test results shows the following error, as we expect:

$ scrip_from_mpas -m thwaites.4km.210608.nc -s scrip.nc 
== Gathering information.  (Invoke with --help for more details. All arguments are optional)
 -- Landice Masks are disabled

Traceback (most recent call last):
  File "/Users/hollyhan/miniconda3/envs/test_mpas_tools/bin/scrip_from_mpas", line 33, in <module>
    sys.exit(load_entry_point('mpas-tools', 'console_scripts', 'scrip_from_mpas')())
  File "/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 177, in main
    scrip_from_mpas(options.mpasFile, options.scripFile, options.landiceMasks)
  File "/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 52, in scrip_from_mpas
    raise ValueError("lonCell is not in the desired range (0, 2pi)")
ValueError: lonCell is not in the desired range (0, 2pi)

[hollyhan]

@matthewhoffman,

@xylar and I cleaned up the previous commits by git-rebasing and made new changes to the code - the code now checks for the proper range [0 2pi] of longitude for MPAS scripfiles and raises a value error if the range is outside of [0 2pi] rather than changing the longitude convention to [0 2pi].

[matthewhoffman]

@hollyhan , these changes make sense based on where the discussion ended up. If you haven't already done so, it's probably a good idea to rerun your workflow of adding lat/lon to a mesh file, running the remap tool, and making sure the map file works properly for regridding inside your SLM code in MALI.

[xylar]

@hollyhan, I'll merge this as soon as you have re-tested your workflow and let us know that it worked as expected.

[hollyhan]

Hi @xylar and @matthewhoffman, I have tested the entire workflow and and confirmed that the new workflow passes the test. The following is how I've tested:

1. Try creating a scripfile from a thwaites mesh that uses undesired longitude convention:

$ scrip_from_mpas -m thwaites.4km.220203.nc -s thwaites.4km.220203.scripfile.nc 

== Gathering information.  (Invoke with --help for more details. All arguments are optional)
 -- Landice Masks are disabled

Traceback (most recent call last):
  File "/Users/hollyhan/miniconda3/envs/test_mpas_tools/bin/scrip_from_mpas", line 33, in <module>
    sys.exit(load_entry_point('mpas-tools', 'console_scripts', 'scrip_from_mpas')())
  File "/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 177, in main
    scrip_from_mpas(options.mpasFile, options.scripFile, options.landiceMasks)
  File "/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mpas_tools/scrip/from_mpas.py", line 52, in scrip_from_mpas
    raise ValueError("lonCell is not in the desired range (0, 2pi)")
ValueError: lonCell is not in the desired range (0, 2pi)

Error appears as expected.

2. Change the longitude convention with the new changes made in this PR

$ set_lat_lon_fields_in_planar_grid.py -f thwaites.4km.220203.nc --proj ais-bedmap2 
== Gathering information.  (Invoke with --help for more details. All arguments are optional)

Using ais-bedmap2 projection, defined as: +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
Input file xCell min/max values: -1623967.0094751006 -1073016.3730225994
Input file yCell min/max values: -828184.9435592396 -161464.9911215813
/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mesh_tools/planar_grid_transformations/set_lat_lon_fields_in_planar_grid.py:84: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  lonCell, latCell = pyproj.transform(projections[options.projection], projections['latlon'], xCell[:], yCell[:], radians=True)
/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mesh_tools/planar_grid_transformations/set_lat_lon_fields_in_planar_grid.py:85: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  lonVertex, latVertex = pyproj.transform(projections[options.projection], projections['latlon'], xVertex[:], yVertex[:], radians=True)
/Users/hollyhan/Desktop/E3SM-Project/MPAS_Tools/fix_MPAS_mesh_tools/conda_package/mesh_tools/planar_grid_transformations/set_lat_lon_fields_in_planar_grid.py:86: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
  lonEdge, latEdge = pyproj.transform(projections[options.projection], projections['latlon'], xEdge[:], yEdge[:], radians=True)
Calculated latCell min/max values (radians): -1.3943921358761688 -1.3008508250402886
Calculated lonCell min/max values (radians): 4.119929571080232 4.586072104970466
Lat/lon calculations completed.  File has been written.

3. Create a scripfile with the new mesh that uses the desired longitude convention

$ scrip_from_mpas -m thwaites.4km.220203.nc -s thwaites.4km.220203.scripfile.nc 
== Gathering information.  (Invoke with --help for more details. All arguments are optional)
 -- Landice Masks are disabled

 -- WARNING: sphereRadius<0 so setting sphereRadius = 6371220.0
 -- WARNING: 'on_a_sphere' attribute is 'NO', which means that there may be some disagreement regarding area between the planar (source) and spherical (target) mesh
Input latCell min/max values (radians): -1.3943921358761688, -1.3008508250402886
Input lonCell min/max values (radians): 4.119929571080232, 4.586072104970466
Calculated grid_center_lat min/max values (radians): -1.3943921358761688, -1.3008508250402886
Calculated grid_center_lon min/max values (radians): 4.119929571080232, 4.586072104970466
Calculated grid_area min/max values (sq radians): 4.728403336933561e-07, 2.4751188489024916e-06
Creation of SCRIP file is complete.

4. Create a mapfile going from the mpas mesh to the gaussian grid using ESMF weight generator.

$ ESMF_RegridWeightGen -d gauss_latlon512x1024_rank2_n2s_scripfile.nc -s thwaites.4km.220203.scripfile.nc -w mpas_to_grid.220203.nc -i --64bit_offset  --src_regional

 Starting weight generation with these inputs: 
   Source File: thwaites.4km.220203.scripfile.nc
   Destination File: gauss_latlon512x1024_rank2_n2s_scripfile.nc
   Weight File: mpas_to_grid.220203.nc
   Source File is in SCRIP format
   Source Grid is a regional grid
   Source Grid is an unstructured grid
   Use the center coordinates of the source grid to do the regrid
   Destination File is in SCRIP format
   Destination Grid is a global grid
   Destination Grid is a logically rectangular grid
   Use the center coordinates of the destination grid to do the regrid
   Regrid Method: bilinear
   Pole option: NONE
   Ignore unmapped destination points
   Output weight file in 64bit offset NetCDF file format
   Line Type: cartesian
   Norm Type: dstarea
   Extrap. Method: none

5. Interpolate ice thickness from the thwaites mesh to the global gaussian grid using ncremap

$ ncremap -m mpas_to_grid.220203.nc -i thwaites.4km.220203.nc -o interp_mpas_to_grid_results_220203.nc -v thickness

The resulting interpolated thwaites ice thickness is correct.

[xylar]

@hollyhan, thanks for documenting your testing carefully. It looks good!

[matthewhoffman]

Yes, thanks @hollyhan for documenting the workflow! I'm sure it will be useful to have this for reference in the PR history.

[xylar]

@matthewhoffman, @hollyhan and @trhille, I'll do a release of MPAS-Tools soon and update the version in compass so you can make use of this