stereographic distortion close to the poles

[tomsail]

This issue is related to what was mentioned in the PR for global meshes about the stereographic distortion.

Close to the poles, the triangles get distorted and don't have the right equilateralness.
As mentioned,

it has to do with the forces being too powerful

To solve this issue, I have implemented a parametric law in the mesh_generator():

def _parametric(lat):
    ones = np.ones(lat.shape)
    res = ((90 - lat) * 2 + 18) / 180 * np.pi
    return np.minimum(res, ones)

and illustrated it in the following python notebook

I will propose a PR with this parametric law. But I am still not convinced this is the best way to go..

The meshes look almost perfect closer to the poles now:

Keen to get another opinion on this.

[WPringle]

I think if the results are good (equilateral and matching the edge length that you want) then it's reasonable approach. What happens for a non-global mesh around the poles, like for Arctic circle only?

[tomsail]

Thanks @WPringle, following your post, I realized it didn't work for other resolutions.

So I have done some sensitivity tests on the arctic circle and ended up using the following parametric law:

def _parametric(lat):
    ones = np.ones(lat.shape)
    res1 = ((90 - lat) * 2 + 16) / 180 * np.pi
    res2 = ((90 - lat) * 4 + 8) / 180 * np.pi
    res3 = ((90 - lat) * 8 + 4) / 180 * np.pi
    res4 = ((90 - lat) * 16 + 2) / 180 * np.pi
    res5 = ((90 - lat) * 32 + 1) / 180 * np.pi
    res6 = ((90 - lat) * 64 + 0.5) / 180 * np.pi
    res7 = ((90 - lat) * 128 + 0.25) / 180 * np.pi
    res8 = ((90 - lat) * 252 + 0.125) / 180 * np.pi
    res = np.minimum(res1, ones)
    res = np.minimum(res2, res)
    res = np.minimum(res3, res)
    res = np.minimum(res4, res)
    res = np.minimum(res5, res)
    res = np.minimum(res6, res)
    res = np.minimum(res7, res)
    res = np.minimum(res8, res)
    return res

I have implemented it and also changed the initial point distribution during the mesh generation.
See the last changes

I have tested up until 3km, with a uniform mesh size or a varying one.

Equilateralness is okayish but I am not sure if I get the right edge length.
I still end up with a singularity around the north pole (or (0,0) in stereo projection).

I guess for now, I will avoid constraining my model with too high resolution at the pole.

[tomsail]

I found out what was the problem.
During the re-computation of the forces, we were over-constraining because there were too many points closer to the poles.

I was actually trying to solve this by the wrong end:

• I was trying to adapt the forces in _compute_forces() although its definition was correct.
• the problem was coming in fact from the _generate_initial_points() function.

In generate_mesh, here is how we generate the points before the mesh iteration:

if stereo:
    bbox = np.array([[-180, 180], [-89, 89]])
p = np.mgrid[tuple(slice(min, max + min_edge_length, min_edge_length) for min, max in bbox)].astype(float)

this distribution is fine if you want to mesh in the WGS84 projection.
However, since we want to generate global meshes, we need to reduce the number of point along the latitude.

Same as for regular gaussian grids, the longitude rows need to decrease towards the poles.

I compared here the density of longitudes along the latitude for the following ECMWF's Reduced Gaussian grids:

• N640
• N512
• N400

Or, in other terms, above is the evolution of reduced_points/regular_points along latitude.
It corresponds to: np.cos

In the same function (source), we define:

r0m = np.min(r0[r0 >= min_edge_length])
p = p[np.random.rand(p.shape[0]) < r0m**2 / r0**2]

to trim out extra points points depending on the mesh size.
Since we need to address that min_edge_length changes along the latitude, I have defined:

def _edge_length_wgs84(lat):
    lrad = np.radians(lat)
    return 1 / np.cos(lrad) ** (1 / 2)

but instead of changing the whole equation of r0m, I implemented the change in r0.

[tomsail]

The sensitivity test on the polar circle are spot on now:
here is the latest version of the notebook

[WPringle]

@tomsail That makes sense, lon spacing definitely should be adjusted according to latitude irrespective of global mesh or not I reckon. If edge length is specified in degrees then we should assume it is meant degrees in latitude or "equator degrees".

[tomsail]

Good point. I could provide an example for Iceland, Greenland and see the differences with/without that correction.
I'll do it later though. I have no rush to merge this right away + I can use my fork for my use cases now + I'd finally also wait for the CI to be fixed and merged (see #78)

[tomsail]

Hi @WPringle and @krober10nd,

It's been a while, but I had some time to look at this.
I have tried to implement the distorsion and see how it impacts the meshing close to the polar circle.

here are the results for area above Greenland.
with

and without

I don't think that using lon/lat projection is the best method to do this.

I have documented this in a notebook where I tested 3 different methods:

1. normal meshing using current config with EPSG 4326 (WGS84)
2. normal meshing using current config with EPSG 3995 (Arctic Polar Stereographic)
3. meshing using distortion (i.e. the _edge_length_wgs84(), both in _generate_initial_points and _compute_forces

I think method 2 is the safest method, because no assumption is made under the radar of the user. A feature to transform DEMs from WGS84 to any other projection would be nice but then again, it can be done separately too.

[tomsail]

In conclusion:

• the PR added parametric law for stereo distortion (north pole) #76 addresses this issue only for global meshes.
• there would be more investigation needed for distortion and meshing close to polar locations.
• I think we can merge added parametric law for stereo distortion (north pole) #76, but keep this issue open, since we haven't completely solved it