Create southern boundary features for each MOC basin

[xylar]

With this merge, setup_ocean_basins first removes small polygons
from MOC basins that are not useful for MOC computations. Then,
it finds the southern boundary of each MOC basin. Finally,
the southern boundaries are subdivided into 0.1 degree segments
to ensure that the MPAS mask creator follows them with high
fidelity.

[xylar]

This PR is a replacement for #75, where the approach I used didn't work out.

[xylar]

@milenaveneziani, I implemented the features I described here: #75 (comment)

I was able to create the masks with the mask creator using:

./MpasMaskCreator.x mesh.nc MOCSouthernBoundaries.nc -f MOCSouthernBoundaries.geojson

mesh points to the QU_240km mesh. Then, I visualized it in paraview after using the following paraview extractor command:

./paraview_vtk_field_extractor/paraview_vtk_field_extractor.py -m mesh.nc \
-f MOCSouthernBoundaries.nc -v all \
-d maxEdges= nTransects=: TWO=: vertexDegree= maxEdges2=

[xylar]

@milenaveneziani, the results described above look reasonable to me. If we don't yet have a way to test them as part of an MOC calculation, it might be best to go ahead and merge this PR, and we can always make corrections in the future. Let me know what you think of this approach overall and if I can clarify what I've done.

[xylar]

I just pushed some changes to make the code PEP8 compliant (at least those aspects that flake8 checks). I also updated the docstring for setup_ocean_basins.py. @pwolfram, you may want to have a look at the code to make sure you're good with the style.

[pwolfram]

@xylar, I took a quick look and didn't notice anything pressing. Eventually we will want to standardize style, e.g., #54, but this probably should not be a priority right now.

[milenaveneziani]

@xylar: thanks! This makes sense to me.
I will test it tomorrow using the 60to30 mesh.
@DieMuhkuh can also test it with the MOC AM, once we have the maskfile.

[milenaveneziani]

@xylar, a side note: I made the mask file using the 60to30 mesh, but when I ran the paraview extractor using your command logic above, I get the following error:
Traceback (most recent call last):
File "./paraview_vtk_field_extractor/paraview_vtk_field_extractor.py", line 779, in
valid_mask, use_32bit, args.combine_output, args.append )
File "./paraview_vtk_field_extractor/paraview_vtk_field_extractor.py", line 534, in build_field_time_series
assert('Time' not in mesh_file.variables[var_name].dimensions)
AssertionError

which in a sense takes me back to the issue of what exactly we want/need to have in a mesh file in order for our tools/analysis to work always.

Anyway, @DieMuhkuh, @mark-petersen : I now have the mask files (basins and transects) for the 60to30 mesh. Do you think you can test the MOC AM with these when IC comes back up?

[xylar]

@milenaveneziani, I think you are probably working with an out-of-date version of the paraview extractor. This issue has been fixed.

which in a sense takes me back to the issue of what exactly we want/need to have in a mesh file in order for our tools/analysis to work always.

You don't really have a choice about what's in the mesh file. It's a file that defines the MPAS mesh, so has nothing specific to MPAS-O in it. It comes from the mesh database as part of setting up a testcase.

The mesh file is not the same thing as an output file, though the output file may contain all the information form the mesh by adding the the mesh to its stream definition. The discussion is that the mesh information is rather large and doesn't need to be in every output file, since we always know that its contents are in the mesh file (usually mesh.nc). The latest version of the paraview extractor should work when you supply it with a mesh file separately from the the time-dependent output file.

[milenaveneziani]

hmm, I was using the version that @pwolfram was working on for facilitating the building of MPAS-Tools on different platforms (because I was using edison, and therefore referred to the comments in MPAS-Dev/MPAS-Tools#136).

[xylar]

@milenaveneziani, I see. That branch is a good choice for the mask creator, etc. but not for the paraview extractor. But it does indicate that we should push to have MPAS-Dev/MPAS-Tools#136 merged sooner rather than later so both can be found in the same branch.

[milenaveneziani]

@xylar, @DieMuhkuh: @vanroekel used the MOCSouthernBoundary mask file that I created from the geojson file of this PR and the 60to30 mesh, to compute offline MOC (from a script that @mark-petersen provided). He finds weird results (Atlantic MOC has a weird structure, with very small positive values in the top 1500 m and large negative values underneath). He then used the same MOCbasin mask, but a transect that Mark used in the past, and the results look much more reasonable. So the problem seems to be the southern boundary mask (MOCBasin mask is fine).
I will try to visualize it with matlab or the latest paraview extractor.

[xylar]

@milenaveneziani, @mark-petersen and @DieMuhkuh, sorry this isn't working as you need it to. It could be that I haven't understood correctly what the requirements of the southern boundary feature, in which case it would be good to have this clarified as part of an issue or a confluence page. It could also be that geometric_features is not the right place to create the feature you need to define the southern boundary because specific information about the MPAS grid is needed for this purpose.

My understanding is that what is needed is a transect that includes all edges on the MPAS grid that make up the southern boundary of each MOC region. Again, if I have understood correctly, this transect would include all edges with one neighboring cell in the MOC region and one cell out of it. Presumably for convenience, the transect could also include edges of cells along the southern boundary that border land but I'm not sure if that's a requirement.

To help you work on this, I would need a better sense of what @mark-petersen's southern boundary transects look like and how these differ from what I provided. While I understand that the resulting MOC was not satisfactory, this information isn't very helpful to me in determining what the southern boundary feature should look like to fix the problem. In short, I would not like to continue to operate in a way where I produce a feature, wait several days, and find out that that's not what you were after, can I try again?

[milenaveneziani]

@xylar, @mark-petersen, @DieMuhkuh: would it help if we visualized the MOC southern boundary feature that resulted from this PR and also the corresponding mask file? I meant to do this yesterday, but didn't make.
At least Mark and Nils can see where the transect is flawed and what we should do (either in geometric_features or in the mask creator) to fix it.

[xylar]

@milenaveneziani, ideally we would visualize the MOC boundaries produced by @mark-petersen's script and the one produced by this PR so I could understand better how they differ. It would also be very helpful to understand how the MOC southern boundary produced by @mark-petersen's script corresponds to the MOC basin feature (i.e. is it exactly the edges that make up the southern boundary as I described above; are those edges complete or are there gaps where the southern boundary intersects land?).

Whatever you can do to help clarify what is needed would be helpful to me.

[milenaveneziani]

ok, so this is the MOC southern boundary feature:
https://gist.github.com/milenaveneziani/aeade438efdd417f9c896d3edf6e89cf

The problem may be in the Australia and Papa New Guinea regions.

[xylar]

@milenaveneziani, could we try to get a more specific requirement? Would it be sufficient if we had a single line that went straight from each westernmost point to each easternmost point without following the land boundary?

[milenaveneziani]

I think that that should be sufficient, but again @mark-petersen and @DieMuhkuh should confirm this.
I have the mask file that Mark used in the past for the offline MOC calculation: I hope to visualize that as well in a moment.

[xylar]

I spoke with @milenaveneziani on the phone and we agreed that the Atlantic MOC should have worked if it were sufficient to have a transect that's a single straight line for each MOC basin, whereas it didn't work. This suggests that maybe the southern boundary can't be found in geometric_features and may need to be extracted after masks have been created by the mask creator. again, we'd benefit from thoughts from @mark-petersen and @DieMuhkuh

[DieMuhkuh]

@xylar As far as I'm concerned, the transect would be fine as long as each edge only appears once and with the correct sign. Since we only consider the transport through each facet below a transect edge, it should not make a difference if the transect goes around land, since the water transport through land should be zero.

[xylar]

Thanks, @DieMuhkuh. @milenaveneziani and @DieMuhkuh (and maybe @mark-petersen), can you look into why the southern boundary for the Atlantic MOC basin didn't work? That one should have no land (or weird double-counting of edges) involved and yet it sounds like it didn't work. Is that because the mask creator just got the wrong edges?

[milenaveneziani]

@xylar: I may have a hint of why the Atlantic MOC boundary didn't work, but I have to verify it. I can try tomorrow (sunday), although I may have to postpone to Tue (Mon is already full for other reasons). So, let us put this on hold for a couple of days. thanks for the huge help so far.

[milenaveneziani]

@xylar and all: I did not have much luck with finding out what the problem may be with the Atlantic_MOC southern boundary generated here. Here are some plots though, to see if other people (especially @DieMuhkuh and @mark-petersen) can make an informed guess.
First, this is a picture showing the mask for the southern Atlantic_MOC boundary obtained by running the mask generator on the corresponding feature computed from the code of this PR (and using the 60to30 grid):

You can see that there are various places where two cells share an edge in the meridional direction: is this going to be an issue?

And this is the southern boundary used by Mark for previous MOC (offline) computations:

I looked closely at the above mask and there are no edges to appear twice (all cells only share one edge with the others).

[xylar]

@milenaveneziani, it sounds like maybe the behavior of the mask creator with a LineString made up of multiple segments isn't going to work. With a single line segment, the mask creator finds a path from the closest edge to the first point to the closest edge to the last point, never reusing an edge, which is the desired behavior. The undesired behavior is that it does so in a zig-zag pattern, as above.

@milenaveneziani, @mark-petersen and @DieMuhkuh, am I right in thinking this mean we need to fix the mask creator (MPAS-Dev/MPAS-Tools#143) and that if we did so, a simple segment (with no subdivision) across the bottom of each MOC basin would be the desired geometric feature?

[milenaveneziani]

And these are Atlantic_MOC plots.
This first one is obtained by using Mark's Atlantic_MOC regionMask and Mark's Atlantic_MOC southern boundary (ignore title weirdness):

This second one is done by using our new Atlantic_MOC regionMask, while keeping Mark's mask for the southern transect (note that the plot is very similar to the above plot, except that it extends further to the north):

And finally, here is the same plot obtained when using both our new Atlantic_MOC regionMask and the new southern boundary:

[milenaveneziani]

@xylar: may I venture a question? Since we have the MOCBasin mask already, could we simply pick the cells that make up the southern border of that? We could then be sure that that southern boundary is really the southern boundary corresponding to that particular MOCBasin mask, as it should be.

[DieMuhkuh]

@milenaveneziani In my opinion, this would be the most correct solution to this problem. The problem I see with that is, that the mask creator would have to know which regions are MOC regions so that it knows when a southern boundary has to be calculated.

[milenaveneziani]

@DieMuhkuh: we do have the regionGroup in the features. The mask creator could be instructed by that.

[xylar]

@milenaveneziani and @DieMuhkuh, I think this requirement is too specialized to try to build it into the mask creator. However, it would probably be quite easy to build a specialized python tool to extract the southern boundary of the MOC regions from the NetCDF file (assuming we could give a satisfactory definition of what is meant by the southern boundary). The python tool could create the same transect fields that the mask creator would make, so the workflow would be essentially unchanged once the tool was run.

Why don't I try to give this a go, and we can see if the results are satisfactory?

[milenaveneziani]

ah! that seems like a great option, @xylar. So we don't have to change the mask_creator at all.
Do you think you have all the elements to give this a try, or are there other clarifications about how the MOC am works that should be made?

[xylar]

@milenaveneziani, I think I have a pretty good idea of what it needed. If the results of my tool aren't quite what's needed, at least we'll have something to tweak. But I'm pretty confident I can get you what you want.

[xylar]

I am closing this PR but keeping the useful parts in #79. As discussed above, we will pursue alternative methods for finding the southern boundaries of MOC features.