Fix AntarcticIceCoverage regions

[xylar]

The regions have been extracted directly from the original masks on
the Bedmap2 grid, rather than after being projected onto a lon/lat
grid. These regions also previously had self-intersections thats
made them invalid for certian shapely operations.

The masks were first turned into a signed-distance functions and
then smoothed slightly (gaussian filter with sigma of 0.5 grid
cells) before finding the zero-contour.

[xylar]

@pwolfram, I noticed when I was trying to remove land from ocean domains that the previous feature for Antarctic ice coverage (either just grounded ice or both grounded and floating ice) included some self-intersecting polygons and was therefore "invalid" according to shapely. This created problems for certain operations like combining this geometry with the land mask from the rest of the ocean.

I have remedied this by recomputing the shapes directly from the Bedmap2 data set, with some intermediate smoothing but no intermediate projection onto a lon/lat grid.

Here is a Gist where I have the scripts for 1) exporting the Bedmap2 geometry including lon/lat to a NetCDF file and 2) for computing the features based on a mild smoothing of the the Bedmap2 masks.

https://gist.github.com/xylar/7a140936ec7c3bc7be0b6fdba19e7771

Please let me know if you have any concerns and merge at your convenience.

[xylar]

@mark-petersen, these changes should be tested on one or more global_ocean test case with land-ice cavities before getting merged. I don't anticipate any complications but it would be best to be sure.

So @pwolfram, go ahead and take a look at this but don't merge it just yet.

[pwolfram]

@xylar, should we be storing the processing scripts in this repo? I think we probably should be putting them in driver_scripts to ensure that the whole workflow is preserved. The downside here is that instructions on obtaining the external data source will be required. The script may also need to be documented to clearly designate its example usage in the script doc string too.

[xylar]

@pwolfram, yes, I agree that's a good idea. I think these probably belong somewhere other than driver scripts to keep things from getting confusing. I'll think about it and put them in another folder. I agree that the docstring needs to include instructions on downloading Bedmap2 in this case.

Obviously hold off on merging this until I can make this update and also until we can make sure this does no harm to global_ocean test cases both with and without land ice.

[xylar]

@pwolfram, see what you think of the docstrings for the two scripts I just added. I don't think you need to be as rigorous with the coding practices or general usability of scripts going into the feature_creation_scripts directory because these will not be intended to be modified or part of the day-to-day workflow of users. For example, I didn't provide a mechanism for pointing to the Bedmap2 dataset at the command line. Even so, their purpose should at least be clear and it should be possible to recreate the associated features if desired.

[pwolfram]

@xylar, the docstrings are good. Should I test the workflow described in them, however? It seems like if I can get identical results this should be good. Would you agree?

[xylar]

@pwolfram, you can do that but it seems like overkill to me. The data set is several GB and the processing isn't super fast. So I definitely invite you to do so but it seems like maybe a waste of your time. More important is that the results are useful, and that's something that Mark and I will need to test over the next few days.

[pwolfram]

Thanks for the info @xylar that I don't necessarily need to test these scripts. Thanks for including them in the directory, however, to fully document the work. I think feature_creation_scripts is a good name for the directory.

[xylar]

This PR still needs to be tested with global_ocean init_step1 test cases in MPAS-O. It should also wait until #52 gets fulfilled because I have rebased to include those changes. The reason for doing so was that the scripts in create_feature_scripts/extract_bedmap2_ice_coverage could be updated with the new syntax.

[xylar]

Testing

I have now tested this branch on init_step1 of a number of MPAS-O global_ocean test cases:

QU_240km/performance_test
QU_240km/restart_test
QU_240km/se_blocks_test
QU_240km/rk4_blocks_test
QU_240km/analysis_test
QU_240km/with_land_ice
QU_120km/with_land_ice
EC_60to30km/spin_up
EC_60to30km/with_land_ice
RRS_30to10km/default
RRS_30to10km/with_land_ice

Since the mask is different, the results are not bit-for-bit compared with using the Antarctic land coverage from master (except for the QU_240km test cases without land-ice cavities, which coincidentally have identical masks to those produced with master). However, the changes are quite small as the following images show, and are not expected to adversely affect simulation results.

EC_60to30km without land-ice cavities

master

this branch

EC_60to30km with land-ice cavities

master

this branch

RRS_30to10km without land-ice cavities

master

this branch

RRS_30to10km with land-ice cavities

master

this branch

[xylar]

cc: @mark-petersen, @toddringler, @jonbob, @milenaveneziani, @pwolfram, @maltrud, @vanroekel

This merge would mean that initial conditions produced in MPAS-O, both with and without land-ice cavities, would change. Importantly, they would have a different number of cells (and edges and vertices). Presumably, this means that this PR can only be merged after initial conditions for ACME v1 are finalized. Would you agree?

I think this merge is worthwhile for two reasons: 1) the two AntarcticIceCoverage masks we are currently using has self-intersections, which makes them problematic to work with (e.g. they can't be used to mask out other features like the Southern Ocean), and 2) I have provided a pair of scripts (in feature_creation_scripts/extract_bedmap2_ice_coverage) that can be used to reproduce the masks directly from the Bedmap2 data set available to the public (i.e. there is now provenance).

I have slightly smoothed the Bedmap2 mask (using a signed-distance function, for those interested in algorithmic details), which has removed the self-intersections that occurred with the mask interpolated to the ETOPO1 grid that I had previously used.

Once this PR is merged, it would be important that everyone running MPAS-O global_ocean test cases updates to the latest master so we do not have issues where different people are producing different initial conditions from identical MPAS-O code.

[pwolfram]

@xylar, just a quick question- what is the current status of this PR? I'm guessing it has just dropped priority but wanted to make sure I'm not holding up this process.

[xylar]

@pwolfram, this really needs to wait until we have an ACME-v1 branch of this repo, since this change would affect ACME initial conditions. So not waiting on you.

[xylar]

@mark-petersen, when you have some time, I'd like your help thinking about this. The problem remains that I raised last October. I have a fix to the Antarctic topography (both with and without ice-shelf cavities) that I'd like to merge to master. It would be good if this were synchronized with creating the next round of initial conditions for ACME (v4 grids?). In fact, it would have been good if this change had already been there when v3 grids were created. But we don't seem to have a workflow in place for this. Maybe something we can discuss in a little over a week when I'm there?

[mark-petersen]

@xylar Thanks, I lost track of this one (obviously). I'm glad you just pinged me on it, as I finally have time to work on PRs for initial conditions. Reproducing the "V3" grids is an issue. Once this is merged, the only way to reproduce the previous horizontal mesh is to point to the geometric features repo just before this merge, correct? It looks like it changes nCells for both standard domains and those with ice shelf cavities.

[xylar]

@mark-petersen, that is correct. We should probably make a tag for the current state of the repo (e.g. acmev3) to indicate that this tag should be used for v3 grids.

[xylar]

@mark-petersen, here's the 1 1/4 year old PR you're thinking of.

[xylar]

@mark-petersen, I just created a tag "E3SMv3grids" that is appropriate for reconstructing the v3 meshes and initial conditions. I would stick with that tag for what you're doing. We can merge this when you're comfortable but I wouldn't use it until you're ready to create some v4 meshes, as I just said.

[mark-petersen]

@xylar I agree it is better to not project to lat/lon coordinates here. I'd like to merge it in to make this next set of RRS meshes, rather than wait, since I'm going through the trouble to create meshes, review them carefully, and make the mapping files. I will make a note in the MPAS repo that the EC60to30V3 was made just before this PR was merged, with the previous hash.

Once we make these meshes, they tend to be around for several years. Let's use the best method, to our knowledge, each time we do it.

[xylar]

Any reason not to call the meshes you create "v4" to indicate that there is an important distinction between these and v3? The Antarctic coastline will be non-negligibly different.

[xylar]

@mark-petersen, I want to remind you of this work. I had delayed merging this into geometric_features because we didn't want to disrupt any work on v3 grids. However, I think with the tag we made (https://github.com/MPAS-Dev/geometric_features/releases/tag/E3SMv3grids), we should now be good to merge this. I'm worried it won't make it into the next generation of grids if we don't.

[mark-petersen]

See testing of this PR, and review at MPAS-Dev/MPAS-Model#169 (comment). It all works.