MOM6-community / sectionate

compute section in ocean models
GNU General Public License v3.0
8 stars 7 forks source link

Robust handling of corner cases (e.g. section crossing ocean_grid boundary) #13

Closed hdrake closed 4 months ago

hdrake commented 1 year ago

It seems like we've just gotten lucky that until now our regions of interest have not wrapped around the periodic zonal coordinate of ocean grids. This is not generally going to be the case, so we need to be able to deal with this robustly.

I see two possible options: 1) use da.roll({'xh':idx, 'xq':idx}) so that the sections no longer wrap around the zonal axis 2) add conditional logic to infer_broken_line to handle the case where the shortest path crosses the grid's periodic x boundary

I think option 2) will end up being more straight-forward and robust. Probably the easiest thing is to just check if np.abs(i2-i1) >= nx/2, in which case the faster direction will be to wrap around the boundary. The current implementation is only able to stay within the bounds of the geolon coordinate (here [-300, 60], as shown by the dashed black line) and so find the very long route that goes all the way around the world in the other direction (see figure below).

Screenshot 2022-12-16 at 11 28 09 AM (Note: I've actually used option 1 to find the correct orange section in the, but in doing so I realized option 2 is a better approach; I just have to code it up now.)

@raphaeldussin, this is a major bug so my priority is to make a new pull request to fix this without any of my other changes so it can be merged ASAP into the master branch as a hotfix.

raphaeldussin commented 1 year ago

my original thought was that users could roll their array if the have to cross the periodic boundary. However I can see how that might be impractical if you want to extract multiple close sections,... I'll be happy to pull your changes in.

hdrake commented 1 year ago

my original thought was that users could roll their array if the have to cross the periodic boundary. However I can see how that might be impractical if you want to extract multiple close sections,...

Yes, also problematic for zonally-circumpolar sections (like Becki's isobath-based section for Dense Shelf Water export off of Antarctica or circum-Arctic sections) like:

Screenshot 2022-12-16 at 3 02 05 PM

If a user is careful they could probably carefully using roll as a workaround even for the circumpolar case, but definitely very complex sections like the one below would be really difficult:

Screenshot 2022-12-16 at 3 04 10 PM

These work intuitively by default with my new implementation of method (2), which just modifies a few lines to your infer_broken_line function. I will make a Pull Request later today.

hdrake commented 1 year ago

Addressed by https://github.com/raphaeldussin/sectionate/pull/16/commits/ff1ec9cc966fea93b96e3c7ff757215fe687fe0f of PR https://github.com/raphaeldussin/sectionate/pull/16 (superceded by https://github.com/raphaeldussin/sectionate/pull/19); can close when the PR is merged.