Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added support for global meshes (stereographic projection) #68

Merged
merged 1 commit into from
Dec 4, 2023

Conversation

tomsail
Copy link
Contributor

@tomsail tomsail commented Oct 9, 2023

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):
    image

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
Copy link
Contributor Author

tomsail commented Oct 9, 2023

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

bSGP = bSGP.difference(pSGP)
# Append polygon segment to mainland
mainland = np.vstack((mainland, poly))
# Clip polygon segment from boubox and regenerate path
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The order of the comments has changed here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Contributor Author

tomsail commented Oct 9, 2023

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
image

@krober10nd
Copy link
Collaborator

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

@tomsail
Copy link
Contributor Author

tomsail commented Oct 10, 2023

@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
Copy link
Contributor Author

tomsail commented Oct 10, 2023

@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
Copy link
Contributor Author

tomsail commented Nov 8, 2023

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
Screenshot from 2023-11-08 17-58-45

Except maybe on the poles:
image

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
Copy link
Collaborator

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

@krober10nd
Copy link
Collaborator

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
Copy link
Collaborator

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

@tomsail
Copy link
Contributor Author

tomsail commented Nov 13, 2023

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
Copy link
Contributor Author

tomsail commented Nov 13, 2023

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 tomsail marked this pull request as ready for review November 13, 2023 15:56
@tomsail
Copy link
Contributor Author

tomsail commented Nov 13, 2023

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

@krober10nd
Copy link
Collaborator

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
Copy link
Contributor Author

tomsail commented Nov 14, 2023

No I couldn't enhance or fix it

@tomsail
Copy link
Contributor Author

tomsail commented Nov 14, 2023

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

@krober10nd
Copy link
Collaborator

krober10nd commented Nov 14, 2023

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

@tomsail
Copy link
Contributor Author

tomsail commented Nov 21, 2023

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 tomsail marked this pull request as draft November 21, 2023 14:37
@tomsail
Copy link
Contributor Author

tomsail commented Nov 21, 2023

@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

@tomsail tomsail marked this pull request as ready for review November 22, 2023 11:29
@krober10nd
Copy link
Collaborator

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

@tomsail
Copy link
Contributor Author

tomsail commented Dec 4, 2023

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

@krober10nd krober10nd merged commit c4c56d2 into CHLNDDEV:master Dec 4, 2023
6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants