Add functions for computing and plotting transects

[xylar]

Checklist

• [x] Developer's Guide has been updated
• [x] API documentation in the Developer's Guide (api.md) has any new or modified class, method and/or functions listed
• [x] Documentation has been built locally and changes look as expected
• [x] Testing comment in the PR documents testing used to verify the changes

[xylar]

Testing

This has been used in testing the branch https://github.com/xylar/polaris/tree/add-isomip-plus-init and produced plots like these:

Sigma coordinate with 10 layers and a thin-film region

Z-star coordinate with 36 layers (and no thin-film region)

[xylar]

@cbegeman, I spent about a week rewriting the transect functionality to use quads. There are a lot of good reasons for us to do this, high among them:

• quads are more intuitive
• we can plot contours with labels
• the matplotlib functionality for plotting fields on quads is a lot more robust, mature, and intuitive than that for triangles This functionality will be used extensively for the MPAS-Analysis port so it's definitely time well spent.

For now, I have the functionality here in Polaris rather than in MPAS-Tools. This is partly to make it easier to test and make bug fixes. It is also partly to prevent breaking anything that relies on the current MPAS-Tools implementation. Since it's a big change to the API, I don't want to put it into MPAS-Tools until everything that relies on the old API (probably just MPAS-Analysis) has either been updated or declared deprecated.

I will post some new plots with this approach shortly. You can try it out on your case(s) as well. I believe it will fix the issue you pointed out to me on Slack but I haven't tried your branch out myself yet.

[cbegeman]

@cbegeman, I spent about a week rewriting the transect functionality to use quads.

Awesome! I'm excited to try this out. I'm having a little bit of trouble because it seems like interp_mpas_to_transect_cells is looking for nVertLevels in ds_transect but my ds_transect does not have that dimension. Can you point me to a branch and task where you have the viz working? I can try to figure it out from there.

[xylar]

I'm having a little bit of trouble because it seems like interp_mpas_to_transect_cells is looking for nVertLevels in ds_transect but my ds_transect does not have that dimension.

Two questions:

• Did you get your ds_transect from the compute_transect() routine here or from something prior to this update?
• Are you trying to call interp_mpas_to_transect_cells() yourself or letting plot_transect() do it for you?

Here is an example of me using compute_transect(): https://github.com/xylar/polaris/blob/add-isomip-plus-init/polaris/ocean/tasks/isomip_plus/init/init.py#L306-L307 and here are 2 different ways I'm using plot_transect(): https://github.com/xylar/polaris/blob/add-isomip-plus-init/polaris/ocean/tasks/isomip_plus/init/init.py#L399-L403 https://github.com/xylar/polaris/blob/add-isomip-plus-init/polaris/ocean/tasks/isomip_plus/init/init.py#L448-L452

[xylar]

I'm having a little bit of trouble because it seems like interp_mpas_to_transect_cells is looking for nVertLevels in ds_transect but my ds_transect does not have that dimension.

Taking a look at the code: https://github.com/E3SM-Project/polaris/blob/c1efd85773b501880019148977ea7600e4a5b0dc/polaris/ocean/viz/transect/vert.py#L202-L231 it seems like maybe you're passing something for da that doesn't have the expected dimensions (nCells, nVertLevels). I don't see where nVertLevels would be expected in ds_transect and indeed it's also not in mine.

[cbegeman]

@xylar Thanks! I forgot that I had vertVelocityTop in there which has dimension nVertLevelsP1. I just changed the dimension and now it works.

I'm a little confused about the vertical lines because I would expect them to change position if they are the boundaries between cells. However, I've tried 3 different line positions so far and I can't get them to change.

[xylar]

I'm a little confused about the vertical lines because I would expect them to change position if they are the boundaries between cells.

These vertical lines are from setting cell_boundary_color='black'? I do have to agree that the bars look pretty unexpected. I don't see why there should be 4 so close together.

However, I've tried 3 different line positions so far and I can't get them to change

You're calling compute_transect() again each time you change the line, I guess? Yes, I would expect the locations of cell boundaries to change each time unless they happen to intersect cells the same way each time (e.g. you're moving by exactly an integer number of cells).

[cbegeman]

I got the line locations to change by moving diagonally across the domain, so I guess the transects did cut across the cells in similar ways (the lines were perpendicular). Sorry for the false alarm!

[xylar]

@cbegeman, thanks. I did a test where I have a transect zigzagging through the ISOMIP+ domain. I do seem to be seeing the right cell boundaries (red) and skipping the internal triangle boundaries (blue dots): So I think the cell boundaries are working as expected at least in my test.

[xylar]

I know what it is. I'm not removing the periodic edges from the mesh and they're contaminating the plot (most likely) and the cell boundaries (for sure)

[cbegeman]

@xylar Should I go ahead and review now?

[xylar]

@cbegeman, yes you can try this again. I believe I have handling of periodic boundaries in here correctly now. I added plots to baroclinic channel to test it out: This seems feasible given the coarse resolution:

[xylar]

Still working on this...

[xylar]

Okay, I think I'm happy with this but open to suggestions:

[cbegeman]

@xylar My understanding is that if you want to plot the interface locations that correspond to the same time slice that the field is derived from, you have to include the mesh stream in the output stream so that those variables are available when you pass the output dataset to compute_transect. For OMEGA, my understanding is that mesh variables will not be included in the output stream. Have you already thought about adding a separate input argument to compute_transect in addition to ds_3d_mesh that is ds.layerThickness.isel(Time=tidx)?

[xylar]

I don't really understand yet what Omega thinks of as mesh variables and what are not. Until I understand the boundary, I can't really design around it. For now, I've been assuming anything that is generic to any MPAS mesh can be shared but variables like maxLevelCell, dimensions like nVertLevels and anything else associated with the vertical coordinate can't be separated out.

I'm happy to design the transect for 2 different input files but it's not clear to me that layerThickness is sufficient to accomplish this.

[cbegeman]

I don't really understand yet what Omega thinks of as mesh variables and what are not. Until I understand the boundary, I can't really design around it.

Sounds good. We can certainly wait on this.

[xylar]

@cbegeman, how would you feel about there being a ds_horiz_mesh parameter and a ds_vert_coord parameter, where the latter has to have:

• minLevelCell
• maxLevelCell
• layerThickness.isel(Time=tidx)
• bottomDepth

and can optionally have landIceFraction (or that can be a separate parameter)?

[cbegeman]

@xylar I think that makes sense but maybe you're right that we should just wait for OMEGA streams to be worked out so we don't waste time.

Your comment triggered another question though: Can you explain how you're using landIceFraction in compute_transect? It seems like if this is a weight on an ocean field for the interpolation we actually need the grounded land ice fraction which is landIceFraction - landIceFloatingFraction. But your plots look fine so I must be misunderstanding what's going on.

[xylar]

I'm actually not using landIceFraction right now. It was used in MPAS-Analysis to make the background light blue where there was ice and brown elsewhere but I didn't try to do that here. So we could remove it for now.

[xylar]

@cbegeman, I added more parameters to compute_transect() and I'm pretty happy with it now. It makes a lot more sense now in baroclinic channel's viz step because layerThickness is used at the final time rather than from the initial condition.

One thing to be careful about is that you now have to subtract 1 from minLevelCell and maxLevelCell as you pass them in.

I removed landIceFraction for now. We can add it back (along with other changes that may be needed) when we want to plot transects for analysis purposes.

[xylar]

One weird thing is that, if I plot the transect as originally requested: it obviously looks quite different from what actually gets plotted: So it might be better if folks pass the actual coordinates of the transect from ds_transect rather than the ones that were originally requested.

[rljacob]

Late but re: separate input for the mesh variables, uxarray has been assuming that the data is in one file and the mesh is in another because we think modeling groups won't put time-invariant things like the mesh info in every output file.

[philipwjones]

Some mesh info will be time-dependent - both for wetting/drying and for when we start exercising more of the Lagrangian vertical mesh. For Omega, I suspect that will mean that much of the horz mesh can sit in separate static mesh file, but that some masks and much of the vertical mesh is likely to end up with other time-dependent prognostic/diagnostic vars (in terms of output). We'll know more once we start putting that design doc together.

[xylar]

@philipwjones, that's helpful and follows what I was thinking. To support diverse vertical coordinates and processes like glacial isostatic adjustment, bottomDepth, minLevelCell and maxLevelCell (or their Omega equivalents) will all need to be allowed to depend on time.

The solution I have here only assumes (for purposes of transects) that the horizontal mesh is in its own dataset so we should be good.

[xylar]

I ran the baroclinic channel default test one more time and results were as expected.

[xylar]

Thanks for an excellent review, @cbegeman. You had a ton of useful suggestions that made the final result much better.