Feature: Snyder Polyhedral Equal Area projections

[felixpalmer]

• Claude Code supported my development of this PR. See our policy about AI tool use. Use of AI tools must be indicated.
• Tests added
• Added clear title that can be used to generate release notes
• Fully documented, including updating docs/source/*.rst for new API

Background

Please read the updated documentation in docs/source/operations/projections/polyhedral.rst first as well as Snyder's original paper: AN EQUAL-AREA MAP PROJECTION FOR POLYHEDRAL GLOBES (CARTOGRAPHICA Vol 29 No 1 SPRING 1992 pp 10–21)

Back-compatibility with existing isea

Gives same results, with comparisons in test/gie/builtins.gie. I would propose we remove the old isea code, replacing with this implementation.

Architecture

To define a projection a polyhedron and net must be defined. The simplest example is that of the 4-sided tetrahedron along with the net definitions:

• tetrahedron.h specifies a list of 3D vertices and faces indexing those vertices
• tsea.h specifies a list of parents, with 0 as the root

These details are not exposed to the user via the PROJ API, instead they reference them by name, e.g. +proj=tsea +net=star

When building the projection, the following process followed:

• Options parsed and populated into pj_polyhedral_data struct
• polyhedral::load_meshes called with the selected polyhedron and matching net, which:
  • Unfolds the net by starting with the root face and then attaching other faces in a 2D plane
  • Applies a Conway meta operation to cut up the polyhedron and net faces into right-triangles
  • Stores the 2 resulting triangle lists on pj_polyhedral_data
• The Snyder forward and inverse projection is then defined independently between each of these paired triangles

The design is generic and adding polyhedra/nets is straightforward and doesn't require large code updates

Change list

• Snyder projection functions: snyder_fwd & snyder_inv
• Unfolding helpers
• Conway helpers
• Spherical geometry helpers
• Basic 3D vector operation and helpers
• 3 polyhedrons: Tetrahedron, Dodecahedron & Icosahedron and some nets for each
• Tests

Example renders ISEA

Example renders TSEA

Example renders DSEA

[felixpalmer]

@jerstlouis @rouault @sahrk @jjimenezshaw @ldesousa @mikima @busstoptaktik @ppKrauss @thangqd

[jerstlouis]

@felixpalmer Lots of cool stuff here :) Thank you for this!

As discussed in #4389 (comment) , I was also planning to propose a replacement for the current ISEA projection based on a generalized ISEA / RTEA / IVEA method, using a closed-form inverse solution, and address some issues with the current set of parameters.

You've taken that a step further with support for different polyhedron nets (and multiple unfolding setups too!!) which is also very cool. Is this using Snyder's original equations, which did not have a closed-form inverse, or the Slice and Dice approach with closed form inverse? We should discuss :)

I was also going to test some further improvements suggested by @brsr which could potentially fully remove the need for Slerp in the inverse.

Roundtrip test (less than 1nm roundtrip error)

That is quite impressive. Was that tested throughout the globe surface or only for some specific test points?

cc. @mpadillaruiz @ppKrauss

[felixpalmer]

@jerstlouis the equations are from A5, DGGAL and BRSR, as noted in the header. I didn't reference Slice & Dice for the code, but it is possible that the 3 sources did.

I'm open to improvements to the code and equations (e.g. removing the slerp()), one thing I would like to keep is the modular framework I've designed, which makes it easy to add new polyhedra/projections/nets (see my follow-up PRs for how easy it is)

The roundtrip error is just based on a static set of points, you're right that perhaps it is worse if this is expanded across the globe

[jerstlouis]

@felixpalmer In particular, to replace the current isea, we should:

• Have orient_lat, orient_lon and azi parameters that determine the position of the first vertex (with geodetic latitude), and azi as the azimuth to the next vertex (based on whichever polyhedron is chosen), and for the icosahedron net unfolding (ISEA / IVEA / RTEA) the default should be authalic latitude at arctan(golden ratio) converted to geodetic latitude, 11.20 East, and possibly 11.25 E in the case where you specifically request the sphere rather than GRS80 as the ellipsoid for backward compatibility with the current +proj=isea projection.

EDIT: It seems that you already support this with polyhedral::set_orient_from_angles(Q, ISEA_STD_LAT_DEG, ISEA_STD_LON_DEG, IVEA_AZ_DEG);, though as I was saying we should use 11.2 E instead of 11.25 E when the GRS80 ellipsoid is used (which should be the new default, latitude should be ~58.397145907431 N (geodetic).

Would specifying the rhombic triacontahedron solid naturally yield the RTEA projection (as in DGGAL, https://www.mdpi.com/2220-9964/9/5/315, and https://www.tandfonline.com/doi/full/10.1080/17538947.2022.2130459?utm_source=researchgate.net&utm_medium=article ?

I see the special names of Hexakis icosahedron and Decakis dodecahedron that you give to the disdyakis triacontahedron, but based on Snyder's paper specifically, the center of the faces is what is used for tracing the great circles, so the names of the solids would be the Dodecahedron and the Icosahedron (as per the DSEA and ISEA acronyms), right? And for Rhombic Triacontahedron that also naturally yields RT(S)EA. It's a very nice feature to be able to unfold them onto the different nets. Would this unfolding also work on polyhedrons without icosahedral symmetry? The unfolding net needs to have same symmetry as the projection solid?

[felixpalmer]

@jerstlouis could we continue the ISEA discussion here? I've given you access, feel free to update the constants - you understand them better than me. As you say it is already supported by the code.

naturally yield the RTEA projection?

I believe so, the code is designed to be a generic spherical-planar triangle mapper so you would just need to generate the relevant config files. The ISEA PR is a good example of what is required to add a new projection & polyhedron

special names of Hexakis icosahedron and Decakis dodecahedron

Yes I thought a lot about the naming. I settled on this approach as the name tells us which vertex to use as the apex and how many triangles to cut the face into. See the docs for details

For brevity (and following the example set by ISEA) I chose to just use the underlying platonic solid when naming the projection, e.g. DSEA, TSEA etc

Would this unfolding also work on polyhedrons without icosahedral symmetry?

Yes, this PR literally implements TSEA, which is on a (hexakis) tetrahedron. It should also work on other more complex nets, like the truncated icosahedron or even a Myriahedron. The unfolding net just needs to have the same number of triangles. Technically, they don't even need to be connected as a net, you could float them in space separately. And implementing your 5x6 net is as simple as adding a new config file

[jerstlouis]

Yes I thought a lot about the naming. I settled on this approach as the name tells us which vertex to use as the apex and how many triangles to cut the face into.

Yes, but technically these all refer to the same polyhedron, and I wonder if there couldn't be a more explicitly way to describe it, especially if we consider the possibility to specify the polyhedron and the choice of radiating vertex as projection parameters.

For the RTEA configuration, would the name specifying that great circle setup mapping for the d120 be kisrhombic triacontahedron ?

Also of course I'm very curious whether this modularity would help test the Archmidean solids configuration like the Truncated Icosahedron (TISEA / TIVEA), and still map that to the icosahedral net and the 5x6 space :)

[felixpalmer]

RTEA configuration, would the name specifying

I think it should be tetrakis (rhombic) triacontahedron - so a triacontahedron with each face cut into 4

It is true that these all have the same vertices, but by naming them this way you can avoid the whole "radiating vertex" concept as the apex is always in the center of the base polyhedron face (the tip of the kis pyramid)

[jerstlouis]

@felixpalmer Thanks, I am starting to understand how these *kis prefixes work :)

tetrakis rhombic triacontahedron makes sense, kis pyramids with 4 faces on each of the 30 faces of the rhombic triacontahedron. kisrhombic triacontahedron is one of the name of the d120 mentioned on https://en.wikipedia.org/wiki/Disdyakis_triacontahedron and likely just a shorthand for that longer name. And apparently "disdya" means twice two-fold (also 4), so disdyakis triacontahedron, tetrakis rhombic triacontahedron, kisrhombic triacontahedron would all naturally mean the RT(S)EA setup.

[kbevers]

Converted to draft PR to safe-guard against a premature merge. There's a lot going on in this PR request. This looks impressive but from my very superficial first glance it seems as if there's reimplementations of already existing functionality. In particular, my spidey sense is triggered by authalic.h and the python code to create figures for the documentation. Is this Claude spewing out tons of code without consideration for the existing code base?

[jerstlouis]

@kbevers The authalic conversion stuff is something that I pointed out before in the codebase was already duplicated in multiple places, and which would be nice to re-use throughout (in multiple other projections). Also the version in this PR in the current stage I think is not taking into consideration the selected ellipsoid's flattening which it ideally should... trying to work with @felixpalmer on improving a few things.

[felixpalmer]

As noted in the changelist, the render scripts are purely to help reviewers recreate the images. The authalic conversion @jerstlouis already pointed out is likely implemented in a number of places in the codebase and I agree we should be leveraging the existing code. I would love to reduce the amount of code added, I'm aware that it is a lot, @jerstlouis is interested in collaborating on this

[felixpalmer]

@jerstlouis @kbevers I've removed the authalic.h file, replacing it with the built-in pj_authalic_... functions

[jerstlouis]

@rouault

Why not removing isea_legacy entirely ?

Sure. @felixpalmer Can we please do that?

it is not clear to me that having that in PROJ as 2D planar projection (without face numbers) is useful.

This PR should first be considered as fixes and improvements (performance and precision, likely fixing #4389) to the existing ISEA projection, and also introducing the underlying capability to support the vertex-oriented equal-area icosahedral projection (IVEA) which avoids the cusps of ISEA and is generally superior to ISEA for all practical purposes (@felixpalmer did you remove the +proj=dsea +net=icosahedron ? would that be re-introduced in a subsequent PR for +proj=ivea?).

There are many reasons for supporting these projections in PROJ:

• The Snyder equal-area map projection for polyhedral globes is generally useful as an equal-area projection with minimimal distortion outside the context of DGGS. While the large number of interruptions between disconnected faces might be seen as problematic, the flexible orientation and the ability to unfold the net in different configurations (I hope to make this more flexible in a subsequent PR by allowing custom parent list for +net) allows to reconfigure the nets to keep oceans or land connected and to represent a particular world perspective,
• DGGS applications can leverage this implementation of the projection. The faces can simply be cut-out from the 2D planar projection. For a DGGS based on the icosahedral symmetry, including hexagonal DGGRSs like ISEA3H / IVEA7H, or pentagonal like A5, this is particularly easy to do after applying an affine transformation to the 5x6 rhombic/icosahedral space to the icosahedral unfolding (which I hope to add as a future PR so that this can be done without +proj=pipeline which is not supported in several tools), where the rhombuses formed by two triangles become squares on 1x1 boundaries,
• By having these projections widely supported in GIS tools, we bridge the gap between DGGS and these tools, by using a traditional 2D CRS. This can notably be used to define a 2D Tile Matrix Set in the 5x6 space and store/exchange data in GeoPackages, WMTS, OGC API - Tiles etc. It also allows to properly reference vector or raster data converted from DGGS-(UB)JSON(-FG) (e.g., returned by OGC API - DGGS) to traditional formats like GeoTIFF / GeoJSON while avoiding re-projection,
• It provides support for these equal-area projections in tools like PostGIS to use them for custom hashing as per @ppKrauss 's use case at addressforall.

Could you please approve the workflows so that any remaining issue can be addressed? (maybe after @felixpalmer amends to remove isea_legacy). Thanks!

[felixpalmer]

@rouault I've removed the legacy ISEA code, just for comparison the new implementation is actually slightly less LoC that the previous one (934 rather than 954). As @jerstlouis says it can thus be viewed as a better refactor of the existing implementation.

what is the use for that in a library like PROJ?

To add to what @jerstlouis wrote, I believe such a family of projections would be useful in PROJ to provide a robust, canonical reference for others to build upon. ISEA is historically the most popular of the bunch, but other variants actually provide superior performance (e.g. DSEA / IVEA) but are not adopted due to lack of good implementations. I certainly would have been very grateful for such a reference when building the A5 DGGS!

I note that there is also good precedent (other than isea) for such unfolded polyhedral projections in PROJ, e.g. airocean and healpix.

PS @jerstlouis yes +proj=dsea +net=icosahedron (a.k.a IVEA) was removed with the view to follow-up, in order to not bloat this already large PR

[kbevers]

This is a massive PR and I am in no way confident in reviewing the projection implementation. The general state of the code looks better than when I last looked at it, but before I sink too much time into this I'd like a have the following clarified or reworked. At this stage I am mostly focused on documentation and overall structure of the code, as well as long-term maintainability.

Documentation and usage

1. Is a +proj=polyhedral generalisation available for users? docs/source/operations/projections/polyhedral.rst sort of hints at it but it's not clarified on the page. If it does some examples are further clarification is needed. If such a projection doesn't exist this doc page goes against the grain of all the others (every doc page in that chapter refers to a specific operation).
2. The figures for the docs are very helpful. I'd like to see them (re-)created using the existing scripts that we use for (almost) all the other similar figures. We might want to change the visuals of the figures at some point in the future, this let's us do that. See https://github.com/OSGeo/PROJ/tree/master/docs/plot.
3. The figures of the newly added projections should be included in docs/source/operations/projections/all_images.rst.

Code

1. I don't like the addition of the src/projections/polyhedral/ folder. I understand that it seperates code into smaller, logical chunks but it differs greatly from the existing practice of keeping everything in one file. As with the docs, one (class of) projection(s) in one file. Including the supporting bits.
2. Parts of the contents of src/projections/polyhedral/ is more generic than the polyhedral code. In particular the vector maths. This would be better suited in src/ so it can be used more easily in other parts of the code. I don't think an equivalent exists today but I am sure some of the same concepts are used throughout the code. Putting this in a central place would provide the basis for a later effort to consolidate vector maths across the code base.
3. test/unit/print_polyhedral_layout.cpp seems like left overs from a previous iteration of this PR

[ldesousa]

@rouault I would like to stress the merit of a reference implementation for this family of map projections. Beyond DGGS, they also provide a natural mechanism to minimise distortion in equal-area projections. This implementation is too an homage to the work of John Snyder.

@mpadillaruiz @brsr @joaocsmanuel @sunayana All of you have proper insight into this projection family. Would be great if some of you could provide an extra review to this PR.

[jjimenezshaw]

3. The figures of the newly added projections should be included in docs/source/operations/projections/all_images.rst.

Read https://github.com/OSGeo/PROJ/tree/master/docs/docbuild to see how to automatically generate documentation and that file in particular.

[jjimenezshaw]

@felixpalmer could you tell what was done with AI (Claude)?
Reviewing code generated with AI is different, and usually much more costly for the reviewer.

[felixpalmer]

@kbevers

Is a +proj=polyhedral generalisation available for users

No, that doc page is there to give a high level explanation of the general Snyder polyhedral approach, while tsea, isea and dsea are the specific projections added. We discussed with @jerstlouis about a generic polyhedral projection but that is out of scope for this PR

The figures for the docs are very helpful

Thanks for the pointer, figures added and regenerated using plotdefs.json. Also added to all_images.rst. I included all the net variants for completeness, I can remove some if you prefer.

I don't like the addition of the src/projections/polyhedral/ folder.

The split wasn't purely stylistic — the headers also delineate where existing open-source code came from (gl-matrix, A5, Snyder's construction). That said, I'm happy to inline them into a single polyhedral.cpp to match the per-projection convention. The total SLOC is already lower than the isea.cpp it replaces (−23, 957 → 934 excluding comments), so it stays reasonably digestible as one file.

is more generic than the polyhedral code. In particular the vector maths. This would be better suited in src

Agreed, I was surprised there was no such common code myself. I would propose moving just the vec3.h file, the vec3_helpers.h are more specific.

test/unit/print_polyhedral_layout.cpp seems like left overs

Removed. It was there to help generate net layout diagrams like this:

[felixpalmer]

@jjimenezshaw there isn't a clean split, I use it as a tool. So for refactors, design discussions etc.

The docs are all handwritten, anything that has a header indicating inclusion from another project was translated using AI but is not worth reviewing (this is well-tested deployed code). The rest I wrote myself, AI was used to assist, not to blindly generate.

[brsr]

I wonder if we can split this pull up into more easily digestible pieces. I think of this pull as having two conceptual parts:

1. Code for a generalized version of the triangular projection described in Snyder's 1988 paper.
2. Code for separating the sphere into specified triangular regions, defining corresponding triangular regions in the plane, and keeping track of which projection object transforms between which pair of spherical and planar triangles.

The projection in 1 can be used on its own, independent from any global partitioning scheme. It would have mandatory parameters +lat_1 +lon_1 +lat_2 +lon_2 +lat_3 +lon_3. This is an interesting projection in its own right and is worth making accessible by itself. (In my opinion. This is the part I worked the most on, so I may be biased.)

Regarding part 2, I notice that airocean contains similar but less flexible code. Sharing the part 2 procedure with it would open up a lot of new triangulations. airocean is fairly new and the author is still active, and I'm also fairly familiar with that projection. I don't want to widen the scope too much by introducing changes to airocean but I also don't want to knowingly create duplication.

This probably needs discussion.

[jerstlouis]

@brsr @felixpalmer @kbevers

The projection in 1 can be used on its own, independent from any global partitioning scheme. It would have mandatory parameters +lat_1 +lon_1 +lat_2 +lon_2 +lat_3 +lon_3. This is an interesting projection in its own right and is worth making accessible by itself.

Perhaps as a first step we could go back to splitting that more general code for a specific triangle into the common shared src/ tree, where e.g. airocean could re-use it? A future PR could then potentially expose that generic triangle projection requiring three anchor points parameters, if that is deemed useful?

(actually I think it is still split in the last version, but inside polyhedral/, as src/projections/polyhedral/snyder.h ?)

[jerstlouis]

@felixpalmer @kbevers

Would it be possible to update the PROJ scripts to be able to optionally show that red dot where the default projected 0,0 happens to be?

The previous images with the red dot were a very significant improvement showing where that happens to be compared to the previous version (and the images currently in this PR after the recent change). Not having this clearly illustrated is a major pain point in trying to use these projections.

[jerstlouis]

@felixpalmer

Applies a Conway meta operation to cut up the polyhedron and net faces into right-triangles

I don't think the meta operation necessarily results in right triangles? (though for the trigonometric formulation, the math is different as it cannot assume one angle is 90 degrees).
Taking the rhombic triacontahedron (RTEA), I think the meta operation does not yield right-triangles, but the triangles before the full split (splitting the RT faces into 4) yields right triangles (the 120 disdyakis triacontahedron triangles).

Wondering if you made the optimization yet where actually for the typical unfolding the full meta operation is not necessary (you only need it if the net explicitly unfolds the two different parts of the fanned triangles from face center to face vertices into a specially cut-out net -- just the kis operation is otherwise enough). The future custom nets would allow that, but none of the the current hardcoded configurations need that I think (only 12 projection triangles rather than 24 are necessary for the tetrahedron, and only 60 rather than 120 are needed for the dodecahedron and icosahedron)

[felixpalmer]

@brsr thanks for the input and your great work this PR has built on.

Long-term agree the partitioning code should be sharable with airocean / s2 / healpix, but pushing back on adding that here. The PR is already a lot and the cross-projection consolidation is its own design conversation. I've moved vec3.h and vec3_helpers.h (renamed sphere.h) to src/ to address @kbevers's "too much code in src/projections/" point. That gives future shared partitioning a place to pull primitives from, and snyder.h is already standalone, so the lift you're describing is partly done.

@jerstlouis red-dot: agree it's useful UX, but it's a plot.py change rather than just plotdefs, and applies to other projections with non-obvious origins too. One for a follow-up.

I don't think the meta operation necessarily results in right triangles?

If the face is a regular polygon it does. You're right that this won't work for RTEA but let's deal with that in the future if we add that projection