[PROJ] ISEA projection: proposed parameters deprecation and input latitude interpretation

[Jérôme St-Louis]

Dear PROJ community,

I recently had the chance to integrate an implementation of the inverse 
transformation for the Icosahedral Snyder Equal-Area (ISEA) projection 
<https://proj.org/en/9.4/operations/projections/isea.html> (PR #4211 
<https://github.com/OSGeo/PROJ/pull/4211> which was recently merged).

Doing so resulted in some initial re-factoring to share code between the 
forward and inverse projection.

*Deprecating grid parameters*

During this process, I realized that several parameters for the 
projection (+proj=isea) really have nothing to do with the projection 
itself, but are parameters for constructing discrete global grids on top 
of an ISEA projection.

In fact, the code applying these parameters 
<https://github.com/OSGeo/PROJ/blob/master/src/projections/isea.cpp#L931> 
happens as a later step after the default planar projection (with the 
exception of the trivial isea_tri_plane() 
<https://github.com/OSGeo/PROJ/blob/master/src/projections/isea.cpp#L656>).

Therefore, I suggest this code really does not belong in the PROJ 
library, but outside, in the implementation of a discrete global grid 
system, such as in libraries like DGGRID 
<https://github.com/sahrk/DGGRID> and geogrid 
<https://github.com/mocnik-science/geogrid/>, and these parameters 
should be deprecated.

In addition:

- there is no documentation at all for the parameters and they are very 
poorly understood,
- there are no existing tests for these parameters,
- the new inverse transformation does not support these parameters.

The three parameters proposed for deprecation are:

*+mode* which can be set to /plane/ (the default), /*di*/, */dd/*, or 
/*hex*/
        The default mode is /plane/ which would become the only thing 
supported.

/hex/ referred to returning coordinates related to a hexagonal grid, but 
it is really not clear.
/di/ and /dd/ seemed to refer to a rhombus grid, in a way very similar 
to the approach we use for an ISEA9R DGGRS 
<https://docs.ogc.org/DRAFTS/21-038.html#isea9r-dggrs> (or ISEA4R for a 
quad-tree).
         The ISEA planar projection can be very easily transformed to a 
space where icosahedron triangles merge into 10 axis-aligned rhombuses 
by applying a 60 degrees rotation and 30 degrees horizontal shearing, 
which can be done with +proj=affine.
         This seems to be what isea_ptdd() 
<https://github.com/OSGeo/PROJ/blob/master/src/projections/isea.cpp#L669> 
is doing for /dd/.
/di/ seems to be referring to an extra step done in isea_ptdi() 
<https://github.com/OSGeo/PROJ/blob/master/src/projections/isea.cpp#L849>. 
I have no idea what this is doing.

*+aperture* which is the aperture (refinement ratio) of the discrete 
global grid for those di, dd and hex modes

*+resolution* which I am not sure how it affects the di, dd, and hex 
modes, but is only used with those modes proposed for deprecation

Systems using these parameters could still extract the code for these 
extra steps in its current state 
<https://github.com/OSGeo/PROJ/blob/2ce6c7a9b1ca44b8a7cfb1acca7ff10c8c616774/src/projections/isea.cpp>, 
and invoke as an extra step after the planar projection (after undoing 
the trivial 180 deg rotation and translation that isea_tri_plane() 
<https://github.com/OSGeo/PROJ/blob/master/src/projections/isea.cpp#L656> 
is doing).

If you make use of these parameters and object to their deprecation, 
please reply with feedback.
I plan to submit a pull request for a staged deprecation, which first 
mentions in the documentation that these parameters will be deprecated 
in a future version.

On the flip side, I hope to eventually contribute support for the 
remaining *+azi*= parameter not yet supported by the inverse projection 
(which allows support for orientations such as the R. Buckminster 
Fuller's Dymaxion Orientation).
Possibly also eventually allowing to specify an arbitrary latitude and 
longitude for the first icosahedron vertex, as opposed to the two 
available options of +orient=pole and +orient=isea with fixed values.

*
Proposing to interpret input latitude as authalic latitude*

As discussed in the PR 
<https://github.com/OSGeo/PROJ/pull/4211#issuecomment-2275206640>, the 
ISEA projection is only defined on a sphere.

Although I clarified this in the documentation 
<https://osgeo-proj--4211.org.readthedocs.build/en/4211/operations/projections/isea.html>, 
the most common use case is having geodetic WGS84 latitude and longitude 
coordinates which you want to project to the ISEA planar projection.

This includes the use case of integrating data into a discrete global 
grid system using a DGGRS based on the ISEA projection, such as ISEA3H 
or ISEA9R.

Before applying the spherical projection with a sphere whose surface 
area is the same as the GRS80 or WGS84 ellipsoid (authalic sphere), the 
geodetic latitude should be converted to an authalic latitude 
<https://en.wikipedia.org/wiki/Latitude#Authalic_latitude>.

This is already done by PROJ in other projections, including the Lamber 
Azimuthal Equal-Area projection of which ISEA is modified form, if I 
understand correctly (see 
https://github.com/OSGeo/PROJ/blob/master/src/projections/laea.cpp#L38 
and https://github.com/OSGeo/PROJ/pull/4211#issuecomment-2284383509 ).

I propose that the geodetic latitude be automatically converted to an 
authalic latitude when using the projection with a non-spherical 
ellipsoid (which is the default, and therefore currently quite prone to 
misuse).

Please reply if you see any issue with this change, for which I also 
plan to submit a pull request.

Thank you very much!

Kind regards,

-Jerome
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240824/ac3eadfc/attachment.htm>