Added support for global meshes (stereographic projection)

[tomsail]

Changes were made in the following routines:

• oceanmesh/init.py: added the new functions in the path
• oceanmesh/region.py: added 4 function to pass from lat/lon to stereographic coordinates and vice versa
• oceanmesh/edgefx.py: added supplementary gradient computation to eliminate warped and distorted values of the mesh size around the north (in stereographic projection)
• oceanmesh/geodata.py: rectified the _AREAMIN test (to suppress islands smaller than the resolution) for stereographic projection.
• oceanmesh/grid.py: added stereo keyword in the arguments for the plot() function
• oceanmesh/mesh_generator.py: added the new function _stereo_distortion and the new keyword stereo to the options (opts) to rectify the following functions :
  • _generate_initial_points(): initial distribution in the (-180.00, 180.00, -89.00, 90.00) bounding bbox and reproject them in stereo coordinates
  • _compute_forces(): rectify the computation of the lengths of the triangles segments, either due to the stereographic projection distortion or the elements accrooss the dateline (-180/+180)
• tests/test_global_stereo.py: a test function that produces the following global mesh (0.5° resolution at the coast ~50km):
  

NB: bounding box for global mesh has to be (-180.00, 180.00, -89.00, 90.00)
NBB: it is necessary to define all the coastlines at once: the Shoreline class will the detect the biggest coastline (antartica and define it as the outside boundary). This is because the Shoreline hasn't been defined has a multi-shapefile reader - which would be a "nice to have"

[tomsail]

I can also push the ocean.shp vector I used for the coastlines, as it is really light (42kB). Let me know

[pmav99]

The order of the comments has changed here.

[tomsail]

true, I guess I could have just added the condition at the end and overwrite bSGP.
Let me know what you think, I'll rectify it

[tomsail]

For some reason (I'll need to justify it), after some sensibility tests (resolution, distance, feature and/or wavelength mesh sizes, this last version of the stereo distortion function gives the correct mesh on the north pole

[krober10nd]

@tomsail This is looking quite good. Could it be due to slight variations in the initial point distribution?

[tomsail]

@tomsail This is looking quite good. Could it be due to slight variations in the initial point distribution?

It could be. After the initial distribution - a grid covering the world with a resolution of min_edge_length - the excess points are randomly trimmed out using the mesh size.
I kept this functionality, but also applied the stereographic distortion to adapt the selection in the transformed space.

[tomsail]

@krober10nd I also have missed the stitching of the polygons on either side of the International Date Line (IDL): you can see clearly on the last picture I posted above the Russian peninsula cut from the rest of the continent.

Summing up here is the checklist to complete this PR:

• validate the stereographic distortion
• evaluate the possibility to put the stereographic pole on land (probably Greenland) to avoid the extra gradient smoothing
• add minimal implementation for the stitching of world polygons (see pyposeidon function global_tag)

[tomsail]

Hi @krober10nd, I corrected and simplified the compute_force() routine to account for the stereographic distortion.

I am quite happy now with the quality of the mesh

Except maybe on the poles:

This validate the first item of the checklist above.

The other 2 items on the other hand would require both quite extensive implementations, which are already available in pyposeidon.

2 options :

1. We could merge it for now and address the other items in a further PR,
2. or wait for at least the 3rd item to merge this PR

[krober10nd]

I will take a look at this hopefully this weekend thanks for your great work!!

[krober10nd]

I've seen that convergence pattern in the mesh shape before and it has to do with the forces being too powerful. Could you try re-making the global meshes using a lower pseudo delta tilmestep smaller value? Also, the ratio of the medians for the scale factor could work better (in general) but that's not related to your changes.

For example https://github.com/CHLNDDEV/OceanMesh2D/blob/f0c9257bc8f1a00836510e251328f07ef288324b/%40meshgen/meshgen.m#L811

I'm on board with merging this in, just please address the flake8 aspects and we're good to go.

[krober10nd]

Seems like scipy's interpolation routine now requires the coordinates to be montonically increasing causing most of the tests to fail.

[tomsail]

Awesome Keith thanks.
I'll fix the lint problems, do some sensitivity with delta_t around the north pole, and try to propose also a workaround for the error, caused by create_vectors():

    y = self.x0y0[1] + np.arange(0, self.ny) * abs(self.dy)
    y = y[::-1]  # descending monotonically

[tomsail]

I implemented 2 fixes and one enhancement, see the last commits:

1. FIX: Monotonically ascending vectors.
2. FIX: For the CIs, fixing the monotonically ascending vectors wasn't enough, there were outdated arguments in the grid.plot() function, I basically copied and adapted what was already there in the DEM.plot() function ). This should fix some of the tests. I tested everything locally only for python3.10 and commented the parts with the bathy gradient (because I knew it wasn't working already).
3. I also implemented the median ratios (it was already implemented in the commented version of compute_forces() )

[tomsail]

sorry I realised I didn't check the last implementations I did with flake8 ..

[krober10nd]

That's not a problem, I can fix that. Did the modifications to the forces improve the mesh quality near the poles (or the experiments with reducing the timestep)?

[tomsail]

No I couldn't enhance or fix it

[tomsail]

All CIs should pass now (I fixed flake8 and lowpass filter for bathy gradient)

[krober10nd]

I will remove the bathymetric gradient test as this function is not functional and needs to be revisited regardless.

[tomsail]

FIY now the bathy gradient function works with the default lowpass filter.
Reason for the function hanging was due to the ESPG problem (conversion degrees<>meters wasn't done at all) #69

[tomsail]

@krober10nd, I will

1. clean all the work in this PR (git rebase),
2. check if the CIs work on my fork (I hadn't activated it before, now I do)
3. open it again for review once all CIs pass, so you can have a last look before merging

Note:
I identified 2 issues that I will address in 2 different PRs for more clarity:

1. The ESPG argument (mentioned above)
2. A bug in the bathy gradient through the grid.values setter

[krober10nd]

@tomsail everything looks good. could you add an example on the README with an image to highlight this new feature?

[tomsail]

I added supporting files for global coastlines. It'll be clearer for users.
Also this partly covers the following TODO:

• add minimal implementation for the stitching of world polygons (see pyposeidon function global_tag)

as I refer to the global_tag() function in pyposeidon. Either in test_global_stereo() or in the README

I think this is a wrap