Update cice.t-test.py to use cartopy instead of basemap.

[daveh150]

For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium,
please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers

PR checklist

• Short (1 sentence) summary of your PR:
  Update cice.t-test.py to use cartopy instead of basemap.

• Developer(s):
  David Hebert

• Suggest PR reviewers from list in the column to the right.

• Please copy the PR test results link or provide a summary of testing completed below.

• How much do the PR code changes differ from the unmodified code?

  • [ x] bit for bit
  • different at roundoff level
  • more substantial

• Does this PR create or have dependencies on Icepack or any other models?

  • Yes
  • [x ] No

• Does this PR add any new test cases?

  • Yes
  • [ x] No

• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/. A test build of the technical docs will be performed as part of the PR testing.)

  • [ x] Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • [ x] No

• Please provide any additional information or relevant details below:

[eclare108213]

@daveh150
These are fairly extensive changes -- it took some work to do this. Thanks! Do you have sample output, before and/or at least after?

For future reference (and you can edit your text above) - deleting the space from the checkbox code when you add an x will make them work... [x] instead of [ x] etc.

[daveh150]

Thank you. The contour plot section was the challenging part, needing both +/- 180 and 0-360 lat/lons.

I admittedly did not test the plot in this exact script, I do not have a perturbation run to test with (I likely missed how to set that up on the testing documentation). I was testing plots outside of this script using output from one the 'smoke' tests. If anyone has suggestion on how to make a perturbation test I can test this.

[daveh150]

Oh, I think I see why I couldn't run a QC test before. Previously I got a message that my test and baseline were bit-for-bit, so no qc test was required. But I did not make any changes to the ice_in, so the runs were identical. I'll try to change the shortwave and ktherm in the 'test' run and see if I can run the QC test to get example plots from this script. Sorry about that.

[daveh150]

Does anyone use the contour option in these plots? That is by far the most challenging, the scatter and pcolor seem to do what we want. I'd like to just remove the contour plots if they aren't used.

[eclare108213]

I didn't know that contour plotting was even available. As far as I'm concerned, color maps should be sufficient and you can comment the contouring part out. (Am I correctly understanding the issue here?) Check the docs to see if it's mentioned, and if so then comment that out too. Thx

[apcraig]

This looks fine to me, I'll try to run a t-test on cheyenne and make sure everything is working to independently confirm. I've struggled at times getting the python environment setup, has anything changed in that realm?

[daveh150]

@apcraig The only update to your python environment I can think of is the need for the cartopy package.

Here are results I obtained with gx3 case using the scatter and pcolor options.

The contour plot is below, I am having a very difficult time getting rid of the space in the Northern Hemisphere plot. I think it is related to contour lines crossing +/- 180 but I can't seem to fix it.

[apcraig]

Thanks @daveh150, those plots look pretty good overall. No idea about the missing data in the contour plots.

[daveh150]

I made an update where if contour is selected, the plots will revert to pcolor. Eventually that whole section of the code could be removed once we are sure no one wants to use it. It cleans up the plotting code significantly. I did not see a need to add the cyclic arrays with the pcolormesh.

[phil-blain]

I took a quick look to the changes and this looks good, thanks @daveh150 !

One thing I would suggest is to add cartopy to the conda environment specification file (environment.yml under configuration/scripts/machines) and remove basemap from there.

[dabail10]

Just got a chance to look at this. This is great! Thanks @daveh150 for doing this. Here are the plots from my gx3 B vs C grid comparison. A couple of notes for enhancements later.

1. This does seem to be in the right sense of the differences I would expect. That is the B grid solution is thinner.
2. However, It would be nice if the colorbar could be reset a bit to bring out the differences. The "red" area is basically zero everywhere and so a red-white-blue colormap might be nice for the difference plot. Then instead of having the maximum distance as the colorbar, maybe a one-sigma difference is better. That is computing the spatial variance of the distances and using that as the max on the colorbar instead of the absolute max difference. Perhaps even two-sigma. This would help bring out the smaller differences better.
3. Also, it would be nice if the colorbar for the two cases was the same. This way one could see the differences pretty clearly in these plots. Or we could always use a max of 5m in the NH and 2m in the SH. Similarly in the difference plot, we could always use a min/max of -1/1m.

[dabail10]

This looks great. I added some additional comments for enhancements.

[dabail10]

I made an update where if contour is selected, the plots will revert to pcolor. Eventually that whole section of the code could be removed once we are sure no one wants to use it. It cleans up the plotting code significantly. I did not see a need to add the cyclic arrays with the pcolormesh.

I agree that it is probably not worth chasing down the issues with the contour plots. Typically we don't see issues like this except on tripolar grids. It also looks like the line of missing values goes along the 90E/W line. It would make more sense if it was related to 180/0. I recall that Matt had designed it so that one could choose in case the pcolor doesn't work for some reason. What happens if you don't use filled contours?

[apcraig]

@dabail10, thanks for reviewing and testing the scripts. @daveh150, Are we ready to merge? It seems like the outstanding requests may be

• update environment.yml, add cartopy, remove basemap
• update the colorbar

I think those are two good suggestions. If we prefer to defer, I'm also OK merging.

[apcraig]

@daveh150, just checking in again to see how we should proceed.

[daveh150]

Sorry for my delays, I’ve been out the past 2 weeks on military leave. I am catching up now and should be able to follow up tomorrow on how to proceed to commit this to the repository. Thanks! David From: Tony Craig ***@***.***> Sent: Tuesday, June 21, 2022 14:49 To: CICE-Consortium/CICE ***@***.***> Cc: David Hebert, Code 7322 ***@***.***>; Mention ***@***.***> Subject: Re: [CICE-Consortium/CICE] Update cice.t-test.py to use cartopy instead of basemap. (PR #727) @daveh150 <https://github.com/daveh150> , just checking in again to see how we should proceed. — Reply to this email directly, view it on GitHub <#727 (comment)> , or unsubscribe <https://github.com/notifications/unsubscribe-auth/AE52VPBLCKGNZPFFUYFLXBLVQIMBZANCNFSM5XWBT6JQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[daveh150]

Thanks for your patience while I catch up from being out the past two weeks.

@dabail10 I took you suggestion and added individual colormaps to each plot. For the NH/SH data I hardwired the limits to 0-5 in NH, 0-2 in SH. happy to change that if needed. Also, for the difference plot I added blue/white/red colormap and determine the range from the absolute value of the data. Could you please see how this works for you?

@phil-blain I modifies environment.yml to remove basemap and add cartopy.

[apcraig]

Thanks @daveh150. @dabail10, do you have some QC data you can run the updated scripts on to verify they work and look good? If not, I will try to generate some data and do the same. thanks!

[dabail10]

I do have some runs I can test this on.

[dabail10]

These are great! It is great to see such a large area of white where the differences are close to zero. We might even be able to set the min/max on the difference plots to even smaller, but this is great for now.

[phil-blain]

That might be something on my side, but with the script in this PR I get plots like this;

I haven't yet looked into why it's all blue...

[phil-blain]

Suggested change

	print("contour plot depreciated. using pcolor.")
	print("contour plot deprecated. using pcolor.")

[dabail10]

This all looks good to me. Not sure why Phillipe is getting blue everywhere. Is the colormap slightly different on your machine?

[phil-blain]

It should not be... I'm looking into it ....

[phil-blain]

This is with Cartopy 0.20.2, Matpltlib 3.5.2. What versions are you using ?

[dabail10]

Cartopy 0.20.2 matplotlib 3.5.1

[phil-blain]

I can investigate on my own as to not delay merging this.

[phil-blain]

my plots are on the gx1 grid, maybe that's why I'm seeing something different?

[dabail10]

Just ran it on one of my gx1 tests and I don't see the blue.

[phil-blain]

In my plot I've just noticed there is also a dark longitudinal line near 30°E...

[phil-blain]

I've spent more time debugging this. I think it's due to the fact that the t_lat and t_lon arrays are masked arrays due to land proc elimination:

t_lat:

t_lon:

If I do something like t_lat=t_lat.filled(0), and the same for t_lon, i.e. filling the masked values with zeros:

Then I get the expected plot:

EDIT the doc for pcolormesh does state that masked arrays is only supported for the data itself, not the coordinates arrays...

EDIT 2 and if we try to do a plain matplotlib pcolormesh (i.e. without cartopy), we get an error saying the coordinate arrays can't be masked arrays. So it's really cartopy here that seems to fail to check its inputs...

---

So @dabail10 I guess my next question would be what was your MPI x OpenMP configuration for the gx1 runs you used to test the new script ? Would it be possible that there was no land procs?

[phil-blain]

So I zoomed out a bit in my plot above and I don't think it's as easy as filling the lat/lon arrays with zeros:

we get spurious polygons, with different shapes depending on the filling value. I think we should use "complete" i.e. non masked tlat, tlon. This leads to the correct plot:

So one way to do it would be for the script to read the location of the input data directory (CICE_data) from the env.machine_compiler file and then read the grid from there to get the full lat/lon arrays, The only "problem" is that we currently do not have the gx1 grid as a NetCDF file.

[daveh150]

@phil-blain @dabail10 if the issue is related to masked locations in the tlon/tlat but NOT the data, could you try this to get the mask the same. Then I think the pcolor will ignore the masked locations:

data.mask = np.logical_or(data.mask,lon.mask)

[phil-blain]

I just tried that and it does not work (gives me the blue plot again). I think that we should use the full lat/lon arrays.

[eclare108213]

In my opinion, tlat and tlon ought to be made 'complete' (no eliminated land blocks) in the output files. One way to do this would be to write the original, global arrays, rather than writing them after the decomposition when the blocks are eliminated. Not sure how easy this would be, or if it would work for PIO etc.

[daveh150]

Thanks @phil-blain. I think you had an idea to read the gx1 grid file instead of from the history file, if it was netcdf? That should be straight forward, we just need to write the TGRID in addition the UGRID.

[daveh150]

I made a grid_gx1.nc file. Hopefully it is attached, can y'all verify it is okay then perhaps add it to the CICE_data grid files for future use?
grid_gx1.nc.gz

Updating the cice.t-test.py to find this file will take a bit of thought, since it is not in the run directory.

[daveh150]

Sorry for quick modification of the gx1 netcdf, I changed the variables to be all lowercase in accordance with popgrid_nc. Also changed the kmt data type to integer and unit to 'none'.

grid_gx1.nc.gz

[phil-blain]

I implemented what we talked about in the previous meeting i.e. disabling land blocks elimination.

@daveh150, you could update your branch to include my commit thusly:

git checkout update_plots
git pull git@github.com:phil-blain/CICE.git update_plots
git push

[apcraig]

Thanks for doing that @phil-blain, I think I volunteered but then dropped the ball. @daveh150, could you pull @phil-blain's change to your sandbox and retest? @phil-blain, have you confirmed that this fixes the issues you were seeing? If so, and everyone agrees, we can problem merge soon. Thanks!

[phil-blain]

I had already verified that using the full TLON/TLAT arrays made the script work correctly. I had done that by creating a gx1 test with -p 1x1, which then changed the block size and resulted in a single block and thus no land blocks elimination. I then manually loaded these TLON/TLAT arrays in Python and used those to create the "correct" maps in #727 (comment).

I just verified now that with this updated qc option, i.e. including distribution_wght = 'blockall', creating a test with -p 44x1 -s qc , i.e. what the doc suggests using, the land blocks are indeed not eliminated and the TLON/TLAT arrays in the output files do not have any missing values. So I think this shows the fix is sufficient :)

[apcraig]

Thanks @phil-blain, it makes sense the blockall would fix the problem. Thanks for confirming.

@daveh150, if there are any mods related to pulling in a separate lon/lat file for plotting, those can be removed again.

@daveh150, @phil-blain's mod is a one line change. If you prefer to implement it manually rather than pull/push, that would be fine too. You can see the mod here, daveh150/CICE@update_plots...phil-blain:CICE:update_plots, as noted by @phil-blain.

[apcraig]

See #742

[daveh150]

I've been on leave since last week. I did not have any other changes for the lat/lon issue, glad that the blockall solution by @phil-blain worked! Thank you.