Add a tool to create transect cell and edge fields with edge sign and distance

[xylar]

See the new README for how to use it.

[xylar]

@alicebarthel, I think this tool will be helpful for the transport plots you're working on.

It needs to be used in conjunction with MPAS-Dev/geometric_features#202 so we'll need to discuss special instructions for how to install that verson of geometric_features before it gets released. Here's the gist.

Create a new environment with the dependences:

conda create -y -n transects python=3.12 geometric_features matplotlib mpas_tools netcdf4 numpy scipy shapely xarray
conda activate transects

Then, go to the location where you've checked out MPAS-Dev/geometric_features#202

python -m pip install -e .
export GEOMETRIC_DATA_DIR=$PWD/geometric_data

(You'll need to set GEOMETRIC_DATA_DIR each time you want to use this conda environment or it won't find the new transect).

To run this code, go into the subdirectory:

cd ocean/transects/python

Then, you will need to find the initial condition you want, e.g.:

ln -s /lcrc/group/e3sm/public_html/inputdata/ocn/mpas-o/SOwISC12to60E2r4/ocean.SOwISC12to60E2r4.230220.nc

and then run:

./create_transect_masks.py -m SOwISC12to60E2r4.230220.nc -f "Antarctic isobath at -1000.0 m" -o SOwISC12to60E2r4_Antarctic_isobath_at_-1000.0_m.nc --process_count 4

You could safely run with more processes if you run on a compute node. The resulting fields on cells and edges should be useful masks, distances and edge sign (on edges, of course).

As an example, here's the distance on edges:

[milenaveneziani]

@xylar: this immediately caught my attention, since I've been doing lots of transects and regional open boundary calculations lately.
@alicebarthel: do you know how this is different from the workflow I usually use (get mask through MpasMaskCreator and find edges of transects to plot/compute quantities)?
Thanks!

[xylar]

@milenaveneziani, the main added value here is that we needed distance along the transect at each edge where the edge sign is non-zero. I don't think we have a way to compute that anywhere else so that's the added value here.

[xylar]

get mask through MpasMaskCreator and find edges of transects to plot/compute quantities

The edge mask and edge sign is found with the python version of the mask creator but it's otherwise pretty similar.

The distance along the transect is computed first at every place that a transect segment intersects with the edges of an MPAS cell (something we already had for MPAS-Analysis), then those distances are averaged back to cells or edges based on all the segments that cross the cell or the diamond shape associated with an edge.

[xylar]

To me, this part is the important new addition.

The transect module is just here because the version in Polaris (which this is a copy of) is cleaner than the version in MPAS-Tools but we don't have a way to use it in QuickViz because Polaris isn't a simple conda package that can be installed from conda-forge.

[alicebarthel]

yes, @milenaveneziani I did a first-order analysis by combining the Antarctic shelf regions and using the openBoundaries script. Ordering the edges by longitude was not completely clean around the peninsula, and I wanted to expand the analysis to plot by distance along the edge. In addition, I had to check that the recombining of regions is continuous across the 2 resolutions (LR and SORRM) I've done analysis over. @xylar offered to develop this transect (PR#202) to have a reference 1000-m circum-Antarctic transect, with the distance calculated using the existing tools.
From today's conversation, @xylar suggested adding the relevant circum-Antarctic shelf regionmask too. I think that will be great to have, esp. with it being neatly consistent with the transect.

[milenaveneziani]

Thank you both. This is great. I think the reason why you still need to go through MpasMaskCreator is that that tool gives you sorted edges (from point A to B): is that correct?
And yes: having a circumantarctic transect along the shelf sounds great. I would go further and make it into the offshore boundary of a circumantarctic shelf region. So, basically adding two features in geometric_features.

[xylar]

@milenaveneziani, yes, we also decided that a feature defining the enclosed region would be useful and I did that in:
MPAS-Dev/geometric_features#202

@alicebarthel and @cbegeman, I made the requested changes to MPAS-Dev/geometric_features#202.

I will add a script tomorrow that can be used to cut out a square from a closed-loop transect (centered at a given lat/lon and with a given size). I'm too tired to do that tonight even though I wanted to keep the momentum going...

[xylar]

Okay, I added the tool for cutting the loop, as promised. The updated instructions to replace the last step in are:

./cut_closed_transect.py -f "Antarctic 1000m isobath" --lat -65.91 --lon 80.38 --size 0.1 -o cut_Antarctic_1000m_isobath.geojson
./create_transect_masks.py -m SOwISC12to60E2r4.230314.nc -g cut_Antarctic_1000m_isobath.geojson -o SOwISC12to60E2r4_Antarctic_1000m_isobath.nc --process_count 4

the instructions before that are all the same.

You may also want to look at the resulting fields in paraview:

 paraview_vtk_field_extractor.py -v all -m SOwISC12to60E2r4.230314.nc -f SOwISC12to60E2r4_Antarctic_1000m_isobath.nc --ignore_time -d nTransects=0
paraview vtk_files/staticFieldsOnEdges.vtp

[xylar]

@alicebarthel, could you test this out on the meshes you plan to analyze (in particular the LR mesh)? It would be important to know if there are gaps in the transect, following @milenaveneziani's concerns. Even if so, we might be okay with that since these are presumably places where there simply isn't any continental shelf according to the model.