IGEO7 ellipsoidal adjustments as optional in/out-processing option

[allixender]

with the latest improvements from #18 we can make IGEO7 compatible with IVEA7H ISEA7H (+Z7).

We use the vertex repositioning and add geographiclib as dependency and provide a post-processing option to make authalic adjustment:

LandscapeGeoinformatics/pydggsapi#11 (comment)

[allixender]

We should also provide additional more explicit documentation on the DGGRIDv8 parameter use and have a separate docs on the IGEO7 adjustments towards OGC DGGS API compliance/compatibility , like rotation and authalic adjustment:

• https://github.com/allixender/dggrid4py/blob/master/docs/source/IGEO7.rst
• https://github.com/allixender/dggrid4py/blob/master/docs/source/dggrid4py/dggrid4py.igeo7.rst
• https://github.com/allixender/dggrid4py/blob/master/docs/source/dggrid4py/dggrid4py.DGGRIDv8.rst

@fmigneault @tik65536

[jerstlouis]

@allixender Did you mean ISEA7H (+Z7)?

As far as I know for IVEA you would need to tweak the icosahedral projection inside DGGRID (which would be nice as well).

[allixender]

Well, for reasonably small enough cells, I thought that IVEA7H just has smoother boundaries, but not completely different DGGH?
But, yes, primarily ISEA7H. (edited above)

[allixender]

However, geographiclib Python package does not expose AuxLatitude functions 😠 https://geographiclib.sourceforge.io/C++/doc/auxlat.html

That's a bit of a bummer, we need to either expose them example, or find an alternative like: pygeodesy

pygeodesy does like so:

>>> import pygeodesy
>>> import pygeodesy.auxilats
>>> from pygeodesy.ellipsoids import Ellipsoids
>>> wgs84 = Ellipsoids.WGS84
>>> test_el_lat = 52.12345678
>>> to_lat_geodetic = lambda lat_authalic: wgs84.auxAuthalic(lat_authalic, inverse=True)
>>> to_lat_authalic = lambda lat_geodetic: wgs84.auxAuthalic(lat_geodetic, inverse=False)
>>> to_lat_authalic(test_el_lat)
51.99904725957745
>>> to_lat_geodetic(to_lat_authalic(test_el_lat))
52.12345678

[jerstlouis]

@allixender

Well, for reasonably small enough cells, I thought that IVEA7H just has smoother boundaries, but not completely different DGGH?

IVEA7H is a completely different DGGH. Even for relatively small zones, the zones end up at a completely different place.

e.g.,

https://maps.gnosis.earth/ogcapi/collections/HRDEM-RedRiver/dggs/ISEA7H_Z7/zones/021366213

https://maps.gnosis.earth/ogcapi/collections/HRDEM-RedRiver/dggs/IVEA7H_Z7/zones/021366213

[tik65536]

@allixender , @fmigneault

I want to clarify one thing. The change in orientation to shift 0.05 degrees to the West is due to the authalic latitude?
My understanding is that when converting the original latitude and longitude of vect_0 to ellipsoidal coordinates using authalic conversion, it falls within the land of Sweden, so it is not "truly" on the sea.

When we manually adjusted the vert_0 by the setting dggs_vert0_lon and dggs_vert0_lat, are we still using the coordinate based on the authalic sphere?

I think we should have an authalic_convert mode setting when we instantiate dggrid4py (calling DGGRIDv7() or DGGRIDv8() ), so that everything involving coordinates should perform the authalic conversion (including the settings) according to the mode is on or off. It will be more consistent.

[jerstlouis]

@tik65536

The exact authalic latitude to achieve an orientation symmetrical about the equator is arctan(golden ratio) = arctan((1 + sqrt(5)) / 2) = 58.28252558853899... (authalic latitude).

When converting that to geodetic latitude, that yields ~58.397145907431 geodetic latitude.

DGGRID essentially always worked with authalic latitude, so the coordinates of dggs_vert0_lat is the authalic latitude (58.28252558853899...).

The geodetic coordinates 58.397145907431 N, 11.25 E is in Sweden.

If we shift it 0.05 degrees west to 58.397145907431 N, 11.20 E, we clear the land (starting from refinement level 9), without any other vertices ending up on land.

Note that the authalic latitude conversion itself does not affect the longitude coordinates.

EDIT: To clarify, the authalic latitude should always be 58.28252558853899..., whether the conversion is performed or not, to maintain the orientation symmetrical about the equator.

So if you want to specify the location of the first vertex in geodetic latitude when performing the conversion, the latitude would need to be specified as 58.397145907431... N.

[fmigneault]

I am planing to regenerate the IGEO7 and H3 manitoba data with the vert0 correction.
If this snippet (#28 (comment)) is the one we decide to integrate, I'll include it as well in that regeneration.

We should also provide additional more explicit documentation on the DGGRIDv8 parameter

I agree about this. Even though the options are now available via dggrid4py, they still have to be provided explicitly to override DGGRID's defaults.

I think we should have an authalic_convert mode

Agreed. It should be done automatically if enabled (and probably be done by default). Users would have to explicitly disable it if they wish to preserve the older behavior producing authalic latitudes. I think most would expect geodetic latitudes to be produced by default for seamless integration is other GIS tools.

[allixender]

@fmigneault @tik65536 please have a look at https://dggrid4py.readthedocs.io/en/latest/dggrid4py/dggrid4py.auxlat.html

There are some draft functions available for authalic conversion. It is really not very fast, several seconds for ca 20000 polygons/hexagons (with 6 vertices, 120000 calculations). Some numpy vectorization or similar will likely be needed to make it faster.

@tik65536 can you review your legacy driver test when you have time, the clip coarse cell example failed (tests/test_legacy_driver_name.py::test_grid_cell_polygons_from_cellids_coarse_cells)