diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json index da0a85d041..faf64e5ebf 100644 --- a/docs/plot/plotdefs.json +++ b/docs/plot/plotdefs.json @@ -269,6 +269,54 @@ "res": "low", "type": "poly" }, + { + "filename": "dsea.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "dsea", + "projstring": "+proj=dsea", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, + { + "filename": "dsea_a5.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "dsea_a5", + "projstring": "+proj=dsea +net=a5", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, + { + "filename": "dsea_crescent.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "dsea_crescent", + "projstring": "+proj=dsea +net=crescent", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, + { + "filename": "dsea_flower.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "dsea_flower", + "projstring": "+proj=dsea +net=flower", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, { "filename": "airocean.png", "latmax": 90, @@ -647,6 +695,18 @@ "type": "line", "delta_cut": 1e6 }, + { + "filename": "isea_pole.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "isea_pole", + "projstring": "+proj=isea +orient=pole", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, { "filename": "kav5.png", "latmax": 90, @@ -1504,6 +1564,30 @@ "res": "low", "type": "poly" }, + { + "filename": "tsea.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "tsea", + "projstring": "+proj=tsea", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, + { + "filename": "tsea_star.png", + "latmax": 90, + "latmin": -90, + "lonmax": 180, + "lonmin": -180, + "name": "tsea_star", + "projstring": "+proj=tsea +net=star", + "res": "low", + "type": "line", + "delta_cut": 1e6 + }, { "filename": "ups.png", "latmax": 90, diff --git a/docs/source/operations/options/azi.rst b/docs/source/operations/options/azi.rst new file mode 100644 index 0000000000..e1f11c555e --- /dev/null +++ b/docs/source/operations/options/azi.rst @@ -0,0 +1,5 @@ +.. option:: +azi= + + Azimuth from the polyhedron's first vertex to its second, in degrees + clockwise. Together with ``+orient_lat`` and ``+orient_lon`` this fully + constrains the 3D pose of the polyhedron on the sphere. diff --git a/docs/source/operations/options/lat_0_polyhedral.rst b/docs/source/operations/options/lat_0_polyhedral.rst new file mode 100644 index 0000000000..40e1cfa0e9 --- /dev/null +++ b/docs/source/operations/options/lat_0_polyhedral.rst @@ -0,0 +1,18 @@ +.. option:: +lat_0= + + Geodetic latitude of the geographic point that should land at the + projected origin ``(0, 0)``. Used together with ``+lon_0``; the pair + translates the projected output without rotating the polyhedron. + ``+x_0`` / ``+y_0`` are applied on top of this translation. + + When unset, the projected origin falls on the unfold's root face + instead (its centroid in general, or its bounding-box centre for + ``+proj=isea``) rather than at a fixed lat / lon. + + .. note:: + The default convention is to interpret this value as decimal degrees. To + specify radians instead, follow the value with the "r" character. + + Example: `+lat_0=1.570796r` + + See :ref:`Projection Units ` for more information. diff --git a/docs/source/operations/options/lon_0_polyhedral.rst b/docs/source/operations/options/lon_0_polyhedral.rst new file mode 100644 index 0000000000..f7fa94c6da --- /dev/null +++ b/docs/source/operations/options/lon_0_polyhedral.rst @@ -0,0 +1,18 @@ +.. option:: +lon_0= + + Geodetic longitude of the geographic point that should land at the + projected origin ``(0, 0)``. Used together with ``+lat_0``; the pair + translates the projected output without rotating the polyhedron. + ``+x_0`` / ``+y_0`` are applied on top of this translation. + + When unset, the projected origin falls on the unfold's root face + instead (its centroid in general, or its bounding-box centre for + ``+proj=isea``) rather than at a fixed lat / lon. + + .. note:: + The default convention is to interpret this value as decimal degrees. To + specify radians instead, follow the value with the "r" character. + + Example: `+lon_0=1.570796r` + + See :ref:`Projection Units ` for more information. diff --git a/docs/source/operations/options/orient_lat.rst b/docs/source/operations/options/orient_lat.rst new file mode 100644 index 0000000000..7bc9606a92 --- /dev/null +++ b/docs/source/operations/options/orient_lat.rst @@ -0,0 +1,5 @@ +.. option:: +orient_lat= + + Geodetic latitude of the polyhedron's first vertex, in degrees. Together + with ``+orient_lon`` and ``+azi`` this fully constrains the 3D pose of + the polyhedron on the sphere. diff --git a/docs/source/operations/options/orient_lon.rst b/docs/source/operations/options/orient_lon.rst new file mode 100644 index 0000000000..b6b365105d --- /dev/null +++ b/docs/source/operations/options/orient_lon.rst @@ -0,0 +1,5 @@ +.. option:: +orient_lon= + + Longitude of the polyhedron's first vertex, in degrees. Together with + ``+orient_lat`` and ``+azi`` this fully constrains the 3D pose of the + polyhedron on the sphere. diff --git a/docs/source/operations/projections/all_images.rst b/docs/source/operations/projections/all_images.rst index d7a5eaa03b..0f6146e5a5 100644 --- a/docs/source/operations/projections/all_images.rst +++ b/docs/source/operations/projections/all_images.rst @@ -292,6 +292,50 @@ List of all projection images ******************************************************************************** +:ref:`dsea` + +.. figure:: ./images/dsea.png + :width: 500 px + :align: center + :alt: dsea + + +******************************************************************************** + + +:ref:`dsea` + +.. figure:: ./images/dsea_a5.png + :width: 500 px + :align: center + :alt: dsea + + +******************************************************************************** + + +:ref:`dsea` + +.. figure:: ./images/dsea_crescent.png + :width: 500 px + :align: center + :alt: dsea + + +******************************************************************************** + + +:ref:`dsea` + +.. figure:: ./images/dsea_flower.png + :width: 500 px + :align: center + :alt: dsea + + +******************************************************************************** + + :ref:`eck1` .. figure:: ./images/eck1.png @@ -655,6 +699,17 @@ List of all projection images ******************************************************************************** +:ref:`isea` + +.. figure:: ./images/isea_pole.png + :width: 500 px + :align: center + :alt: isea + + +******************************************************************************** + + :ref:`kav5` .. figure:: ./images/kav5.png @@ -1480,6 +1535,28 @@ List of all projection images ******************************************************************************** +:ref:`tsea` + +.. figure:: ./images/tsea.png + :width: 500 px + :align: center + :alt: tsea + + +******************************************************************************** + + +:ref:`tsea` + +.. figure:: ./images/tsea_star.png + :width: 500 px + :align: center + :alt: tsea + + +******************************************************************************** + + :ref:`ups` .. figure:: ./images/ups.png diff --git a/docs/source/operations/projections/dsea.rst b/docs/source/operations/projections/dsea.rst new file mode 100644 index 0000000000..908cb79a41 --- /dev/null +++ b/docs/source/operations/projections/dsea.rst @@ -0,0 +1,118 @@ +.. _dsea: + +******************************************************************************** +Dodecahedral Snyder Equal Area +******************************************************************************** + +.. versionadded:: 9.9 + +Snyder's equal-area mapping :cite:`Snyder1992` applied to the twelve pentagonal +faces of a regular dodecahedron and unfolded into a planar net. + +The dodecahedron is subdivided into 12 × 10 = 120 right sub-triangles, +and each sub-triangle is mapped independently using the area-preserving +Snyder construction. + +See :ref:`polyhedral` for the shared theory. + ++---------------------+----------------------------------------------------------+ +| **Classification** | Polyhedral, equal area | ++---------------------+----------------------------------------------------------+ +| **Available forms** | Forward and inverse, spherical and ellipsoidal | ++---------------------+----------------------------------------------------------+ +| **Defined area** | Global | ++---------------------+----------------------------------------------------------+ +| **Alias** | dsea | ++---------------------+----------------------------------------------------------+ +| **Domain** | 2D | ++---------------------+----------------------------------------------------------+ +| **Input type** | Geodetic coordinates | ++---------------------+----------------------------------------------------------+ +| **Output type** | Projected coordinates | ++---------------------+----------------------------------------------------------+ + + +.. figure:: ./images/dsea.png + :width: 500 px + :align: center + :alt: Dodecahedral Snyder Equal Area + + proj-string: ``+proj=dsea`` + + +Nets +################################################################################ + +The default net (shown above) follows Snyder's Figure 11. Three alternative +nets are available via ``+net=``: + +a5 +-------------------------------------------------------------------------------- + +Layout used by the `A5 index `_. The first 8 faces contain +the majority of the populated land mass. + +.. figure:: ./images/dsea_a5.png + :width: 500 px + :align: center + :alt: DSEA A5 net + + proj-string: ``+proj=dsea +net=a5`` + +crescent +-------------------------------------------------------------------------------- + +.. figure:: ./images/dsea_crescent.png + :width: 500 px + :align: center + :alt: DSEA crescent net + + proj-string: ``+proj=dsea +net=crescent`` + +flower +-------------------------------------------------------------------------------- + +.. figure:: ./images/dsea_flower.png + :width: 500 px + :align: center + :alt: DSEA flower net + + proj-string: ``+proj=dsea +net=flower`` + + +Parameters +################################################################################ + +.. note:: + All parameters are optional. + +.. option:: +net= + + Selects the planar unfolding. Accepted values: ``dsea``, ``a5``, + ``crescent``, ``flower``. + + *Defaults to* ``dsea``. + +.. include:: ../options/orient_lat.rst + +*Defaults to* ``atan((1 + 2·cos(36°))/2) ≈ 52.6226°``. + +.. include:: ../options/orient_lon.rst + +*Defaults to −36.0 (or −129.0 when* ``+net=a5`` *).* + +.. include:: ../options/azi.rst + +*Defaults to 240.0.* + +.. include:: ../options/lat_0_polyhedral.rst + +.. include:: ../options/lon_0_polyhedral.rst + +.. include:: ../options/x_0.rst + +.. include:: ../options/y_0.rst + +.. include:: ../options/ellps.rst + +.. include:: ../options/R.rst diff --git a/docs/source/operations/projections/generate_all_images.py b/docs/source/operations/projections/generate_all_images.py index d2cf4fefb8..e5140df1e6 100644 --- a/docs/source/operations/projections/generate_all_images.py +++ b/docs/source/operations/projections/generate_all_images.py @@ -9,6 +9,11 @@ "peirce_q_square" : "peirce_q", "peirce_q_diamond" : "peirce_q", "peirce_q_horizontal" : "peirce_q", + "dsea_a5" : "dsea", + "dsea_crescent" : "dsea", + "dsea_flower" : "dsea", + "isea_pole" : "isea", + "tsea_star" : "tsea", } with open(all_images_path, 'w') as all_images: diff --git a/docs/source/operations/projections/images/dsea.png b/docs/source/operations/projections/images/dsea.png new file mode 100644 index 0000000000..d302ff94db Binary files /dev/null and b/docs/source/operations/projections/images/dsea.png differ diff --git a/docs/source/operations/projections/images/dsea_a5.png b/docs/source/operations/projections/images/dsea_a5.png new file mode 100644 index 0000000000..b968f11474 Binary files /dev/null and b/docs/source/operations/projections/images/dsea_a5.png differ diff --git a/docs/source/operations/projections/images/dsea_crescent.png b/docs/source/operations/projections/images/dsea_crescent.png new file mode 100644 index 0000000000..089eec52a7 Binary files /dev/null and b/docs/source/operations/projections/images/dsea_crescent.png differ diff --git a/docs/source/operations/projections/images/dsea_flower.png b/docs/source/operations/projections/images/dsea_flower.png new file mode 100644 index 0000000000..dfb0856c93 Binary files /dev/null and b/docs/source/operations/projections/images/dsea_flower.png differ diff --git a/docs/source/operations/projections/images/isea.png b/docs/source/operations/projections/images/isea.png index f99c3d013e..440c34344f 100644 Binary files a/docs/source/operations/projections/images/isea.png and b/docs/source/operations/projections/images/isea.png differ diff --git a/docs/source/operations/projections/images/isea_pole.png b/docs/source/operations/projections/images/isea_pole.png new file mode 100644 index 0000000000..0920823cfc Binary files /dev/null and b/docs/source/operations/projections/images/isea_pole.png differ diff --git a/docs/source/operations/projections/images/tsea.png b/docs/source/operations/projections/images/tsea.png new file mode 100644 index 0000000000..e62ac8e2b2 Binary files /dev/null and b/docs/source/operations/projections/images/tsea.png differ diff --git a/docs/source/operations/projections/images/tsea_star.png b/docs/source/operations/projections/images/tsea_star.png new file mode 100644 index 0000000000..daecd643ed Binary files /dev/null and b/docs/source/operations/projections/images/tsea_star.png differ diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst index d206ef9123..298390a6f8 100644 --- a/docs/source/operations/projections/index.rst +++ b/docs/source/operations/projections/index.rst @@ -39,6 +39,7 @@ Projections map the spherical 3D space to a flat 2D space. comill crast denoy + dsea eck1 eck2 eck3 @@ -116,6 +117,7 @@ Projections map the spherical 3D space to a flat 2D space. pconic peirce_q poly + polyhedral putp1 putp2 putp3 @@ -147,6 +149,7 @@ Projections map the spherical 3D space to a flat 2D space. tobmerc tpeqd tpers + tsea ups urm5 urmfps diff --git a/docs/source/operations/projections/isea.rst b/docs/source/operations/projections/isea.rst index 686fe8fe73..fb4bdb5544 100644 --- a/docs/source/operations/projections/isea.rst +++ b/docs/source/operations/projections/isea.rst @@ -4,16 +4,19 @@ Icosahedral Snyder Equal Area ******************************************************************************** -Snyder's Icosahedral Equal Area map projection on an icosahedron polyhedral globe -offers relatively low scale and angular distortion. The equations involved are -relatively straight-forward. The interruptions remain a disadvantage, as with -any low-error projection system applied to the entire globe :cite:`Snyder1992`. -This projection is used as a basis for defining discrete global grid hierarchies. +Snyder's equal-area mapping :cite:`Snyder1992` applied to the twenty triangular +faces of a regular icosahedron and unfolded into a planar net. + +The icosahedron is subdivided into 20 × 6 = 120 right sub-triangles, and each +sub-triangle is mapped independently using the area-preserving Snyder +construction. + +See :ref:`polyhedral` for the shared theory. +---------------------+----------------------------------------------------------+ | **Classification** | Polyhedral, equal area | +---------------------+----------------------------------------------------------+ -| **Available forms** | Forward and inverse, spherical | +| **Available forms** | Forward and inverse, spherical and ellipsoidal | +---------------------+----------------------------------------------------------+ | **Defined area** | Global | +---------------------+----------------------------------------------------------+ @@ -34,59 +37,56 @@ This projection is used as a basis for defining discrete global grid hierarchies proj-string: ``+proj=isea`` -.. note:: - As the projection is only defined on a sphere, it should only be used - with a spherical ellipsoid e.g., ``+R=6371007.18091875`` for a sphere with the - authalic radius of the WGS84 ellipsoid. For mapping coordinates on the WGS84 - ellipsoid to the authalic sphere, the input latitude should be converted - from geodetic latitude to authalic latitude. A future version may - automatically perform this conversion when using a non-spherical ellipsoid. -Parameters +Orientations ################################################################################ -.. note:: All parameters are optional for the projection. +``isea`` ships a single net (Snyder's Figure 12, shown above). A second named +orientation is available via ``+orient=pole``, which places one icosahedron +vertex on the geographic north pole: -.. option:: +orient= +.. figure:: ./images/isea_pole.png + :width: 500 px + :align: center + :alt: ISEA pole orientation - Can be set to either ``isea`` or ``pole``. See Snyder's Figure 12 for pole orientation :cite:`Snyder1992`. + proj-string: ``+proj=isea +orient=pole`` - *Defaults to isea* -.. option:: +azi= +Parameters +################################################################################ - Azimuth. +.. note:: + All parameters are optional. - Not supported by the inverse. +.. option:: +orient= - *Defaults to 0.0* + Shorthand for two named orientations. Accepted values: ``isea``, ``pole``. + Equivalent to setting ``+orient_lat`` / ``+orient_lon`` / ``+azi`` + explicitly; individual ``+orient_*`` parameters still override. -.. option:: +aperture= + *Defaults to* ``isea``. - Not supported by the inverse. +.. include:: ../options/orient_lat.rst - *Defaults to 3.0* +*Defaults to ~58.40° geodetic (arctan(φ) ≈ 58.2825° authalic).* -.. option:: +resolution= +.. include:: ../options/orient_lon.rst - Not supported by the inverse. +*Defaults to 11.25° on a sphere, 11.20° on an ellipsoid.* - *Defaults to 4.0* +.. include:: ../options/azi.rst -.. option:: +mode= +*Defaults to 0.0.* - Can be either ``plane``, ``di``, ``dd`` or ``hex``. +.. include:: ../options/lat_0_polyhedral.rst - Only ``plane`` supported by the inverse. +.. include:: ../options/lon_0_polyhedral.rst - *Defaults to plane* +.. include:: ../options/x_0.rst -.. include:: ../options/lon_0.rst +.. include:: ../options/y_0.rst -.. include:: ../options/lat_0.rst +.. include:: ../options/ellps.rst .. include:: ../options/R.rst - -.. include:: ../options/x_0.rst - -.. include:: ../options/y_0.rst diff --git a/docs/source/operations/projections/polyhedral.rst b/docs/source/operations/projections/polyhedral.rst new file mode 100644 index 0000000000..1688aee0ce --- /dev/null +++ b/docs/source/operations/projections/polyhedral.rst @@ -0,0 +1,110 @@ +.. _polyhedral: + +******************************************************************************** +Polyhedral Snyder Equal Area +******************************************************************************** + +.. versionadded:: 9.9 + +.. figure:: ./images/dsea_a5.png + :width: 800 px + :align: center + :alt: Polyhedral projection (DSEA A5 net) + +The polyhedral Snyder equal-area projection maps the sphere onto the faces of +a polyhedron using the method described in :cite:`Snyder1992`. The projection +is area-preserving and supports both forward and inverse transforms on either +a sphere or an ellipsoid (input geodetic latitudes are converted to authalic +latitudes internally so the equal-area property is preserved). + +The implementation is generic: each polyhedral face is subdivided into right +sub-triangles, Snyder's equal-area mapping is applied independently on each +sub-triangle, and the planar results are laid out according to a *net* — a +2D unfolding of the polyhedron. + +Three projection aliases are registered: + +- :ref:`tsea` — Tetrahedral Snyder Equal Area +- :ref:`dsea` — Dodecahedral Snyder Equal Area +- :ref:`isea` — Icosahedral Snyder Equal Area + +Triangles +******************************************************************************** + +At its core, the projection is an area-preserving mapping from a spherical +right triangle ABC to a planar right triangle XYZ. Both triangles have one +distinguished vertex — the *apex* — which is listed first in either triple: +A on the sphere, X in the plane. The apex is shared by the sub-triangles +that fan out from a single face centre, so the sub-triangles meeting at the +apex tile a complete polyhedral face. + +The planar triangles can be chosen with any consistent shape, as long as the +edges around the shared apex match the way the corresponding edges meet on +the sphere. The freedom in choosing this planar shape — together with the +freedom in choosing how to glue the faces together — is what defines a +*net*. + +Polyhedra and nets +******************************************************************************** + +Each polyhedral projection is built from two compact definitions: + +- A *polyhedron*: a list of 3D vertices on the unit sphere plus a list of + face vertex indices (CCW from outside). By convention the first vertex is + always `{0, 0, 1}` +- A *net*: a parent tree over the faces. The indexing is 1-based, with `0` + reserved as the marker for the root face (that which has no parent). + For example, `{2, 3, 0, 3}` means: + + - the 1st face has the **2nd face** as its parent + - the 2nd face has the **3rd face** as its parent + - the 3rd face has the **no face** as its parent (root) + - the 4th face has the **3rd face** as its parent + +Projection specifications +******************************************************************************** + +Because the planar layout is decoupled from the spherical triangulation, the +same polyhedron can be unfolded into several different nets, which can be +selected by name using the `+net` parameter, e.g. `+proj=tsea +net=star` + + +The defaults for ``+orient_lat`` / ``+orient_lon`` / ``+azi`` rotate each +polyhedron into a sensible default pose (e.g. centred on Europe / Africa for +``dsea``; matching Snyder's published Figure 12 for ``isea``). When +``+lat_0`` / ``+lon_0`` are not supplied, the projected origin (0, 0) lands +on the centroid of the unfold's root face — except for ``isea``, which +uses the bounding-box centre of the root face instead so the unfolded net +sits centred in the output. + +Parameters +******************************************************************************** + +.. note:: + All parameters are optional. Each polyhedral projection ships with a + default orientation that places its vertices in symmetric, well-known + positions; the parameters below override that default. + +.. option:: +net= + + Selects the planar unfolding (net). + + The accepted values and defaults depend on the projection + +.. include:: ../options/orient_lat.rst + +.. include:: ../options/orient_lon.rst + +.. include:: ../options/azi.rst + +.. include:: ../options/lat_0_polyhedral.rst + +.. include:: ../options/lon_0_polyhedral.rst + +.. include:: ../options/x_0.rst + +.. include:: ../options/y_0.rst + +.. include:: ../options/ellps.rst + +.. include:: ../options/R.rst diff --git a/docs/source/operations/projections/tsea.rst b/docs/source/operations/projections/tsea.rst new file mode 100644 index 0000000000..59b4b6ef5c --- /dev/null +++ b/docs/source/operations/projections/tsea.rst @@ -0,0 +1,94 @@ +.. _tsea: + +******************************************************************************** +Tetrahedral Snyder Equal Area +******************************************************************************** + +.. versionadded:: 9.9 + +Snyder's equal-area mapping :cite:`Snyder1992` applied to the four faces of a +regular tetrahedron and unfolded into a planar net. + +The tetrahedron is subdivided into 4 × 6 = 24 right sub-triangles, and each +sub-triangle is mapped independently using the area-preserving Snyder +construction. + +See :ref:`polyhedral` for the shared theory. + ++---------------------+----------------------------------------------------------+ +| **Classification** | Polyhedral, equal area | ++---------------------+----------------------------------------------------------+ +| **Available forms** | Forward and inverse, spherical and ellipsoidal | ++---------------------+----------------------------------------------------------+ +| **Defined area** | Global | ++---------------------+----------------------------------------------------------+ +| **Alias** | tsea | ++---------------------+----------------------------------------------------------+ +| **Domain** | 2D | ++---------------------+----------------------------------------------------------+ +| **Input type** | Geodetic coordinates | ++---------------------+----------------------------------------------------------+ +| **Output type** | Projected coordinates | ++---------------------+----------------------------------------------------------+ + + +.. figure:: ./images/tsea.png + :width: 500 px + :align: center + :alt: Tetrahedral Snyder Equal Area + + proj-string: ``+proj=tsea`` + + +Nets +################################################################################ + +The default net (shown above) follows Snyder's Figure 8. An alternative +``+net=star`` layout is also available: + +.. figure:: ./images/tsea_star.png + :width: 500 px + :align: center + :alt: TSEA star net + + proj-string: ``+proj=tsea +net=star`` + +In the ``star`` layout the south-cap face sits at the centre, with its three +neighbours fanned around it. + + +Parameters +################################################################################ + +.. note:: + All parameters are optional. + +.. option:: +net= + + Selects the planar unfolding. Accepted values: ``tsea``, ``star``. + + *Defaults to* ``tsea``. + +.. include:: ../options/orient_lat.rst + +*Defaults to 90.0.* + +.. include:: ../options/orient_lon.rst + +*Defaults to 0.0.* + +.. include:: ../options/azi.rst + +*Defaults to 0.0.* + +.. include:: ../options/lat_0_polyhedral.rst + +.. include:: ../options/lon_0_polyhedral.rst + +.. include:: ../options/x_0.rst + +.. include:: ../options/y_0.rst + +.. include:: ../options/ellps.rst + +.. include:: ../options/R.rst diff --git a/docs/source/spelling_wordlist.txt b/docs/source/spelling_wordlist.txt index 3cfc298c3d..bdb1876a80 100644 --- a/docs/source/spelling_wordlist.txt +++ b/docs/source/spelling_wordlist.txt @@ -270,6 +270,7 @@ doxygennamespace doxygenpage doxygenstruct doxygentypedef +dsea duplicateBreak d'utilisation dx @@ -668,6 +669,7 @@ naptrans natearth natively nd +neighbours nell netcdf netCDF @@ -768,6 +770,7 @@ poder Poder polewards Polyconic +Polyhedra polynormal pre Pre @@ -957,6 +960,7 @@ tradeoff tran triangulations triaxial +tsea Trimetric Tripel Triptychial diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 400b9e38a1..75cdf42bac 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -12327,7 +12327,6 @@ PROJStringParser::Private::buildProjectedCRS(int iStep, "q," // urm5 "path,lsat," // lsat "W,M," // hammer - "aperture,resolution," // isea )) { double value = getNumericValue(param.value, &hasError); values.push_back(ParameterValue::create( diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index cfde386a61..dcddfbb907 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -35,7 +35,7 @@ set(SRC_LIBPROJ_PROJECTIONS projections/bipc.cpp projections/bonne.cpp projections/eqdc.cpp - projections/isea.cpp + projections/polyhedral.cpp projections/ccon.cpp projections/imw_p.cpp projections/krovak.cpp diff --git a/src/pj_list.h b/src/pj_list.h index b813f9899e..6c895de34c 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -37,6 +37,7 @@ PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)") PROJ_HEAD(defmodel, "Deformation model") PROJ_HEAD(deformation, "Kinematic grid shift") PROJ_HEAD(denoy, "Denoyer Semi-Elliptical") +PROJ_HEAD(dsea, "Dodecahedral Snyder Equal Area") PROJ_HEAD(airocean, "Airocean Fuller") PROJ_HEAD(eck1, "Eckert I") PROJ_HEAD(eck2, "Eckert II") @@ -169,6 +170,7 @@ PROJ_HEAD(tobmerc, "Tobler-Mercator") PROJ_HEAD(topocentric, "Geocentric/Topocentric conversion") PROJ_HEAD(tpeqd, "Two Point Equidistant") PROJ_HEAD(tpers, "Tilted perspective") +PROJ_HEAD(tsea, "Tetrahedral Snyder Equal Area") PROJ_HEAD(unitconvert, "Unit conversion") PROJ_HEAD(ups, "Universal Polar Stereographic") PROJ_HEAD(urm5, "Urmaev V") diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp deleted file mode 100644 index 144bedd53c..0000000000 --- a/src/projections/isea.cpp +++ /dev/null @@ -1,1428 +0,0 @@ -/* - The public domain code for the forward direction was initially - written by Nathan Wagner. - - The inverse projection was adapted from Java and eC by - Jérôme Jacovella-St-Louis, originally from the Franz-Benjamin Mocnik's ISEA - implementation found at - https://github.com/mocnik-science/geogrid/blob/master/ - src/main/java/org/giscience/utils/geogrid/projections/ISEAProjection.java - with the following license: - -------------------------------------------------------------------------- - MIT License - - Copyright (c) 2017-2019 Heidelberg University - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in - all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -/* The ISEA projection a projects a sphere on the icosahedron. Thereby the size - * of areas mapped to the icosahedron are preserved. Angles and distances are - * however slightly distorted. The angular distortion is below 17.27 degree, and - * the scale variation is less than 16.3 per cent. - * - * The projection has been proposed and has been described in detail by: - * - * John P. Snyder: An equal-area map projection for polyhedral globes. - * Cartographica, 29(1), 10–21, 1992. doi:10.3138/27H7-8K88-4882-1752 - * - * Another description and improvements can be found in: - * - * Erika Harrison, Ali Mahdavi-Amiri, and Faramarz Samavati: Optimization of - * inverse Snyder polyhedral projection. International Conference on Cyberworlds - * 2011. doi:10.1109/CW.2011.36 - * - * Erika Harrison, Ali Mahdavi-Amiri, and Faramarz Samavati: Analysis of inverse - * Snyder optimizations. In: Marina L. Gavrilova, and C. J. Kenneth Tan (Eds): - * Transactions on Computational Science XVI. Heidelberg, Springer, 2012. pp. - * 134–148. doi:10.1007/978-3-642-32663-9_8 - */ - -#include -#include -#include -#include -#include -#include - -#include -#include - -#include "proj.h" -#include "proj_internal.h" -#include - -#define DEG36 0.62831853071795864768 -#define DEG72 1.25663706143591729537 -#define DEG90 M_PI_2 -#define DEG108 1.88495559215387594306 -#define DEG120 2.09439510239319549229 -#define DEG144 2.51327412287183459075 -#define DEG180 M_PI - -/* sqrt(5)/M_PI */ -#define ISEA_SCALE 0.8301572857837594396028083 - -/* 26.565051177 degrees */ -#define V_LAT 0.46364760899944494524 - -// Latitude of center of top icosahedron faces -// atan((3 + sqrt(5)) / 4) = 52.6226318593487 degrees -#define E_RAD 0.91843818701052843323 - -// Latitude of center of faces mirroring top icosahedron face -// atan((3 - sqrt(5)) / 4) = 10.8123169635739 degrees -#define F_RAD 0.18871053078356206978 - -// #define phi ((1 + sqrt(5)) / 2) -// #define atanphi 1.01722196789785136772 - -// g: Spherical distance from center of polygon face to -// any of its vertices on the sphere -// g = F + 2 * atan(phi) - 90 deg -- sdc2vos -#define sdc2vos 0.6523581397843681859886783 - -#define tang 0.76393202250021030358019673567 // tan(sdc2vos) - -// theta (30 degrees) is plane angle between radius -// vector to center and adjacent edge of plane polygon - -#define tan30 0.57735026918962576450914878 // tan(DEG_TO_RAD * 30) -#define cotTheta (1.0 / tan30) - -// G: spherical angle between radius vector to center and adjacent edge -// of spherical polygon on the globe (36 degrees) -// cos(DEG_TO_RAD * 36) -#define cosG 0.80901699437494742410229341718281905886 -// sin(DEG_TO_RAD * 36) -#define sinG 0.587785252292473129168705954639072768597652 -// cos(g) -#define cosSDC2VoS 0.7946544722917661229596057297879189448539 - -#define sinGcosSDC2VoS (sinG * cosSDC2VoS) // sin G cos g - -#define SQRT3 1.73205080756887729352744634150587236694280525381038 -#define sin60 (SQRT3 / 2.0) -#define cos30 (SQRT3 / 2.0) - -// tang * sin(60 deg) -#define TABLE_G (tang * sin60) - -// (1 / (2 * sqrt(5)) + 1 / 6.0) * sqrt(M_PI * sqrt(3)) -#define RprimeOverR 0.9103832815095032 // R' / R - -/* H = 0.25 R tan g = */ -#define TABLE_H (0.25 * tang) - -/* in radians */ -#define ISEA_STD_LAT 1.01722196792335072101 -#define ISEA_STD_LONG .19634954084936207740 - -namespace { // anonymous namespace - -struct GeoPoint { - double lat, lon; -}; // In radians - -struct hex { - int iso; - long x, y, z; -}; -} // anonymous namespace - -/* y *must* be positive down as the xy /iso conversion assumes this */ -static void hex_xy(struct hex *h) { - if (!h->iso) - return; - if (h->x >= 0) { - h->y = -h->y - (h->x + 1) / 2; - } else { - /* need to round toward -inf, not toward zero, so x-1 */ - h->y = -h->y - h->x / 2; - } - h->iso = 0; -} - -static void hex_iso(struct hex *h) { - if (h->iso) - return; - - if (h->x >= 0) { - h->y = (-h->y - (h->x + 1) / 2); - } else { - /* need to round toward -inf, not toward zero, so x-1 */ - h->y = (-h->y - (h->x) / 2); - } - - h->z = -h->x - h->y; - h->iso = 1; -} - -static void hexbin2(double width, double x, double y, long *i, long *j) { - double z, rx, ry, rz; - double abs_dx, abs_dy, abs_dz; - long ix, iy, iz, s; - struct hex h; - - x = x / cos(30 * M_PI / 180.0); /* rotated X coord */ - y = y - x / 2.0; /* adjustment for rotated X */ - - /* adjust for actual hexwidth */ - if (width == 0) { - throw "Division by zero"; - } - x /= width; - y /= width; - - z = -x - y; - - rx = floor(x + 0.5); - ix = lround(rx); - ry = floor(y + 0.5); - iy = lround(ry); - rz = floor(z + 0.5); - iz = lround(rz); - if (fabs((double)ix + iy) > std::numeric_limits::max() || - fabs((double)ix + iy + iz) > std::numeric_limits::max()) { - throw "Integer overflow"; - } - - s = ix + iy + iz; - - if (s) { - abs_dx = fabs(rx - x); - abs_dy = fabs(ry - y); - abs_dz = fabs(rz - z); - - if (abs_dx >= abs_dy && abs_dx >= abs_dz) { - ix -= s; - } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) { - iy -= s; - } else { - iz -= s; - } - } - h.x = ix; - h.y = iy; - h.z = iz; - h.iso = 1; - - hex_xy(&h); - *i = h.x; - *j = h.y; -} - -#define numIcosahedronFaces 20 - -namespace { // anonymous namespace -enum isea_address_form { ISEA_PLANE, ISEA_Q2DI, ISEA_Q2DD, ISEA_HEX }; - -struct isea_sincos { - double s, c; -}; - -struct isea_pt { - double x, y; -}; - -} // anonymous namespace - -// distortion -// static double maximumAngularDistortion = 17.27; -// static double maximumScaleVariation = 1.163; -// static double minimumScaleVariation = .860; - -// Vertices of dodecahedron centered in icosahedron triangular faces -static const GeoPoint facesCenterDodecahedronVertices[numIcosahedronFaces] = { - {E_RAD, DEG_TO_RAD * -144}, {E_RAD, DEG_TO_RAD * -72}, - {E_RAD, DEG_TO_RAD * 0}, {E_RAD, DEG_TO_RAD * 72}, - {E_RAD, DEG_TO_RAD * 144}, {F_RAD, DEG_TO_RAD * -144}, - {F_RAD, DEG_TO_RAD * -72}, {F_RAD, DEG_TO_RAD * 0}, - {F_RAD, DEG_TO_RAD * 72}, {F_RAD, DEG_TO_RAD * 144}, - {-F_RAD, DEG_TO_RAD * -108}, {-F_RAD, DEG_TO_RAD * -36}, - {-F_RAD, DEG_TO_RAD * 36}, {-F_RAD, DEG_TO_RAD * 108}, - {-F_RAD, DEG_TO_RAD * 180}, {-E_RAD, DEG_TO_RAD * -108}, - {-E_RAD, DEG_TO_RAD * -36}, {-E_RAD, DEG_TO_RAD * 36}, - {-E_RAD, DEG_TO_RAD * 108}, {-E_RAD, DEG_TO_RAD * 180}}; - -// NOTE: Very similar to ISEAPlanarProjection::faceOrientation(), -// but the forward projection sometimes is returning a negative M_PI -static inline double az_adjustment(int triangle) { - if ((triangle >= 5 && triangle <= 9) || triangle == 15 || triangle == 16) - return M_PI; - else if (triangle >= 17) - return -M_PI; - return 0; -} - -static struct isea_pt isea_triangle_xy(int triangle) { - struct isea_pt c; - - triangle %= numIcosahedronFaces; - - c.x = TABLE_G * ((triangle % 5) - 2) * 2.0; - if (triangle > 9) { - c.x += TABLE_G; - } - - // REVIEW: This is likely related to - // pj_isea_data::yOffsets - switch (triangle / 5) { - case 0: - c.y = 5.0 * TABLE_H; - break; - case 1: - c.y = TABLE_H; - break; - case 2: - c.y = -TABLE_H; - break; - case 3: - c.y = -5.0 * TABLE_H; - break; - default: - /* should be impossible */ - exit(EXIT_FAILURE); - } - c.x *= RprimeOverR; - c.y *= RprimeOverR; - - return c; -} - -namespace { // anonymous namespace - -class ISEAPlanarProjection; - -struct pj_isea_data { - double o_lat, o_lon, o_az; /* orientation, radians */ - int aperture; /* valid values depend on partitioning method */ - int resolution; - isea_address_form output; /* an isea_address_form */ - int triangle; /* triangle of last transformed point */ - int quad; /* quad of last transformed point */ - isea_sincos vertexLatSinCos[numIcosahedronFaces]; - - double R2; - double Rprime; - double Rprime2X; - double RprimeTang; - double Rprime2Tan2g; - double triTang; - double centerToBase; - double triWidth; - double yOffsets[4]; - double xo, yo; - double sx, sy; - ISEAPlanarProjection *p; - - void initialize(const PJ *P); -}; -} // anonymous namespace - -#ifdef _MSC_VER -#pragma warning(push) -/* disable unreachable code warning for return 0 */ -#pragma warning(disable : 4702) -#endif - -#define SAFE_ARC_EPSILON 1E-15 - -static inline double safeArcSin(double t) { - return fabs(t) < SAFE_ARC_EPSILON ? 0 - : fabs(t - 1.0) < SAFE_ARC_EPSILON ? M_PI / 2 - : fabs(t + 1.0) < SAFE_ARC_EPSILON ? -M_PI / 2 - : asin(t); -} - -static inline double safeArcCos(double t) { - return fabs(t) < SAFE_ARC_EPSILON ? M_PI / 2 - : fabs(t + 1) < SAFE_ARC_EPSILON ? M_PI - : fabs(t - 1) < SAFE_ARC_EPSILON ? 0 - : acos(t); -} - -#undef SAFE_ARC_EPSILON - -/* coord needs to be in radians */ -static int isea_snyder_forward(const struct pj_isea_data *data, - const struct GeoPoint *ll, struct isea_pt *out) { - int i; - double sinLat = sin(ll->lat), cosLat = cos(ll->lat); - - /* - * TODO by locality of reference, start by trying the same triangle - * as last time - */ - for (i = 0; i < numIcosahedronFaces; i++) { - /* additional variables from snyder */ - double q, H, Ag, Azprime, Az, dprime, f, rho, x, y; - - /* variables used to store intermediate results */ - double az_offset; - - /* how many multiples of 60 degrees we adjust the azimuth */ - int Az_adjust_multiples; - - const struct GeoPoint *center = &facesCenterDodecahedronVertices[i]; - const struct isea_sincos *centerLatSinCos = &data->vertexLatSinCos[i]; - double dLon = ll->lon - center->lon; - double cosLat_cosLon = cosLat * cos(dLon); - double cosZ = - centerLatSinCos->s * sinLat + centerLatSinCos->c * cosLat_cosLon; - double sinAz, cosAz; - - /* step 1 */ - double z = safeArcCos(cosZ); - - /* not on this triangle */ - if (z > sdc2vos /*g*/ + 0.000005) { /* TODO DBL_EPSILON */ - continue; - } - - /* snyder eq 14 */ - Az = atan2(cosLat * sin(dLon), centerLatSinCos->c * sinLat - - centerLatSinCos->s * cosLat_cosLon); - - /* step 2 */ - - /* This calculates "some" vertex coordinate */ - az_offset = az_adjustment(i); - - Az -= az_offset; - - /* TODO I don't know why we do this. It's not in snyder */ - /* maybe because we should have picked a better vertex */ - if (Az < 0.0) { - Az += 2.0 * M_PI; - } - /* - * adjust Az for the point to fall within the range of 0 to - * 2(90 - theta) or 60 degrees for the hexagon, by - * and therefore 120 degrees for the triangle - * of the icosahedron - * subtracting or adding multiples of 60 degrees to Az and - * recording the amount of adjustment - */ - - Az_adjust_multiples = 0; - while (Az < 0.0) { - Az += DEG120; - Az_adjust_multiples--; - } - while (Az > DEG120 + DBL_EPSILON) { - Az -= DEG120; - Az_adjust_multiples++; - } - - /* step 3 */ - - /* Calculate q from eq 9. */ - cosAz = cos(Az); - sinAz = sin(Az); - q = atan2(tang, cosAz + sinAz * cotTheta); - - /* not in this triangle */ - if (z > q + 0.000005) { - continue; - } - /* step 4 */ - - /* Apply equations 5-8 and 10-12 in order */ - - /* eq 5 */ - /* R' in the paper is for the truncated (icosahedron?) */ - - /* eq 6 */ - H = acos(sinAz * sinGcosSDC2VoS /* sin(G) * cos(g) */ - cosAz * cosG); - - /* eq 7 */ - /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */ - Ag = Az + DEG_TO_RAD * 36 /* G */ + H - DEG180; - - /* eq 8 */ - Azprime = atan2(2.0 * Ag, RprimeOverR * RprimeOverR * tang * tang - - 2.0 * Ag * cotTheta); - - /* eq 10 */ - /* cot(theta) = 1.73205080756887729355 */ - dprime = RprimeOverR * tang / (cos(Azprime) + sin(Azprime) * cotTheta); - - /* eq 11 */ - f = dprime / (2.0 * RprimeOverR * sin(q / 2.0)); - - /* eq 12 */ - rho = 2.0 * RprimeOverR * f * sin(z / 2.0); - - /* - * add back the same 60 degree multiple adjustment from step - * 2 to Azprime - */ - - Azprime += DEG120 * Az_adjust_multiples; - - /* calculate rectangular coordinates */ - - x = rho * sin(Azprime); - y = rho * cos(Azprime); - - /* - * TODO - * translate coordinates to the origin for the particular - * hexagon on the flattened polyhedral map plot - */ - - out->x = x; - out->y = y; - - return i; - } - - /* - * should be impossible, this implies that the coordinate is not on - * any triangle - */ - - fprintf(stderr, "impossible transform: %f %f is not on any triangle\n", - PJ_TODEG(ll->lon), PJ_TODEG(ll->lat)); - - exit(EXIT_FAILURE); -} - -#ifdef _MSC_VER -#pragma warning(pop) -#endif - -/* - * return the new coordinates of any point in original coordinate system. - * Define a point (newNPold) in original coordinate system as the North Pole in - * new coordinate system, and the great circle connect the original and new - * North Pole as the lon0 longitude in new coordinate system, given any point - * in original coordinate system, this function return the new coordinates. - */ - -/* formula from Snyder, Map Projections: A working manual, p31 */ -/* - * old north pole at np in new coordinates - * could be simplified a bit with fewer intermediates - * - * TODO take a result pointer - */ -static struct GeoPoint snyder_ctran(const struct GeoPoint &np, - const struct GeoPoint &pt) { - struct GeoPoint result; - double phi = pt.lat, lambda = pt.lon; - double alpha = np.lat, beta = np.lon; - double dlambda = lambda - beta /* lambda0 */; - double cos_p = cos(phi), sin_p = sin(phi); - double cos_a = cos(alpha), sin_a = sin(alpha); - double cos_dlambda = cos(dlambda), sin_dlambda = sin(dlambda); - - /* mpawm 5-7 */ - double sin_phip = sin_a * sin_p - cos_a * cos_p * cos_dlambda; - - /* mpawm 5-8b */ - - /* use the two argument form so we end up in the right quadrant */ - double lp_b = /* lambda prime minus beta */ - atan2(cos_p * sin_dlambda, sin_a * cos_p * cos_dlambda + cos_a * sin_p); - double lambdap = lp_b + beta; - - /* normalize longitude */ - /* TODO can we just do a modulus ? */ - lambdap = fmod(lambdap, 2 * M_PI); - while (lambdap > M_PI) - lambdap -= 2 * M_PI; - while (lambdap < -M_PI) - lambdap += 2 * M_PI; - - result.lat = safeArcSin(sin_phip); - result.lon = lambdap; - return result; -} - -static struct GeoPoint isea_ctran(const struct GeoPoint *np, - const struct GeoPoint *pt, double lon0) { - struct GeoPoint cnp = {np->lat, np->lon + M_PI}; - struct GeoPoint npt = snyder_ctran(cnp, *pt); - - npt.lon -= (/* M_PI */ -lon0 + np->lon); - /* - * snyder is down tri 3, isea is along side of tri1 from vertex 0 to - * vertex 1 these are 180 degrees apart - */ - // npt.lon += M_PI; - - /* normalize lon */ - npt.lon = fmod(npt.lon, 2 * M_PI); - while (npt.lon > M_PI) - npt.lon -= 2 * M_PI; - while (npt.lon < -M_PI) - npt.lon += 2 * M_PI; - - return npt; -} - -/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */ - -static int isea_grid_init(struct pj_isea_data *g) { - int i; - - if (!g) - return 0; - - g->o_lat = ISEA_STD_LAT; - g->o_lon = ISEA_STD_LONG; - g->o_az = 0.0; - g->aperture = 4; - g->resolution = 6; - - for (i = 0; i < numIcosahedronFaces; i++) { - const GeoPoint *c = &facesCenterDodecahedronVertices[i]; - g->vertexLatSinCos[i].s = sin(c->lat); - g->vertexLatSinCos[i].c = cos(c->lat); - } - return 1; -} - -static void isea_orient_isea(struct pj_isea_data *g) { - if (!g) - return; - g->o_lat = ISEA_STD_LAT; - g->o_lon = ISEA_STD_LONG; - g->o_az = 0.0; -} - -static void isea_orient_pole(struct pj_isea_data *g) { - if (!g) - return; - g->o_lat = M_PI / 2.0; - g->o_lon = 0.0; - g->o_az = 0; -} - -static int isea_transform(struct pj_isea_data *g, struct GeoPoint *in, - struct isea_pt *out) { - struct GeoPoint i, pole; - int tri; - - pole.lat = g->o_lat; - pole.lon = g->o_lon; - - i = isea_ctran(&pole, in, g->o_az); - - tri = isea_snyder_forward(g, &i, out); - g->triangle = tri; - - return tri; -} - -#define DOWNTRI(tri) ((tri / 5) % 2 == 1) - -static void isea_rotate(struct isea_pt *pt, double degrees) { - double rad; - - double x, y; - - rad = -degrees * M_PI / 180.0; - while (rad >= 2.0 * M_PI) - rad -= 2.0 * M_PI; - while (rad <= -2.0 * M_PI) - rad += 2.0 * M_PI; - - x = pt->x * cos(rad) + pt->y * sin(rad); - y = -pt->x * sin(rad) + pt->y * cos(rad); - - pt->x = x; - pt->y = y; -} - -static void isea_tri_plane(int tri, struct isea_pt *pt) { - struct isea_pt tc; /* center of triangle */ - - if (DOWNTRI(tri)) { - pt->x *= -1; - pt->y *= -1; - } - tc = isea_triangle_xy(tri); - pt->x += tc.x; - pt->y += tc.y; -} - -/* convert projected triangle coords to quad xy coords, return quad number */ -static int isea_ptdd(int tri, struct isea_pt *pt) { - int downtri, quadz; - - downtri = ((tri / 5) % 2 == 1); - quadz = (tri % 5) + (tri / 10) * 5 + 1; - - // NOTE: This would always be a 60 degrees rotation if the flip were - // already done as in isea_tri_plane() - isea_rotate(pt, downtri ? 240.0 : 60.0); - if (downtri) { - pt->x += 0.5; - /* pt->y += cos(30.0 * M_PI / 180.0); */ - pt->y += cos30; - } - return quadz; -} - -static int isea_dddi_ap3odd(struct pj_isea_data *g, int quadz, - struct isea_pt *pt, struct isea_pt *di) { - struct isea_pt v; - double hexwidth; - double sidelength; /* in hexes */ - long d, i; - long maxcoord; - struct hex h; - - /* This is the number of hexes from apex to base of a triangle */ - sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0; - - /* apex to base is cos(30deg) */ - hexwidth = cos(M_PI / 6.0) / sidelength; - - /* TODO I think sidelength is always x.5, so - * (int)sidelength * 2 + 1 might be just as good - */ - maxcoord = lround((sidelength * 2.0)); - - v = *pt; - hexbin2(hexwidth, v.x, v.y, &h.x, &h.y); - h.iso = 0; - hex_iso(&h); - - d = h.x - h.z; - i = h.x + h.y + h.y; - - /* - * you want to test for max coords for the next quad in the same - * "row" first to get the case where both are max - */ - if (quadz <= 5) { - if (d == 0 && i == maxcoord) { - /* north pole */ - quadz = 0; - d = 0; - i = 0; - } else if (i == maxcoord) { - /* upper right in next quad */ - quadz += 1; - if (quadz == 6) - quadz = 1; - i = maxcoord - d; - d = 0; - } else if (d == maxcoord) { - /* lower right in quad to lower right */ - quadz += 5; - d = 0; - } - } else /* if (quadz >= 6) */ { - if (i == 0 && d == maxcoord) { - /* south pole */ - quadz = 11; - d = 0; - i = 0; - } else if (d == maxcoord) { - /* lower right in next quad */ - quadz += 1; - if (quadz == 11) - quadz = 6; - d = maxcoord - i; - i = 0; - } else if (i == maxcoord) { - /* upper right in quad to upper right */ - quadz = (quadz - 4) % 5; - i = 0; - } - } - - di->x = d; - di->y = i; - - g->quad = quadz; - return quadz; -} - -static int isea_dddi(struct pj_isea_data *g, int quadz, struct isea_pt *pt, - struct isea_pt *di) { - struct isea_pt v; - double hexwidth; - long sidelength; /* in hexes */ - struct hex h; - - if (g->aperture == 3 && g->resolution % 2 != 0) { - return isea_dddi_ap3odd(g, quadz, pt, di); - } - /* todo might want to do this as an iterated loop */ - if (g->aperture > 0) { - double sidelengthDouble = pow(g->aperture, g->resolution / 2.0); - if (fabs(sidelengthDouble) > std::numeric_limits::max()) { - throw "Integer overflow"; - } - sidelength = lround(sidelengthDouble); - } else { - sidelength = g->resolution; - } - - if (sidelength == 0) { - throw "Division by zero"; - } - hexwidth = 1.0 / sidelength; - - v = *pt; - isea_rotate(&v, -30.0); - hexbin2(hexwidth, v.x, v.y, &h.x, &h.y); - h.iso = 0; - hex_iso(&h); - - /* we may actually be on another quad */ - if (quadz <= 5) { - if (h.x == 0 && h.z == -sidelength) { - /* north pole */ - quadz = 0; - h.z = 0; - h.y = 0; - h.x = 0; - } else if (h.z == -sidelength) { - quadz = quadz + 1; - if (quadz == 6) - quadz = 1; - h.y = sidelength - h.x; - h.z = h.x - sidelength; - h.x = 0; - } else if (h.x == sidelength) { - quadz += 5; - h.y = -h.z; - h.x = 0; - } - } else /* if (quadz >= 6) */ { - if (h.z == 0 && h.x == sidelength) { - /* south pole */ - quadz = 11; - h.x = 0; - h.y = 0; - h.z = 0; - } else if (h.x == sidelength) { - quadz = quadz + 1; - if (quadz == 11) - quadz = 6; - h.x = h.y + sidelength; - h.y = 0; - h.z = -h.x; - } else if (h.y == -sidelength) { - quadz -= 4; - h.y = 0; - h.z = -h.x; - } - } - di->x = h.x; - di->y = -h.z; - - g->quad = quadz; - return quadz; -} - -static int isea_ptdi(struct pj_isea_data *g, int tri, struct isea_pt *pt, - struct isea_pt *di) { - struct isea_pt v; - int quadz; - - v = *pt; - quadz = isea_ptdd(tri, &v); - quadz = isea_dddi(g, quadz, &v, di); - return quadz; -} - -/* TODO just encode the quad in the d or i coordinate - * quad is 0-11, which can be four bits. - * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf - */ -/* convert a q2di to global hex coord */ -static int isea_hex(struct pj_isea_data *g, int tri, struct isea_pt *pt, - struct isea_pt *hex) { - struct isea_pt v; -#ifdef FIXME - long sidelength; - long d, i, x, y; -#endif - int quadz; - - quadz = isea_ptdi(g, tri, pt, &v); - - if (v.x < (INT_MIN >> 4) || v.x > (INT_MAX >> 4)) { - throw "Invalid shift"; - } - hex->x = ((int)v.x * 16) + quadz; - hex->y = v.y; - - return 1; -#ifdef FIXME - d = lround(floor(v.x)); - i = lround(floor(v.y)); - - /* Aperture 3 odd resolutions */ - if (g->aperture == 3 && g->resolution % 2 != 0) { - long offset = lround((pow(3.0, g->resolution - 1) + 0.5)); - - d += offset * ((g->quadz - 1) % 5); - i += offset * ((g->quadz - 1) % 5); - - if (quadz == 0) { - d = 0; - i = offset; - } else if (quadz == 11) { - d = 2 * offset; - i = 0; - } else if (quadz > 5) { - d += offset; - } - - x = (2 * d - i) / 3; - y = (2 * i - d) / 3; - - hex->x = x + offset / 3; - hex->y = y + 2 * offset / 3; - return 1; - } - - /* aperture 3 even resolutions and aperture 4 */ - sidelength = lround((pow(g->aperture, g->resolution / 2.0))); - if (g->quad == 0) { - hex->x = 0; - hex->y = sidelength; - } else if (g->quad == 11) { - hex->x = sidelength * 2; - hex->y = 0; - } else { - hex->x = d + sidelength * ((g->quad - 1) % 5); - if (g->quad > 5) - hex->x += sidelength; - hex->y = i + sidelength * ((g->quad - 1) % 5); - } - - return 1; -#endif -} - -static struct isea_pt isea_forward(struct pj_isea_data *g, - struct GeoPoint *in) { - isea_pt out; - int tri = isea_transform(g, in, &out); - - if (g->output == ISEA_PLANE) - isea_tri_plane(tri, &out); - else { - isea_pt coord; - - /* convert to isea standard triangle size */ - out.x *= ISEA_SCALE; // / g->radius; - out.y *= ISEA_SCALE; // / g->radius; - out.x += 0.5; - out.y += 2.0 * .14433756729740644112; - - switch (g->output) { - case ISEA_PLANE: - /* already handled above -- GCC should not be complaining */ - case ISEA_Q2DD: - /* Same as above, we just don't print as much */ - g->quad = isea_ptdd(tri, &out); - break; - case ISEA_Q2DI: - g->quad = isea_ptdi(g, tri, &out, &coord); - return coord; - case ISEA_HEX: - isea_hex(g, tri, &out, &coord); - return coord; - } - } - return out; -} - -/* - * Proj 4 integration code follows - */ - -PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph"; - -static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ - PJ_XY xy = {0.0, 0.0}; - struct pj_isea_data *Q = static_cast(P->opaque); - struct isea_pt out; - struct GeoPoint in; - - // TODO: Convert geodetic latitude to authalic latitude if not - // spherical as in eqearth, healpix, laea, etc. - in.lat = lp.phi; - in.lon = lp.lam; - - try { - out = isea_forward(Q, &in); - } catch (const char *) { - proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); - return proj_coord_error().xy; - } - - xy.x = out.x; - xy.y = out.y; - - return xy; -} - -static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P); - -PJ *PJ_PROJECTION(isea) { - char *opt; - struct pj_isea_data *Q = static_cast( - calloc(1, sizeof(struct pj_isea_data))); - if (nullptr == Q) - return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); - P->opaque = Q; - - // NOTE: if a inverse was needed, there is some material at - // https://brsr.github.io/2021/08/31/snyder-equal-area.html - P->fwd = isea_s_forward; - P->inv = isea_s_inverse; - isea_grid_init(Q); - Q->output = ISEA_PLANE; - - /* P->radius = P->a; / * otherwise defaults to 1 */ - /* calling library will scale, I think */ - - opt = pj_param(P->ctx, P->params, "sorient").s; - if (opt) { - if (!strcmp(opt, "isea")) { - isea_orient_isea(Q); - } else if (!strcmp(opt, "pole")) { - isea_orient_pole(Q); - } else { - proj_log_error( - P, - _("Invalid value for orient: only isea or pole are supported")); - return pj_default_destructor(P, - PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); - } - } - - if (pj_param(P->ctx, P->params, "tazi").i) { - Q->o_az = pj_param(P->ctx, P->params, "razi").f; - } - - if (pj_param(P->ctx, P->params, "tlon_0").i) { - Q->o_lon = pj_param(P->ctx, P->params, "rlon_0").f; - } - - if (pj_param(P->ctx, P->params, "tlat_0").i) { - Q->o_lat = pj_param(P->ctx, P->params, "rlat_0").f; - } - - opt = pj_param(P->ctx, P->params, "smode").s; - if (opt) { - if (!strcmp(opt, "plane")) { - Q->output = ISEA_PLANE; - } else if (!strcmp(opt, "di")) { - Q->output = ISEA_Q2DI; - } else if (!strcmp(opt, "dd")) { - Q->output = ISEA_Q2DD; - } else if (!strcmp(opt, "hex")) { - Q->output = ISEA_HEX; - } else { - proj_log_error(P, _("Invalid value for mode: only plane, di, dd or " - "hex are supported")); - return pj_default_destructor(P, - PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); - } - } - - /* REVIEW: Was this an undocumented +rescale= parameter? - if (pj_param(P->ctx, P->params, "trescale").i) { - Q->radius = ISEA_SCALE; - } - */ - - if (pj_param(P->ctx, P->params, "tresolution").i) { - Q->resolution = pj_param(P->ctx, P->params, "iresolution").i; - } else { - Q->resolution = 4; - } - - if (pj_param(P->ctx, P->params, "taperture").i) { - Q->aperture = pj_param(P->ctx, P->params, "iaperture").i; - } else { - Q->aperture = 3; - } - - Q->initialize(P); - - return P; -} - -#define Min std::min -#define Max std::max - -#define inf std::numeric_limits::infinity() - -// static define precision = DEG_TO_RAD * 1e-9; -#define precision (DEG_TO_RAD * 1e-11) -#define precisionPerDefinition (DEG_TO_RAD * 1e-5) - -#define AzMax (DEG_TO_RAD * 120) - -#define westVertexLon (DEG_TO_RAD * -144) - -namespace { // anonymous namespace - -struct ISEAFacePoint { - int face; - double x, y; -}; - -class ISEAPlanarProjection { - public: - explicit ISEAPlanarProjection(const GeoPoint &value) - : orientation(value), cosOrientationLat(cos(value.lat)), - sinOrientationLat(sin(value.lat)) {} - - bool cartesianToGeo(const PJ_XY &inPosition, const pj_isea_data *params, - GeoPoint &result) { - bool r = false; - static const double epsilon = 1E-11; - int face = 0; - PJ_XY position = inPosition; - -#define sr -sin60 // sin(-60) -#define cr 0.5 // cos(-60) - if (position.x < 0 || - (position.x < params->triWidth / 2 && position.y < 0 && - position.y * cr < position.x * sr)) - position.x += 5 * params->triWidth; // Wrap around -// Rotate and shear to determine face if not stored in position.z -#define shearX (1.0 / SQRT3) - double yp = -(position.x * sr + position.y * cr); - double x = - (position.x * cr - position.y * sr + yp * shearX) * params->sx; - double y = yp * params->sy; -#undef shearX -#undef sr -#undef cr - - if (x < 0 || (y > x && x < 5 - epsilon)) - x += epsilon; - else if (x > 5 || (y < x && x > 0 + epsilon)) - x -= epsilon; - - if (y < 0 || (x > y && y < 6 - epsilon)) - y += epsilon; - else if (y > 6 || (x < y && y > 0 + epsilon)) - y -= epsilon; - - if (x >= 0 && x <= 5 && y >= 0 && y <= 6) { - int ix = Max(0, Min(4, (int)x)); - int iy = Max(0, Min(5, (int)y)); - - if (iy == ix || iy == ix + 1) { - int rhombus = ix + iy; - bool top = x - ix > y - iy; - face = -1; - - switch (rhombus) { - case 0: - face = top ? 0 : 5; - break; - case 2: - face = top ? 1 : 6; - break; - case 4: - face = top ? 2 : 7; - break; - case 6: - face = top ? 3 : 8; - break; - case 8: - face = top ? 4 : 9; - break; - - case 1: - face = top ? 10 : 15; - break; - case 3: - face = top ? 11 : 16; - break; - case 5: - face = top ? 12 : 17; - break; - case 7: - face = top ? 13 : 18; - break; - case 9: - face = top ? 14 : 19; - break; - } - face++; - } - } - - if (face) { - int fy = (face - 1) / 5, fx = (face - 1) - 5 * fy; - double rx = - position.x - (2 * fx + fy / 2 + 1) * params->triWidth / 2; - double ry = - position.y - (params->yOffsets[fy] + 3 * params->centerToBase); - GeoPoint dst; - - r = icosahedronToSphere({face - 1, rx, ry}, params, dst); - - if (dst.lon < -M_PI - epsilon) - dst.lon += 2 * M_PI; - else if (dst.lon > M_PI + epsilon) - dst.lon -= 2 * M_PI; - - result = {dst.lat, dst.lon}; - } - return r; - } - - // Converts coordinates on the icosahedron to geographic coordinates - // (inverse projection) - bool icosahedronToSphere(const ISEAFacePoint &c, const pj_isea_data *params, - GeoPoint &r) { - if (c.face >= 0 && c.face < numIcosahedronFaces) { - double Az = atan2(c.x, c.y); // Az' - double rho = sqrt(c.x * c.x + c.y * c.y); - double AzAdjustment = faceOrientation(c.face); - - Az += AzAdjustment; - while (Az < 0) { - AzAdjustment += AzMax; - Az += AzMax; - } - while (Az > AzMax) { - AzAdjustment -= AzMax; - Az -= AzMax; - } - { - double sinAz = sin(Az), cosAz = cos(Az); - double cotAz = cosAz / sinAz; - double area = params->Rprime2Tan2g / - (2 * (cotAz + cotTheta)); // A_G or A_{ABD} - double deltaAz = 10 * precision; - double degAreaOverR2Plus180Minus36 = - area / params->R2 - westVertexLon; - double Az_earth = Az; - - while (fabs(deltaAz) > precision) { - double sinAzEarth = sin(Az_earth), - cosAzEarth = cos(Az_earth); - double H = - acos(sinAzEarth * sinGcosSDC2VoS - cosAzEarth * cosG); - double FAz_earth = degAreaOverR2Plus180Minus36 - H - - Az_earth; // F(Az) or g(Az) - double F2Az_earth = - (cosAzEarth * sinGcosSDC2VoS + sinAzEarth * cosG) / - sin(H) - - 1; // F'(Az) or g'(Az) - deltaAz = -FAz_earth / F2Az_earth; // Delta Az^0 or Delta Az - Az_earth += deltaAz; - } - { - double sinAz_earth = sin(Az_earth), - cosAz_earth = cos(Az_earth); - double q = - atan2(tang, (cosAz_earth + sinAz_earth * cotTheta)); - double d = - params->RprimeTang / (cosAz + sinAz * cotTheta); // d' - double f = d / (params->Rprime2X * sin(q / 2)); // f - double z = 2 * asin(rho / (params->Rprime2X * f)); - - Az_earth -= AzAdjustment; - { - const isea_sincos *latSinCos = - ¶ms->vertexLatSinCos[c.face]; - double sinLat0 = latSinCos->s, cosLat0 = latSinCos->c; - double sinZ = sin(z), cosZ = cos(z); - double cosLat0SinZ = cosLat0 * sinZ; - double latSin = - sinLat0 * cosZ + cosLat0SinZ * cos(Az_earth); - double lat = safeArcSin(latSin); - double lon = - facesCenterDodecahedronVertices[c.face].lon + - atan2(sin(Az_earth) * cosLat0SinZ, - cosZ - sinLat0 * sin(lat)); - - revertOrientation({lat, lon}, r); - } - } - } - return true; - } - r = {inf, inf}; - return false; - } - - private: - GeoPoint orientation; - double cosOrientationLat, sinOrientationLat; - - inline void revertOrientation(const GeoPoint &c, GeoPoint &r) { - double lon = (c.lat < DEG_TO_RAD * -90 + precisionPerDefinition || - c.lat > DEG_TO_RAD * 90 - precisionPerDefinition) - ? 0 - : c.lon; - if (orientation.lat != 0.0 || orientation.lon != 0.0) { - double sinLat = sin(c.lat), cosLat = cos(c.lat); - double sinLon = sin(lon), cosLon = cos(lon); - double cosLonCosLat = cosLon * cosLat; - r = {asin(sinLat * cosOrientationLat - - cosLonCosLat * sinOrientationLat), - atan2(sinLon * cosLat, cosLonCosLat * cosOrientationLat + - sinLat * sinOrientationLat) - - orientation.lon}; - } else - r = {c.lat, lon}; - } - - static inline double faceOrientation(int face) { - return (face <= 4 || (10 <= face && face <= 14)) ? 0 : DEG_TO_RAD * 180; - } -}; - -// Orientation symmetric to equator (+proj=isea) -/* Sets the orientation of the icosahedron such that the north and the south - * poles are mapped to the edge midpoints of the icosahedron. The equator is - * thus mapped symmetrically. - */ -static ISEAPlanarProjection standardISEA( - /* DEG_TO_RAD * (90 - 58.282525589) = 31.7174744114613 */ - {(E_RAD + F_RAD) / 2, DEG_TO_RAD * -11.25}); - -// Polar orientation (+proj=isea +orient=pole) -/* - * One corner of the icosahedron is, by default, facing to the north pole, and - * one to the south pole. The provided orientation is relative to the default - * orientation. - * - * The orientation shifts every location by the angle orientation.lon in - * direction of positive longitude, and thereafter by the angle orientation.lat - * in direction of positive latitude. - */ -static ISEAPlanarProjection polarISEA({0, 0}); - -void pj_isea_data::initialize(const PJ *P) { - struct pj_isea_data *Q = static_cast(P->opaque); - // Only supporting default planar options for now - if (Q->output == ISEA_PLANE && Q->o_az == 0.0 && Q->aperture == 3.0 && - Q->resolution == 4.) { - // Only supporting +orient=isea and +orient=pole for now - if (Q->o_lat == ISEA_STD_LAT && Q->o_lon == ISEA_STD_LONG) - p = &standardISEA; - else if (Q->o_lat == M_PI / 2.0 && Q->o_lon == 0) - p = &polarISEA; - else - p = nullptr; - } - - if (p != nullptr) { - if (P->e > 0) { - double a2 = P->a * P->a, c2 = P->b * P->b; - double log1pe_1me = log((1 + P->e) / (1 - P->e)); - double S = M_PI * (2 * a2 + c2 / P->e * log1pe_1me); - R2 = S / (4 * M_PI); // [WGS84] R = 6371007.1809184747 m - Rprime = RprimeOverR * sqrt(R2); // R' - } else { - R2 = P->a * P->a; // R^2 - Rprime = RprimeOverR * P->a; // R' - } - - Rprime2X = 2 * Rprime; - RprimeTang = Rprime * tang; // twice the center-to-base distance - centerToBase = RprimeTang / 2; - triWidth = RprimeTang * SQRT3; - Rprime2Tan2g = RprimeTang * RprimeTang; - - yOffsets[0] = -2 * centerToBase; - yOffsets[1] = -4 * centerToBase; - yOffsets[2] = -5 * centerToBase; - yOffsets[3] = -7 * centerToBase; - - xo = 2.5 * triWidth; - yo = -1.5 * centerToBase; - sx = 1.0 / triWidth; - sy = 1.0 / (3 * centerToBase); - } -} - -} // anonymous namespace - -static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P) { - const struct pj_isea_data *Q = - static_cast(P->opaque); - ISEAPlanarProjection *p = Q->p; - - if (p) { - // Default origin of +proj=isea is different (OGC:1534 is - // +x_0=19186144.870934911 +y_0=-3323137.7717836285) - PJ_XY input{xy.x * P->a + Q->xo, xy.y * P->a + Q->yo}; - GeoPoint result; - - if (p->cartesianToGeo(input, Q, result)) - // TODO: Convert authalic latitude to geodetic latitude if not - // spherical as in eqearth, healpix, laea, etc. - return {result.lon, result.lat}; - else - return {inf, inf}; - } else - return {inf, inf}; -} - -#undef ISEA_STD_LAT -#undef ISEA_STD_LONG - -#undef numIcosahedronFaces -#undef precision -#undef precisionPerDefinition - -#undef AzMax -#undef sdc2vos -#undef tang -#undef cotTheta -#undef cosG -#undef sinGcosSDC2VoS -#undef westVertexLon - -#undef RprimeOverR - -#undef Min -#undef Max - -#undef inf - -#undef E_RAD -#undef F_RAD - -#undef DEG36 -#undef DEG72 -#undef DEG90 -#undef DEG108 -#undef DEG120 -#undef DEG144 -#undef DEG180 -#undef ISEA_SCALE -#undef V_LAT -#undef TABLE_G -#undef TABLE_H diff --git a/src/projections/polyhedral.cpp b/src/projections/polyhedral.cpp new file mode 100644 index 0000000000..0648b3eb23 --- /dev/null +++ b/src/projections/polyhedral.cpp @@ -0,0 +1,230 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Polyhedral Snyder equal-area map projections. + * Author: Felix Palmer + * + ****************************************************************************** + * Derived from A5 (Apache-2.0). + * https://github.com/felixpalmer/a5 + * + * Copyright (c) A5 contributors + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + ****************************************************************************/ + +#include "polyhedral/conway.h" +#include "polyhedral/nets/dsea/dsea.h" +#include "polyhedral/nets/isea/isea.h" +#include "polyhedral/nets/tsea/tsea.h" +#include "polyhedral/polyhedra/dodecahedron.h" +#include "polyhedral/polyhedra/icosahedron.h" +#include "polyhedral/polyhedra/tetrahedron.h" +#include "polyhedral/state.h" + +#include "proj.h" +#include "proj_internal.h" + +#include +#include +#include + +using polyhedral::pj_polyhedral_data; + +static PJ_XY polyhedral_fwd(PJ_LP lp, PJ *P) { + PJ_XY xy = {0.0, 0.0}; + auto *Q = static_cast(P->opaque); + + Vec3 v = polyhedral::lonlat_to_cart(Q, P, lp.lam, lp.phi); + + int tri_idx = polyhedral::find_triangle(Q, v); + if (tri_idx < 0) { + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); + return proj_coord_error().xy; + } + + polyhedral::Face2D result = + polyhedral::snyder_fwd(v, Q->sph_tris[tri_idx], Q->face_tris[tri_idx]); + + xy.x = result.x - Q->x_offset; + xy.y = result.y - Q->y_offset; + + return xy; +} + +static PJ_LP polyhedral_inv(PJ_XY xy, PJ *P) { + PJ_LP lp = {0.0, 0.0}; + auto *Q = static_cast(P->opaque); + + const double x = xy.x + Q->x_offset; + const double y = xy.y + Q->y_offset; + + for (int i = 0; i < Q->n_triangles; i++) { + polyhedral::Barycentric bary = + polyhedral::face_to_barycentric({x, y}, Q->face_tris[i]); + if (bary.u >= -1e-10 && bary.v >= -1e-10 && bary.w >= -1e-10) { + Vec3 v = + polyhedral::snyder_inv({x, y}, Q->face_tris[i], Q->sph_tris[i]); + polyhedral::cart_to_lonlat(Q, P, v, lp.lam, lp.phi); + return lp; + } + } + + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); + return proj_coord_error().lp; +} + +static PJ *polyhedral_destructor(PJ *P, int errlev) { + if (nullptr == P) + return nullptr; + if (nullptr == P->opaque) + return pj_default_destructor(P, errlev); + free(static_cast(P->opaque)->apa); + return pj_default_destructor(P, errlev); +} + +static PJ *polyhedral_setup(PJ *P, const polyhedral::PolyhedralDefaults &d) { + auto *Q = static_cast(P->opaque); + if (P->es != 0.0) { + Q->apa = pj_authalic_lat_compute_coeffs(P->n); + if (nullptr == Q->apa) + return polyhedral_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); + Q->qp = pj_authalic_lat_q(1.0, P); + } + if (!polyhedral::apply_polyhedral_params(Q, P, d)) + return polyhedral_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); + P->fwd = polyhedral_fwd; + P->inv = polyhedral_inv; + return P; +} + +PROJ_HEAD(tsea, "Tetrahedral Snyder Equal Area") +"\n\tSph&Ell\n\tnet= orient_lat= orient_lon= azi= lat_0= lon_0="; +PJ *PJ_PROJECTION(tsea) { + auto *Q = static_cast( + calloc(1, sizeof(pj_polyhedral_data))); + if (nullptr == Q) + return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); + P->opaque = Q; + P->destructor = polyhedral_destructor; + + const char *net_name = pj_param(P->ctx, P->params, "snet").s; + const int(*parents)[4] = &nets::tsea::tsea; + if (net_name != nullptr) { + if (strcmp(net_name, "tsea") == 0) { + parents = &nets::tsea::tsea; + } else if (strcmp(net_name, "star") == 0) { + parents = &nets::tsea::star; + } else { + proj_log_error(P, _("invalid +net (expected tsea or star)")); + return polyhedral_destructor(P, + PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); + } + } + + const polyhedral::PolyhedralDefaults d = {90.0, 0.0, 0.0}; + polyhedral::load_meshes(Q, polyhedra::tetrahedron, *parents, d); + return polyhedral_setup(P, d); +} + +PROJ_HEAD(dsea, "Dodecahedral Snyder Equal Area") +"\n\tSph&Ell\n\tnet= orient_lat= orient_lon= azi= lat_0= lon_0="; +PJ *PJ_PROJECTION(dsea) { + auto *Q = static_cast( + calloc(1, sizeof(pj_polyhedral_data))); + if (nullptr == Q) + return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); + P->opaque = Q; + P->destructor = polyhedral_destructor; + + // +net= selects the unfolding tree. The `a5` net shifts orient_lon by + // -93° to match the LONGITUDE_OFFSET used by the A5 reference codebase. + const char *net_name = pj_param(P->ctx, P->params, "snet").s; + const int(*parents)[12] = &nets::dsea::dsea; + double default_orient_lon = -36.0; + double default_azi = 240.0; + if (net_name != nullptr) { + if (strcmp(net_name, "dsea") == 0) { + parents = &nets::dsea::dsea; + } else if (strcmp(net_name, "a5") == 0) { + parents = &nets::dsea::a5; + default_orient_lon = -36.0 - 93.0; + } else if (strcmp(net_name, "crescent") == 0) { + parents = &nets::dsea::crescent; + } else if (strcmp(net_name, "flower") == 0) { + parents = &nets::dsea::flower; + } else { + proj_log_error(P, _("invalid +net (expected dsea, a5, crescent, " + "or flower)")); + return polyhedral_destructor(P, + PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); + } + } + + // Defaults rotate the polyhedron so F1's centroid sits at the geographic + // north pole. orient_lat = atan((1+2A)/2) places V0 on the F1 axis. + using polyhedra::dodecahedron_constants::A; + const double default_orient_lat = + std::atan((1 + 2 * A) / 2.0) * 180.0 / M_PI; + const polyhedral::PolyhedralDefaults d = {default_orient_lat, + default_orient_lon, default_azi}; + polyhedral::load_meshes(Q, polyhedra::dodecahedron, *parents, d); + return polyhedral_setup(P, d); +} + +PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") +"\n\tSph&Ell\n\torient= orient_lat= orient_lon= azi= lat_0= lon_0="; +PJ *PJ_PROJECTION(isea) { + auto *Q = static_cast( + calloc(1, sizeof(pj_polyhedral_data))); + if (nullptr == Q) + return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); + P->opaque = Q; + P->destructor = polyhedral_destructor; + + // Standard Snyder ISEA orientation: V0 at authalic latitude arctan(φ) ≈ + // 58.28°N (φ = golden ratio). On WGS84 this corresponds to ~58.40° + // geodetic — PolyhedralDefaults stores the authalic value; user + // +orient_lat is geodetic and converted internally. orient_lon = 11.25° + // on a sphere and 11.20° on an ellipsoid matches Snyder's Figure 12. + const bool sphere = (P->es == 0.0); + double orient_lat = 58.282525588539; + double orient_lon = sphere ? 11.25 : 11.20; + double azi = 0.0; + + // +orient= shorthand for the two named orientations from Snyder's paper. + // Individual +orient_lat / +orient_lon / +azi still override. + const char *orient_name = pj_param(P->ctx, P->params, "sorient").s; + if (orient_name != nullptr) { + if (strcmp(orient_name, "isea") == 0) { + // Defaults above already match. + } else if (strcmp(orient_name, "pole") == 0) { + orient_lat = 90.0; + orient_lon = 0.0; + azi = 0.0; + } else { + proj_log_error(P, _("invalid +orient (expected isea or pole)")); + return polyhedral_destructor(P, + PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); + } + } + + // Default projected origin is the bbox center (mid-width, mid-height) of + // the unfold's root face — the "South Africa" upper-band face F8 — rather + // than its centroid. This places the geometric centre of the 4×5 strip + // layout at (0, 0). Other polyhedra default to the centroid of the root + // face. + const polyhedral::PolyhedralDefaults d = { + orient_lat, orient_lon, azi, polyhedral::DefaultOrigin::FaceBboxCenter}; + polyhedral::load_meshes(Q, polyhedra::icosahedron, nets::isea::isea, d); + return polyhedral_setup(P, d); +} diff --git a/src/projections/polyhedral/conway.h b/src/projections/polyhedral/conway.h new file mode 100644 index 0000000000..7b941e4537 --- /dev/null +++ b/src/projections/polyhedral/conway.h @@ -0,0 +1,63 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Conway meta polyhedron operation + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRAL_CONWAY_H +#define POLYHEDRAL_CONWAY_H + +#include "vec3.h" + +namespace polyhedral { + +// Mesh struct, both for representing polyhedra & flat nets +template struct Mesh { + // List of vertices as 3D vectors + Vec3 vertices[NumVertices]; + // Faces, indexing above vertices + int faces[NumFaces][NumFaceVertices]; +}; + +// ConwayMode: +// - Sphere: project centroid + fan vertices onto unit sphere (for polyhedra) +// - Plane: scale centroid to the arithmetic mean (for nets) +enum class ConwayMode { Sphere, Plane }; + +template +inline void conway_meta(const Vec3 (&vertices)[NumVertices], + const int (&faces)[NumFaces][NumFaceVertices], + Vec3 (&out)[2 * NumFaceVertices * NumFaces][3]) { + constexpr int FanSize = 2 * NumFaceVertices; + for (int i = 0; i < NumFaces; i++) { + Vec3 fan[FanSize]; + Vec3 center{0.0, 0.0, 0.0}; + for (int k = 0; k < NumFaceVertices; k++) { + const Vec3 va = vertices[faces[i][k]]; + const Vec3 vb = vertices[faces[i][(k + 1) % NumFaceVertices]]; + fan[2 * k] = va; + fan[2 * k + 1] = vec3_lerp(va, vb, 0.5); + center = vec3_add(center, va); + } + + if constexpr (Mode == ConwayMode::Sphere) { + for (int k = 0; k < FanSize; k++) + fan[k] = vec3_normalize(fan[k]); + center = vec3_normalize(center); + } else { + center = vec3_scale(center, 1.0 / NumFaceVertices); + } + + for (int k = 0; k < FanSize; k++) { + out[FanSize * i + k][0] = center; + out[FanSize * i + k][1] = fan[(k + 1) % FanSize]; + out[FanSize * i + k][2] = fan[k]; + } + } +} + +} // namespace polyhedral + +#endif // POLYHEDRAL_CONWAY_H diff --git a/src/projections/polyhedral/nets/dsea/dsea.h b/src/projections/polyhedral/nets/dsea/dsea.h new file mode 100644 index 0000000000..94e088ed3d --- /dev/null +++ b/src/projections/polyhedral/nets/dsea/dsea.h @@ -0,0 +1,61 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Net layouts for the DSEA projection. Each constant is a parents + * tree over the 12 faces of the dodecahedron: parents[i] stores + * the 1-based ID of face F(i+1)'s parent in the unfolding tree, + * or 0 if F(i+1) is the root. Face IDs are 1-based to match the + * Snyder paper; 0 means "no parent". Every parent–child pair + * must share an edge on the dodecahedron. + * + * The dodecahedron face order is top-to-bottom, then left-to-right + * within each band: + * F1 north polar cap + * F2..F6 upper band N.Pacific, N.America, Europe, + * India, Asia + * F7..F11 lower band S.Pacific, S.Atlantic, S.Africa, + * Indian, S.Pacific E + * F12 south polar cap + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef NETS_DSEA_DSEA_H +#define NETS_DSEA_DSEA_H + +namespace nets { +namespace dsea { + +// Snyder's two-rows-of-5 layout (Figure 6): F9 (Africa) is the central face +// of the lower row, with the row extending left and right via zig-zag +// connections through the upper row — directly hinging adjacent lower-row +// faces would curl the band into a crescent, so each next-lower-row face is +// reached through its upper-row neighbour instead. +// Right: F9 → F5 → F10 → F6 → F11 (L → U → L → U → L) +// Left: F9 → F4 → F8 → F3 → F7 → F2 (L → U → L → U → L → U) +// The two polar caps hang symmetrically above and below the central column: +// F1 (north pole) off F4 (Europe), F12 (south pole) off F9 (Africa). +constexpr int dsea[12] = {4, 7, 8, 0, 9, 10, 3, 4, 4, 5, 6, 9}; + +// A5 net (re-rooted to match the A5 reference codebase's coordinate origin). +// The visual shape is the same as the published A5 layout; only the parent +// pointers along the path to the new root are reversed. +constexpr int a5[12] = {4, 6, 8, 9, 6, 0, 2, 12, 10, 5, 7, 11}; + +// Crescent: central spine F1 → F4 → F9 → F12 (north pole → UK/Europe → +// Africa → south pole). Two symmetric "crescent" arcs extend out from F4 +// and F9: one through the upper row (F4 → F5 → F6 → F2 → F3) and one +// through the lower row (F9 → F10 → F11 → F7 → F8). +constexpr int crescent[12] = {4, 3, 4, 0, 4, 5, 8, 9, 4, 9, 10, 9}; + +// Flower: F4 (Europe) at the centre, with F1 (north pole) as one petal and +// F9 (Africa) as the link down to the south pole. F1's four other petals +// (F2, F3, F5, F6) hang off F1; F12's four other petals (F7, F8, F10, F11) +// hang off F12. This routes the polar caps through Africa rather than +// through the Americas. +constexpr int flower[12] = {4, 1, 1, 0, 1, 1, 12, 12, 4, 12, 12, 9}; + +} // namespace dsea +} // namespace nets + +#endif // NETS_DSEA_DSEA_H diff --git a/src/projections/polyhedral/nets/isea/isea.h b/src/projections/polyhedral/nets/isea/isea.h new file mode 100644 index 0000000000..9b2a7aa16e --- /dev/null +++ b/src/projections/polyhedral/nets/isea/isea.h @@ -0,0 +1,52 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Net layout for the ISEA projection — the regular icosahedron + * unfolded into the canonical Snyder 4×5 strip. + * + * Face IDs are 1-based (matching Snyder); 0 means "no parent" + * (the unfold root). + * + * The icosahedron face order in `polyhedra::icosahedron` is + * grouped into 4 rows of 5: + * F1..F5 north cap (around V0) + * F6..F10 upper band (down-pointing on the sphere) + * F11..F15 lower band (up-pointing on the sphere) + * F16..F20 south cap (around V11) + * + * F8 is the centred upper-band face — the "South Africa" face in + * Snyder's Figure 12 — and is chosen as the root of the unfolding + * tree so that the default projected origin is its centroid / + * bbox centre, matching the geometric centre of the legacy ISEA + * strip. + * + * The unfold tree zig-zags from the root F8 outward across the + * upper and lower bands, + * + * F1 ← F6 ← F11 ← F7 ← F12 ← F8 → F13 → F9 → F14 → F10 → F15 + * + * with each remaining north-cap face hung off its upper-band + * neighbour and each south-cap face hung off its lower-band + * neighbour. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef NETS_ISEA_ISEA_H +#define NETS_ISEA_ISEA_H + +namespace nets { +namespace isea { + +constexpr int isea[20] = { + 6, 7, 8, 9, 10, // F1..F5 north cap, hung off upper band + 11, 12, 0, 13, 14, // F6, F7 climb back to F1; F8 is root; F9, F10 chain on + 7, 8, 8, 9, 10, // F11..F15 lower band — left of F8 climbs back, right + // chains forward + 11, 12, 13, 14, 15, // F16..F20 south cap, hung off lower band +}; + +} // namespace isea +} // namespace nets + +#endif // NETS_ISEA_ISEA_H diff --git a/src/projections/polyhedral/nets/tsea/tsea.h b/src/projections/polyhedral/nets/tsea/tsea.h new file mode 100644 index 0000000000..b8c5f0761a --- /dev/null +++ b/src/projections/polyhedral/nets/tsea/tsea.h @@ -0,0 +1,34 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Net layouts for the TSEA projection — the regular tetrahedron + * unfolded into a 2D net. + * + * Face IDs are 1-based (matching Snyder); 0 means "no parent" + * (the unfold root). + * + * The default `tsea` net is Snyder's Figure 8 horizontal strip + * F1 — F2 — F3 — F4 (North America — Europe/Africa — south cap — + * Asia/Australia), rerooted at F3 so the projected (0, 0) lands + * on the south-cap face (Antarctica). + * + * The `star` net is the alternative layout — F3 (south cap) at + * the centre with its three V0-touching neighbours (F1, F2, F4) + * hung off each of its edges in a three-way fan. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef NETS_TSEA_TSEA_H +#define NETS_TSEA_TSEA_H + +namespace nets { +namespace tsea { + +constexpr int tsea[4] = {2, 3, 0, 3}; +constexpr int star[4] = {3, 3, 0, 3}; + +} // namespace tsea +} // namespace nets + +#endif // NETS_TSEA_TSEA_H diff --git a/src/projections/polyhedral/polyhedra/dodecahedron.h b/src/projections/polyhedral/polyhedra/dodecahedron.h new file mode 100644 index 0000000000..8ccd4e93fb --- /dev/null +++ b/src/projections/polyhedral/polyhedra/dodecahedron.h @@ -0,0 +1,90 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Regular dodecahedron — 20 vertices, 12 pentagonal faces, + * inscribed in the unit sphere with V0 at (0, 0, 1) and V1 along + * the +x meridian. The 20 vertices form six layers along z: + * + * z = 1 V0 (1 vertex) + * z = √5/3 = (4A−1)/3 V1, V4, V5 (layer 1, 3 vertices) + * z = 1/3 6 vertices (layer 2) + * z = −1/3 6 vertices (layer 3) + * z = −√5/3 3 vertices (layer 4) + * z = −1 V17 (1 vertex) + * + * where A = (1+√5)/4 = cos(36°) and B = √3/3 = 1/√3 are the only + * algebraic constants needed (every other coordinate is a rational + * combination of A, B and 1/3). + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRA_DODECAHEDRON_H +#define POLYHEDRA_DODECAHEDRON_H + +#include "../conway.h" +#include "validation.h" + +namespace polyhedra { +namespace dodecahedron_constants { + +constexpr double A = 0.809016994374947424102293417183; // (1+√5)/4 = cos(36°) +constexpr double B = 0.577350269189625764509148780502; // √3/3 = 1/√3 + +constexpr polyhedral::Mesh<20, 12, 5> mesh = { + // Vertices (unit sphere; V0 at +z, V1 in +x meridian) + { + {0.0, 0.0, 1.0}, + {2.0 / 3.0, 0.0, (4 * A - 1) / 3.0}, + {(4 * A - 1) / 3.0, B, 1.0 / 3.0}, + {(2 - 2 * A) / 3.0, 2 * A *B, 1.0 / 3.0}, + {-1.0 / 3.0, B, (4 * A - 1) / 3.0}, + {-1.0 / 3.0, -B, (4 * A - 1) / 3.0}, + {(4 * A - 1) / 3.0, -B, 1.0 / 3.0}, + {(1 + 2 * A) / 3.0, (2 * A - 1) * B, -1.0 / 3.0}, + {(2 * A - 2) / 3.0, 2 * A *B, -1.0 / 3.0}, + {-(1 + 2 * A) / 3.0, (2 * A - 1) * B, 1.0 / 3.0}, + {(2 - 2 * A) / 3.0, -2 * A *B, 1.0 / 3.0}, + {(1 + 2 * A) / 3.0, -(2 * A - 1) * B, -1.0 / 3.0}, + {1.0 / 3.0, B, -(4 * A - 1) / 3.0}, + {-(4 * A - 1) / 3.0, B, -1.0 / 3.0}, + {-(1 + 2 * A) / 3.0, -(2 * A - 1) * B, 1.0 / 3.0}, + {(2 * A - 2) / 3.0, -2 * A *B, -1.0 / 3.0}, + {1.0 / 3.0, -B, -(4 * A - 1) / 3.0}, + {0.0, 0.0, -1.0}, + {-2.0 / 3.0, 0.0, -(4 * A - 1) / 3.0}, + {-(4 * A - 1) / 3.0, -B, -1.0 / 3.0}, + }, + // Faces (CCW from outside). Ordered top-to-bottom, then left-to-right + // within each band, so the unfolded DSEA net comes out as: + // F1 (north polar cap, around V0) + // F2 F3 F4 F5 F6 (upper band: N.Pacific, N.America, + // Europe, India, Asia) + // F7 F8 F9 F10 F11 (lower band: S.Pacific, S.Atlantic, + // S.Africa, Indian, S.Pacific E) + // F12 (south polar cap, around V17) + { + {0, 1, 2, 3, 4}, + {3, 8, 13, 9, 4}, + {0, 4, 9, 14, 5}, + {0, 5, 10, 6, 1}, + {1, 6, 11, 7, 2}, + {2, 7, 12, 8, 3}, + {9, 13, 18, 19, 14}, + {5, 14, 19, 15, 10}, + {6, 10, 15, 16, 11}, + {7, 11, 16, 17, 12}, + {8, 12, 17, 18, 13}, + {15, 19, 18, 17, 16}, + }, +}; + +[[maybe_unused]] constexpr validate_polyhedron _validate{}; + +} // namespace dodecahedron_constants + +constexpr auto &dodecahedron = dodecahedron_constants::mesh; + +} // namespace polyhedra + +#endif // POLYHEDRA_DODECAHEDRON_H diff --git a/src/projections/polyhedral/polyhedra/icosahedron.h b/src/projections/polyhedral/polyhedra/icosahedron.h new file mode 100644 index 0000000000..0e687ab4e3 --- /dev/null +++ b/src/projections/polyhedral/polyhedra/icosahedron.h @@ -0,0 +1,70 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Regular icosahedron inscribed in the unit sphere, oriented with + * V0 at the north pole. The remaining 11 vertices form two rings + * of 5 (at z = ±1/√5) staggered by 36° in longitude, plus V11 at + * the south pole. Faces are listed in 4 groups of 5 — north cap, + * upper band (down-pointing on the sphere), lower band (up-pointing + * on the sphere), south cap — chosen to match the canonical ISEA + * net unfolding. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRA_ICOSAHEDRON_H +#define POLYHEDRA_ICOSAHEDRON_H + +#include "../conway.h" +#include "validation.h" + +namespace polyhedra { +namespace icosahedron_constants { + +constexpr double A = 0.894427190999915878564; // 2/√5 (ring xy radius) +constexpr double B = 0.447213595499957939282; // 1/√5 (ring |z|) + +// cos/sin of multiples of 72° and 36° (used for the two rings). +constexpr double C72 = 0.309016994374947424102; // cos 72° +constexpr double S72 = 0.951056516295153572116; // sin 72° +constexpr double C36 = 0.809016994374947424102; // cos 36° +constexpr double S36 = 0.587785252292473129169; // sin 36° + +constexpr polyhedral::Mesh<12, 20, 3> mesh = { + // Vertices (unit sphere; V0 at +z, V1 in +x meridian) + { + {0.0, 0.0, 1.0}, + {-A, 0.0, B}, + {-A * C72, -A *S72, B}, + {A * C36, -A *S36, B}, + {A * C36, A *S36, B}, + {-A * C72, A *S72, B}, + {-A * C36, -A *S36, -B}, + {A * C72, -A *S72, -B}, + {A, 0.0, -B}, + {A * C72, A *S72, -B}, + {-A * C36, A *S36, -B}, + {0.0, 0.0, -1.0}, + }, + // Faces (CCW from outside). Ordered north to south, in groups of 5 + // F1..F5 north cap (around V0) + // F6..F10 upper band (down-pointing) + // F11..F15 lower band (up-pointing) + // F16..F20 south cap (around V11) + { + {0, 1, 2}, {0, 2, 3}, {0, 3, 4}, {0, 4, 5}, {0, 5, 1}, + {1, 6, 2}, {2, 7, 3}, {3, 8, 4}, {4, 9, 5}, {5, 10, 1}, + {6, 7, 2}, {7, 8, 3}, {8, 9, 4}, {9, 10, 5}, {10, 6, 1}, + {11, 7, 6}, {11, 8, 7}, {11, 9, 8}, {11, 10, 9}, {11, 6, 10}, + }, +}; + +[[maybe_unused]] constexpr validate_polyhedron _validate{}; + +} // namespace icosahedron_constants + +constexpr auto &icosahedron = icosahedron_constants::mesh; + +} // namespace polyhedra + +#endif // POLYHEDRA_ICOSAHEDRON_H diff --git a/src/projections/polyhedral/polyhedra/tetrahedron.h b/src/projections/polyhedral/polyhedra/tetrahedron.h new file mode 100644 index 0000000000..442a27e036 --- /dev/null +++ b/src/projections/polyhedral/polyhedra/tetrahedron.h @@ -0,0 +1,62 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Regular tetrahedron inscribed in the unit sphere, oriented with + * V0 at the north pole. The remaining vertices sit at latitude + * -arcsin(1/3) ≈ -19.47°, separated by 120° in longitude + * (V1 at +60°, V2 at -60°, V3 at 180°). The three V0-touching + * faces are listed west-to-east, with the south-cap face slotted + * in between (so the default unfold strip F1 — F2 — F3 — F4 reads + * North America — Europe/Africa — south cap — Asia/Australia): + * F1 = V0 + V3 + V2 covers North America (and the Pacific) + * F2 = V0 + V2 + V1 covers Europe + Africa (and Atlantic) + * F3 = V1 + V2 + V3 south-cap face opposite V0 + * F4 = V0 + V1 + V3 covers Asia + Australia (and W Pacific) + * The F1–F2 edge V0-V2 lies on the lon=-60° meridian, splitting + * South America between the two faces (matching Snyder Figure 8). + * F3's vertices are listed V1, V2, V3 so the unfold's polar-face + * fallback (which places face[0] along +y) puts V1 — the +60° + * meridian — "up" in the Antarctica-centred tsea net. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRA_TETRAHEDRON_H +#define POLYHEDRA_TETRAHEDRON_H + +#include "../conway.h" +#include "validation.h" + +namespace polyhedra { +namespace tetrahedron_constants { + +constexpr double A = 0.9428090415820635; // sqrt(8/9) +constexpr double B = 0.8164965809277260; // sqrt(2/3) + +constexpr polyhedral::Mesh<4, 4, 3> mesh = { + // Vertices (unit sphere; V0 at +z, V1 at lon=+60°) + { + {0.0, 0.0, 1.0}, + {0.5 * A, B, -1.0 / 3.0}, + {0.5 * A, -B, -1.0 / 3.0}, + {-A, 0.0, -1.0 / 3.0}, + }, + // Faces (CCW from outside). F1, F2, F4 share V0; F3 is the south-cap + // face opposite V0, listed in V1, V2, V3 order so the unfold places V1 up. + { + {0, 3, 2}, + {0, 2, 1}, + {1, 2, 3}, + {0, 1, 3}, + }, +}; + +[[maybe_unused]] constexpr validate_polyhedron _validate{}; + +} // namespace tetrahedron_constants + +constexpr auto &tetrahedron = tetrahedron_constants::mesh; + +} // namespace polyhedra + +#endif // POLYHEDRA_TETRAHEDRON_H diff --git a/src/projections/polyhedral/polyhedra/validation.h b/src/projections/polyhedral/polyhedra/validation.h new file mode 100644 index 0000000000..9cd5d48f4a --- /dev/null +++ b/src/projections/polyhedral/polyhedra/validation.h @@ -0,0 +1,68 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Compile-time validation helper for polyhedron seed meshes. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRA_VALIDATION_H +#define POLYHEDRA_VALIDATION_H + +#include "../conway.h" + +namespace polyhedra { + +template +constexpr bool +all_vertices_on_unit_sphere(const polyhedral::Mesh &m) { + constexpr double tol = 1e-12; + for (int i = 0; i < NV; i++) { + const auto &v = m.vertices[i]; + const double err = v.x * v.x + v.y * v.y + v.z * v.z - 1.0; + if (err < -tol || err > tol) + return false; + } + return true; +} + +template +constexpr bool +face_indices_well_formed(const polyhedral::Mesh &m) { + for (int f = 0; f < NF; f++) { + for (int k = 0; k < NFV; k++) { + if (m.faces[f][k] < 0 || m.faces[f][k] >= NV) + return false; + for (int j = k + 1; j < NFV; j++) { + if (m.faces[f][k] == m.faces[f][j]) + return false; + } + } + } + return true; +} + +// Two topological conditions: +// V - E + F = 2 for any closed convex polyhedron, +// E = NF * NFV / 2 (each edge shared by exactly two NFV-gon faces). +template +constexpr bool euler_holds(const polyhedral::Mesh &) { + return (NF * NFV) % 2 == 0 && NV - (NF * NFV) / 2 + NF == 2; +} + +template struct validate_polyhedron { + static_assert(M.vertices[0].x == 0.0 && M.vertices[0].y == 0.0 && + M.vertices[0].z == 1.0, + "V0 must be at (0, 0, 1)"); + static_assert(all_vertices_on_unit_sphere(M), + "all vertices must lie on the unit sphere"); + static_assert(face_indices_well_formed(M), + "face vertex indices must be in [0, NV) " + "and unique within each face"); + static_assert(euler_holds(M), + "polyhedron must satisfy Euler's formula V - E + F = 2"); +}; + +} // namespace polyhedra + +#endif // POLYHEDRA_VALIDATION_H diff --git a/src/projections/polyhedral/snyder.h b/src/projections/polyhedral/snyder.h new file mode 100644 index 0000000000..d29cb19ac3 --- /dev/null +++ b/src/projections/polyhedral/snyder.h @@ -0,0 +1,167 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Snyder equal-area polyhedral projection — forward and inverse. + * Author: Felix Palmer + * + ****************************************************************************** + * Derived from A5 (Apache-2.0). + * https://github.com/felixpalmer/a5 + * (modules/projections/polyhedral.ts) + * + * Copyright (c) A5 contributors + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + **************************************************************************** + * Derived from DGGAL (BSD 3-Clause License) + * + * Copyright (c) 2014-2025, Ecere Corporation + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, + * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN + * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + **************************************************************************** + * + * Closed-form equations adapted from: + * https://brsr.github.io/2021/08/31/snyder-equal-area.html + * + ****************************************************************************/ + +#ifndef SNYDER_H +#define SNYDER_H + +#include "sphere.h" +#include "vec3.h" + +namespace polyhedral { + +struct Face2D { + double x, y; +}; + +struct Barycentric { + double u, v, w; +}; + +struct FaceTriangle { + Face2D a, b, c; +}; + +struct SphericalTriangle { + Vec3 a, b, c; +}; + +inline Barycentric face_to_barycentric(Face2D p, FaceTriangle tri) { + double d31[2] = {tri.a.x - tri.c.x, tri.a.y - tri.c.y}; + double d23[2] = {tri.c.x - tri.b.x, tri.c.y - tri.b.y}; + double d3p[2] = {p.x - tri.c.x, p.y - tri.c.y}; + + double det = d23[0] * d31[1] - d23[1] * d31[0]; + double b0 = (d23[0] * d3p[1] - d23[1] * d3p[0]) / det; + double b1 = (d31[0] * d3p[1] - d31[1] * d3p[0]) / det; + double b2 = 1.0 - (b0 + b1); + return {b0, b1, b2}; +} + +inline Face2D barycentric_to_face(Barycentric b, FaceTriangle tri) { + return {b.u * tri.a.x + b.v * tri.b.x + b.w * tri.c.x, + b.u * tri.a.y + b.v * tri.b.y + b.w * tri.c.y}; +} + +inline Face2D snyder_fwd(const Vec3 &v, const SphericalTriangle &stri, + const FaceTriangle &ftri) { + Vec3 a = stri.a, b = stri.b, c = stri.c; + + Vec3 z = vec3_normalize(vec3_subtract(v, a)); + Vec3 p = vec3_normalize(quadruple_product(a, z, b, c)); + + // Ensure p is on the same hemisphere as b and c (the quadruple product + // may return the antipodal intersection depending on vertex winding) + if (vec3_dot(p, b) + vec3_dot(p, c) < 0) { + p = {-p.x, -p.y, -p.z}; + } + + double h = vector_difference(a, v) / vector_difference(a, p); + double area_abc = spherical_triangle_area(a, b, c); + double scaled_area = h / area_abc; + double area_apc = spherical_triangle_area(a, p, c); + double area_abp = spherical_triangle_area(a, b, p); + + Barycentric bary = {1.0 - h, scaled_area * area_apc, + scaled_area * area_abp}; + return barycentric_to_face(bary, ftri); +} + +inline Vec3 snyder_inv(const Face2D &face_point, const FaceTriangle &ftri, + const SphericalTriangle &stri) { + Vec3 a = stri.a, b = stri.b, c = stri.c; + Barycentric bcoords = face_to_barycentric(face_point, ftri); + + constexpr double threshold = 1.0 - 1e-14; + if (bcoords.u > threshold) + return a; + if (bcoords.v > threshold) + return b; + if (bcoords.w > threshold) + return c; + + Vec3 c1 = vec3_cross(b, c); + double area_abc = spherical_triangle_area(a, b, c); + double h = 1.0 - bcoords.u; + double r = bcoords.w / h; + double alpha = r * area_abc; + double s = std::sin(alpha); + double half_c = std::sin(alpha / 2.0); + double cc = 2.0 * half_c * half_c; // 1 - cos(alpha) + + double c01 = vec3_dot(a, b); + double c12 = vec3_dot(b, c); + double c20 = vec3_dot(c, a); + double s12 = vec3_length(c1); + + double tp = vec3_dot(a, c1); // triple product A·(B×C) + double f = s * tp + cc * (c01 * c12 - c20); + double g = cc * s12 * (1.0 + c01); + double q = (2.0 / std::acos(c12)) * std::atan2(g, f); + Vec3 pp = slerp(b, c, q); + double k = vector_difference(a, pp); + double t = safe_acos(h * k) / safe_acos(k); + return slerp(a, pp, t); +} + +} // namespace polyhedral + +#endif // SNYDER_H diff --git a/src/projections/polyhedral/state.h b/src/projections/polyhedral/state.h new file mode 100644 index 0000000000..6443628347 --- /dev/null +++ b/src/projections/polyhedral/state.h @@ -0,0 +1,295 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Polyhedral projection runtime state — `pj_polyhedral_data` + * struct, mesh loading, orientation, triangle containment, and + * lon/lat ↔ cartesian conversion with orientation applied. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRAL_STATE_H +#define POLYHEDRAL_STATE_H + +#include "conway.h" +#include "snyder.h" +#include "unfold.h" + +#include "proj.h" +#include "proj_internal.h" + +#include +#include + +namespace polyhedral { + +constexpr int MAX_TRIANGLES = 120; + +struct pj_polyhedral_data { + int n_triangles; + int face_vertex_count; // NFV: 3 for triangular faces, 5 for pentagonal + int num_faces; + + SphericalTriangle sph_tris[MAX_TRIANGLES]; + FaceTriangle face_tris[MAX_TRIANGLES]; + + double orient[3][3]; + double orient_inv[3][3]; + + int root_face_index; // unfold-tree root; default reference face for + // lat_0/lon_0 + + // Translation in projected space so the geographic anchor named by + // +lat_0 / +lon_0 (or its dynamic default) lands at (0, 0). + // +x_0/+y_0 stack on top via PROJ core. + double x_offset; + double y_offset; + + // Authalic latitude state (ellipsoidal case only; apa is null on spheres). + double *apa; + double qp; +}; + +// Projected origin location on chosen reference face when the user does not +// supply +lat_0 / +lon_0. +enum class DefaultOrigin { + FaceCentroid, // arithmetic mean of the face's outer vertices in the net + FaceBboxCenter, // axis-aligned bbox center of the face's outer vertices +}; + +// Defaults supplied per-projection; user params override these where present. +// +// orient_lat_deg is in *authalic* latitude — this is the frame the rotation +// operates in (the unit sphere is the authalic sphere for equal-area +// projections). User-supplied +orient_lat is geodetic and gets converted to +// authalic before use. On a sphere the two are identical. +struct PolyhedralDefaults { + double orient_lat_deg; + double orient_lon_deg; + double azi_deg; + DefaultOrigin default_origin = DefaultOrigin::FaceCentroid; +}; + +// Geographic +z (the orient frame's "up") expressed in polyhedron +// coordinates — the third column of the orient rotation matrix. The unfold +// uses this so the root face comes out with true geographic north up. +inline Vec3 orient_up_in_poly(const PolyhedralDefaults &d) { + const double lat = d.orient_lat_deg * DEG_TO_RAD; + const double az = d.azi_deg * DEG_TO_RAD; + return {-std::cos(az) * std::cos(lat), -std::sin(az) * std::cos(lat), + std::sin(lat)}; +} + +// Build triangle data from a polyhedron and the unfold-tree `parents` array. +// The unfold orients the root face so that geographic north (derived from +// the orient defaults) lands on +y in the net. +template +inline void load_meshes(pj_polyhedral_data *Q, + const Mesh &polyhedron, + const int (&parents)[NF], const PolyhedralDefaults &d) { + const auto net = unfold_net(polyhedron, parents, orient_up_in_poly(d)); + constexpr int N = 2 * NFV * NF; + Vec3 sph[N][3]; + Vec3 face[N][3]; + conway_meta(polyhedron.vertices, polyhedron.faces, sph); + conway_meta(net.vertices, net.faces, face); + + Q->n_triangles = N; + Q->face_vertex_count = NFV; + Q->num_faces = NF; + Q->root_face_index = 0; + while (Q->root_face_index < NF && parents[Q->root_face_index] != 0) + Q->root_face_index++; + for (int i = 0; i < N; i++) { + Q->sph_tris[i] = {sph[i][0], sph[i][1], sph[i][2]}; + Q->face_tris[i] = { + {face[i][0].x, face[i][0].y}, + {face[i][1].x, face[i][1].y}, + {face[i][2].x, face[i][2].y}, + }; + } +} + +// Reference point on face_idx in the unfolded planar net (centroid or bbox +// center) +inline Face2D face_reference_point(const pj_polyhedral_data *Q, int face_idx, + DefaultOrigin kind) { + const int nfv = Q->face_vertex_count; + const int fan_size = 2 * nfv; + double sum_x = 0.0, sum_y = 0.0; + double min_x = std::numeric_limits::infinity(); + double max_x = -std::numeric_limits::infinity(); + double min_y = std::numeric_limits::infinity(); + double max_y = -std::numeric_limits::infinity(); + for (int k = 0; k < nfv; k++) { + const Face2D &v = Q->face_tris[fan_size * face_idx + 2 * k].c; + sum_x += v.x; + sum_y += v.y; + if (v.x < min_x) + min_x = v.x; + if (v.x > max_x) + max_x = v.x; + if (v.y < min_y) + min_y = v.y; + if (v.y > max_y) + max_y = v.y; + } + if (kind == DefaultOrigin::FaceBboxCenter) + return {0.5 * (min_x + max_x), 0.5 * (min_y + max_y)}; + return {sum_x / nfv, sum_y / nfv}; +} + +// Build rotation matrix that takes the (lat, lon) point on the authalic unit +// sphere to the north pole, followed by an azimuthal rotation about the +// z-axis. All angles are in degrees; lat is authalic latitude (callers must +// convert from geodetic if needed). +inline void set_orient_from_angles(pj_polyhedral_data *Q, double lat_deg, + double lon_deg, double az_deg) { + const double lat = lat_deg * DEG_TO_RAD; + const double lon = lon_deg * DEG_TO_RAD; + const double az = az_deg * DEG_TO_RAD; + const double alpha = M_HALFPI - lat; + const double ca = std::cos(alpha), sa = std::sin(alpha); + const double cl = std::cos(lon), sl = std::sin(lon); + const double cz = std::cos(az), sz = std::sin(az); + + const double r00 = ca * cl, r01 = ca * sl, r02 = -sa; + const double r10 = -sl, r11 = cl, r12 = 0.0; + const double r20 = sa * cl, r21 = sa * sl, r22 = ca; + + double m[3][3]; + m[0][0] = cz * r00 - sz * r10; + m[0][1] = cz * r01 - sz * r11; + m[0][2] = cz * r02 - sz * r12; + m[1][0] = sz * r00 + cz * r10; + m[1][1] = sz * r01 + cz * r11; + m[1][2] = sz * r02 + cz * r12; + m[2][0] = r20; + m[2][1] = r21; + m[2][2] = r22; + + for (int r = 0; r < 3; r++) + for (int c = 0; c < 3; c++) { + Q->orient[r][c] = m[r][c]; + Q->orient_inv[r][c] = m[c][r]; + } +} + +// Find the triangle containing point v by linear search. +inline int find_triangle(const pj_polyhedral_data *Q, const Vec3 &v) { + for (int i = 0; i < Q->n_triangles; i++) { + if (vec3_dot(v, Q->sph_tris[i].a) <= 0) + continue; + double t1 = triple_product(v, Q->sph_tris[i].a, Q->sph_tris[i].b); + double t2 = triple_product(v, Q->sph_tris[i].b, Q->sph_tris[i].c); + double t3 = triple_product(v, Q->sph_tris[i].c, Q->sph_tris[i].a); + if ((t1 >= 0 && t2 >= 0 && t3 >= 0) || (t1 <= 0 && t2 <= 0 && t3 <= 0)) + return i; + } + return -1; +} + +// Convert geodetic lon/lat to a unit vector and apply the orientation +// rotation. On an ellipsoid, latitude is first mapped to the authalic +// latitude so that equal-area is preserved on the sphere. +inline Vec3 lonlat_to_cart(const pj_polyhedral_data *Q, const PJ *P, double lam, + double phi) { + double auth_lat = phi; + if (P->es != 0.0) { + auth_lat = pj_authalic_lat(phi, std::sin(phi), std::cos(phi), Q->apa, P, + Q->qp); + } + double ca = std::cos(auth_lat), sa = std::sin(auth_lat); + double cl = std::cos(lam), sl = std::sin(lam); + double gx = ca * cl, gy = ca * sl, gz = sa; + return {Q->orient[0][0] * gx + Q->orient[0][1] * gy + Q->orient[0][2] * gz, + Q->orient[1][0] * gx + Q->orient[1][1] * gy + Q->orient[1][2] * gz, + Q->orient[2][0] * gx + Q->orient[2][1] * gy + Q->orient[2][2] * gz}; +} + +// Convert unit vector back to lon/lat, applying inverse orientation. +inline void cart_to_lonlat(const pj_polyhedral_data *Q, const PJ *P, + const Vec3 &v, double &lam, double &phi) { + double gx = Q->orient_inv[0][0] * v.x + Q->orient_inv[0][1] * v.y + + Q->orient_inv[0][2] * v.z; + double gy = Q->orient_inv[1][0] * v.x + Q->orient_inv[1][1] * v.y + + Q->orient_inv[1][2] * v.z; + double gz = Q->orient_inv[2][0] * v.x + Q->orient_inv[2][1] * v.y + + Q->orient_inv[2][2] * v.z; + + double auth_lat = aasin(P->ctx, gz); + phi = (P->es != 0.0) ? pj_authalic_lat_inverse(auth_lat, Q->apa, P, Q->qp) + : auth_lat; + lam = adjlon(aatan2(gy, gx)); +} + +// Read +orient_lat / +orient_lon / +azi / +lat_0 / +lon_0 (falling back to +// derived defaults) and build the orientation matrix plus the projected-space +// offset so the chosen geographic anchor lands at (0, 0). +inline bool apply_polyhedral_params(pj_polyhedral_data *Q, PJ *P, + const PolyhedralDefaults &d) { + // Default is authalic; user input is geodetic and converted below. + double orient_lat = d.orient_lat_deg; + double orient_lon = d.orient_lon_deg; + double azi = d.azi_deg; + if (pj_param(P->ctx, P->params, "torient_lat").i) { + const double user_lat_deg = + pj_param(P->ctx, P->params, "dorient_lat").f; + if (P->es != 0.0) { + const double phi = user_lat_deg * DEG_TO_RAD; + const double auth = pj_authalic_lat( + phi, std::sin(phi), std::cos(phi), Q->apa, P, Q->qp); + orient_lat = auth / DEG_TO_RAD; + } else { + orient_lat = user_lat_deg; + } + } + if (pj_param(P->ctx, P->params, "torient_lon").i) + orient_lon = pj_param(P->ctx, P->params, "dorient_lon").f; + if (pj_param(P->ctx, P->params, "tazi").i) + azi = pj_param(P->ctx, P->params, "dazi").f; + + set_orient_from_angles(Q, orient_lat, orient_lon, azi); + + // Bypass PROJ core's automatic +lon_0 subtraction. +lat_0/+lon_0 names a + // geographic point that must be applied through orient, not a plain + // longitude shift on the input. + P->lam0 = 0.0; + P->phi0 = 0.0; + Q->x_offset = 0.0; + Q->y_offset = 0.0; + + const bool user_lat_0 = pj_param(P->ctx, P->params, "tlat_0").i != 0; + const bool user_lon_0 = pj_param(P->ctx, P->params, "tlon_0").i != 0; + + if (user_lat_0 || user_lon_0) { + // User-supplied geographic anchor: project it through the full orient. + const double lat_0 = + user_lat_0 ? pj_param(P->ctx, P->params, "rlat_0").f : 0.0; + const double lon_0 = + user_lon_0 ? pj_param(P->ctx, P->params, "rlon_0").f : 0.0; + Vec3 v = lonlat_to_cart(Q, P, lon_0, lat_0); + int tri_idx = find_triangle(Q, v); + if (tri_idx < 0) { + proj_log_error(P, + _("+lat_0 / +lon_0 lies outside the polyhedron's " + "coverage; cannot be used as projected origin")); + return false; + } + Face2D r = snyder_fwd(v, Q->sph_tris[tri_idx], Q->face_tris[tri_idx]); + Q->x_offset = r.x; + Q->y_offset = r.y; + return true; + } + + // Default anchor: a planar reference point on the unfold's root face. + const Face2D ref = + face_reference_point(Q, Q->root_face_index, d.default_origin); + Q->x_offset = ref.x; + Q->y_offset = ref.y; + return true; +} + +} // namespace polyhedral + +#endif // POLYHEDRAL_STATE_H diff --git a/src/projections/polyhedral/unfold.h b/src/projections/polyhedral/unfold.h new file mode 100644 index 0000000000..3d4155de5a --- /dev/null +++ b/src/projections/polyhedral/unfold.h @@ -0,0 +1,216 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Unfold a 3D polyhedron Mesh into a planar net Mesh by laying each + * face flat in its own local 2D frame and chaining 2D rigid + * transforms over a parents-tree of faces. + * Author: Felix Palmer + * + ****************************************************************************/ + +#ifndef POLYHEDRAL_UNFOLD_H +#define POLYHEDRAL_UNFOLD_H + +#include "conway.h" +#include "vec3.h" + +#include "proj_internal.h" + +#include +#include + +namespace polyhedral { + +// Number of net vertices produced by unfolding NF faces with NFV vertices each +template constexpr int unfolded_vertex_count() { + return NFV + (NFV - 2) * (NF - 1); +} + +// For a regular face inscribed on the unit sphere centred at the origin, the +// face's outward unit normal is just the centroid direction, so we don't +// need a separate face_normal helper. +template +inline Vec3 face_centroid(const Vec3 (&vertices)[NV], const int (&face)[NFV]) { + Vec3 c{0.0, 0.0, 0.0}; + for (int k = 0; k < NFV; k++) + c = vec3_add(c, vertices[face[k]]); + return vec3_scale(c, 1.0 / NFV); +} + +// Rigid transform from a polyhedron face's plane to the 2D plane (z=0), +// sending poly_a → net_a and orienting poly_a→poly_b along net_a→net_b. +// poly_a and poly_b must lie on the face plane; only the direction of +// net_b − net_a is used (the 2D distance is derived from the 3D one). +inline Mat4 mat4_rigid_transform(const Vec3 &normal, const Vec3 &poly_a, + const Vec3 &net_a, const Vec3 &poly_b, + const Vec3 &net_b) { + const Vec3 u = vec3_normalize(vec3_subtract(poly_b, poly_a)); + const Vec3 v = vec3_cross(normal, u); + const Vec3 edge_direction = vec3_normalize(vec3_subtract(net_b, net_a)); + const Vec3 r0 = vec3_subtract(vec3_scale(u, edge_direction.x), + vec3_scale(v, edge_direction.y)); + const Vec3 r1 = vec3_add(vec3_scale(u, edge_direction.y), + vec3_scale(v, edge_direction.x)); + const double t0 = net_a.x - vec3_dot(r0, poly_a); + const double t1 = net_a.y - vec3_dot(r1, poly_a); + return {{{r0.x, r0.y, r0.z, t0}, + {r1.x, r1.y, r1.z, t1}, + {0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 1.0}}}; +} + +template +inline double polygon_area_2d(const Vec3 *verts, const int (&face)[NFV]) { + double a = 0.0; + for (int k = 0; k < NFV; k++) + a += vec3_cross(verts[face[k]], verts[face[(k + 1) % NFV]]).z; + return std::fabs(a) * 0.5; +} + +// Unit vector in Z +inline constexpr Vec3 z{0.0, 0.0, 1.0}; + +// Unfold a polyhedron Mesh into a planar net Mesh. +// +// Faces are referenced by 1-based IDs F1..F_NF (matching Snyder's paper +// numbering). `parents[i]` holds the parent of face F(i+1) in the unfolding +// tree as a 1-based face ID, or 0 if F(i+1) is the root. Exactly one entry +// must be 0. The C array indexing stays 0-based; only the values it stores +// are 1-based. +// +// The polyhedron faces must have consistent CCW winding from the outside. +// Each child face is hinged about its shared edge with its parent; the +// algorithm relies on shared edges being traversed in opposite directions +// in the two adjacent faces (which is the case for any closed polyhedron +// with consistent outward winding). +// +// The resulting net is auto-scaled so its total area equals 4π (i.e. the +// area of the unit sphere) and translated so the root face's centroid lies +// at the origin. The root face is laid out "north-up": geographic north +// (+z), projected onto its face plane, points along +y; for a root face +// centred on the rotation axis the slot[0] vertex is placed up instead. +template +inline Mesh(), NF, NFV> +unfold_net(const Mesh &polyhedron, const int (&parents)[NF], + const Vec3 &up_dir = z) { + constexpr int NV_n = unfolded_vertex_count(); + Mesh net{}; + int n_verts = 0; + + // Find and place root face (the entry whose 1-based parent ID is 0). + int root = 0; + while (parents[root] != 0) + root++; + + // Place the root face "up-direction-up": project up_dir onto the root + // face's plane, then orient so that direction lands on +y. up_dir is + // typically geographic north expressed in polyhedron coordinates (the + // third column of the orient matrix), so the root face comes out with + // true geographic north up. Polar fallback (face on the rotation + // axis): aim from the centroid out toward slot[0] instead. + + // Destination in 2d net (origin and unit step in +y) + const Vec3 origin_net = {0.0, 0.0, 0.0}; + const Vec3 offset_net = {0.0, 1.0, 0.0}; + + const Vec3 origin_poly = + face_centroid(polyhedron.vertices, polyhedron.faces[root]); + const Vec3 normal = vec3_normalize(origin_poly); + + // Vector in up_dir direction constrained to the face plane + const Vec3 step_poly = + vec3_subtract(up_dir, vec3_scale(normal, vec3_dot(up_dir, normal))); + Vec3 offset_poly; + if (vec3_length(step_poly) < 1e-12) + // Fallback to 1st vertex at pole + offset_poly = polyhedron.vertices[polyhedron.faces[root][0]]; + else + offset_poly = vec3_add(origin_poly, step_poly); + + const Mat4 root_face_to_plane = mat4_rigid_transform( + normal, origin_poly, origin_net, offset_poly, offset_net); + for (int k = 0; k < NFV; k++) { + Vec3 face_vertex = polyhedron.vertices[polyhedron.faces[root][k]]; + net.vertices[n_verts] = + vec3_apply_mat4(root_face_to_plane, face_vertex); + net.faces[root][k] = n_verts++; + } + + // Initialize queue with direct children of parent face. `root` is a + // 0-based array index; `parents[i]` is a 1-based face ID. + int queue[NF]; + int q_head = 0, q_tail = 0; + for (int i = 0; i < NF; i++) { + if (parents[i] == root + 1) + queue[q_tail++] = i; + } + + // Hinge each child onto its parent at the shared edge. + while (q_head < q_tail) { + const int c = queue[q_head++]; + const int p = parents[c] - 1; + + // Find shared edge + int ia_c = -1, ib_c = -1, ia_p = -1, ib_p = -1; + for (int kc = 0; kc < NFV; kc++) { + for (int kp = 0; kp < NFV; kp++) { + if (polyhedron.faces[c][kc] == polyhedron.faces[p][kp]) { + if (ia_c == -1) { + ia_c = kc; + ia_p = kp; + } else if (ib_c == -1) { + ib_c = kc; + ib_p = kp; + } + } + } + } + assert(ia_c >= 0 && ib_c >= 0 && ia_p >= 0 && ib_p >= 0); + + // Place face onto net using shared edge for orientation + const int na = net.faces[p][ia_p]; + const int nb = net.faces[p][ib_p]; + const Vec3 normal_c = vec3_normalize( + face_centroid(polyhedron.vertices, polyhedron.faces[c])); + const Mat4 face_to_plane = mat4_rigid_transform( + normal_c, polyhedron.vertices[polyhedron.faces[c][ia_c]], + net.vertices[na], polyhedron.vertices[polyhedron.faces[c][ib_c]], + net.vertices[nb]); + for (int k = 0; k < NFV; k++) { + if (k == ia_c) { + net.faces[c][k] = na; + } else if (k == ib_c) { + net.faces[c][k] = nb; + } else { + Vec3 face_vertex = polyhedron.vertices[polyhedron.faces[c][k]]; + net.vertices[n_verts] = + vec3_apply_mat4(face_to_plane, face_vertex); + net.faces[c][k] = n_verts++; + } + } + + // Append children of this face to queue (parents[i] is 1-based). + for (int i = 0; i < NF; i++) { + if (parents[i] == c + 1) + queue[q_tail++] = i; + } + } + + // Rescale so the net's total area equals the unit sphere's (4π). + double total_area = 0.0; + for (int f = 0; f < NF; f++) + total_area += polygon_area_2d(net.vertices, net.faces[f]); + const double scale = std::sqrt(4.0 * M_PI / total_area); + + // Translate so the root face's centroid is at origin, then scale. + const Vec3 centroid = face_centroid(net.vertices, net.faces[root]); + for (int i = 0; i < NV_n; i++) + net.vertices[i] = + vec3_scale(vec3_subtract(net.vertices[i], centroid), scale); + + return net; +} + +} // namespace polyhedral + +#endif // POLYHEDRAL_UNFOLD_H diff --git a/src/sphere.h b/src/sphere.h new file mode 100644 index 0000000000..37c934c16e --- /dev/null +++ b/src/sphere.h @@ -0,0 +1,100 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Spherical geometry helpers on 3D unit vectors — + * slerp, spherical triangle area, scalar triple/quadruple + * products for great-circle math. + * Author: Felix Palmer + * + ****************************************************************************** + * Derived from A5 (Apache-2.0). + * https://github.com/felixpalmer/a5 + * + * modules/utils/vector.ts + * + * Copyright (c) A5 contributors + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + ****************************************************************************/ + +#ifndef SPHERE_H +#define SPHERE_H + +#include "vec3.h" + +#include +#include + +// Numerically stable angular distance measure. +// Returns a value proportional to sin(angle/2) between a and b. +inline double vector_difference(const Vec3 &a, const Vec3 &b) { + Vec3 midpoint_ab = vec3_normalize(vec3_lerp(a, b, 0.5)); + Vec3 cross_result = vec3_cross(a, midpoint_ab); + double d = vec3_length(cross_result); + + if (d < 1e-8) { + // When A and B are very close, sin(x/2) ≈ x/2 + return 0.5 * vec3_length(vec3_subtract(a, b)); + } + return d; +} + +inline double triple_product(const Vec3 &a, const Vec3 &b, const Vec3 &c) { + return vec3_dot(a, vec3_cross(b, c)); +} + +// Great-circle intersection helper. +// Returns dot(b,cross(c,d))*a - dot(a,cross(c,d))*b +inline Vec3 quadruple_product(const Vec3 &a, const Vec3 &b, const Vec3 &c, + const Vec3 &d) { + Vec3 cross_cd = vec3_cross(c, d); + double tp_acd = vec3_dot(a, cross_cd); + double tp_bcd = vec3_dot(b, cross_cd); + return vec3_subtract(vec3_scale(b, tp_acd), vec3_scale(a, tp_bcd)); +} + +inline Vec3 slerp(const Vec3 &a, const Vec3 &b, double t) { + double gamma = vec3_angle(a, b); + if (gamma < 1e-12) { + return vec3_lerp(a, b, t); + } + double sin_gamma = std::sin(gamma); + double wa = std::sin((1.0 - t) * gamma) / sin_gamma; + double wb = std::sin(t * gamma) / sin_gamma; + return vec3_add(vec3_scale(a, wa), vec3_scale(b, wb)); +} + +// Spherical triangle area via midpoint triple product method. +// More numerically stable than vertex-based spherical excess. +inline double spherical_triangle_area(const Vec3 &v1, const Vec3 &v2, + const Vec3 &v3) { + Vec3 mid_a = vec3_normalize(vec3_lerp(v2, v3, 0.5)); + Vec3 mid_b = vec3_normalize(vec3_lerp(v3, v1, 0.5)); + Vec3 mid_c = vec3_normalize(vec3_lerp(v1, v2, 0.5)); + double s = triple_product(mid_a, mid_b, mid_c); + s = std::clamp(s, -1.0, 1.0); + if (std::fabs(s) < 1e-8) { + return 2.0 * s; + } + return 2.0 * std::asin(s); +} + +// Computes acos(1 - 2*x*x) without precision loss for small x. +inline double safe_acos(double x) { + if (x < 1e-3) { + return 2.0 * x + x * x * x / 3.0; + } + return std::acos(1.0 - 2.0 * x * x); +} + +#endif // SPHERE_H diff --git a/src/vec3.h b/src/vec3.h new file mode 100644 index 0000000000..a58d3de5b4 --- /dev/null +++ b/src/vec3.h @@ -0,0 +1,98 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: 3D vector operations for polyhedral projections. + * Author: Felix Palmer + * + ****************************************************************************** + * Derived from gl-matrix (MIT License). + * + * gl-matrix: Copyright (c) 2015-2021, Brandon Jones, Colin MacKenzie IV. + * https://github.com/toji/gl-matrix + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#ifndef VEC3_H +#define VEC3_H + +#include +#include + +struct Vec3 { + double x, y, z; +}; + +inline double vec3_dot(const Vec3 &a, const Vec3 &b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +inline Vec3 vec3_cross(const Vec3 &a, const Vec3 &b) { + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x}; +} + +inline double vec3_length(const Vec3 &v) { + return std::sqrt(v.x * v.x + v.y * v.y + v.z * v.z); +} + +inline Vec3 vec3_normalize(const Vec3 &v) { + double len = vec3_length(v); + if (len == 0.0) + return v; + return {v.x / len, v.y / len, v.z / len}; +} + +inline Vec3 vec3_subtract(const Vec3 &a, const Vec3 &b) { + return {a.x - b.x, a.y - b.y, a.z - b.z}; +} + +inline Vec3 vec3_add(const Vec3 &a, const Vec3 &b) { + return {a.x + b.x, a.y + b.y, a.z + b.z}; +} + +inline Vec3 vec3_scale(const Vec3 &v, double s) { + return {v.x * s, v.y * s, v.z * s}; +} + +inline Vec3 vec3_lerp(const Vec3 &a, const Vec3 &b, double t) { + return {a.x + t * (b.x - a.x), a.y + t * (b.y - a.y), + a.z + t * (b.z - a.z)}; +} + +inline double vec3_angle(const Vec3 &a, const Vec3 &b) { + double d = vec3_dot(a, b) / (vec3_length(a) * vec3_length(b)); + return std::acos(std::clamp(d, -1.0, 1.0)); +} + +struct Mat4 { + double m[4][4]; +}; + +// Apply an affine 4x4 to a 3D point (treated as (v.x, v.y, v.z, 1)). Assumes +// the bottom row is [0, 0, 0, 1] (true affine, no perspective divide). +inline Vec3 vec3_apply_mat4(const Mat4 &m, const Vec3 &v) { + return { + m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z + m.m[0][3], + m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z + m.m[1][3], + m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z + m.m[2][3], + }; +} + +#endif // VEC3_H diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 886bd0ebb8..5fd68aabed 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -3125,106 +3125,176 @@ expect 0.000898315284 0.000000000000 =============================================================================== -# Icosahedral Snyder Equal Area -# Sph +# Tetrahedral Snyder Equal Area +# Sph&Ell =============================================================================== ------------------------------------------------------------------------------- -operation +proj=isea +a=6400000 +operation +proj=tsea +R=6371007.18091875 ------------------------------------------------------------------------------- tolerance 0.2 mm -accept 2 1 -expect -1097074.948153475765139 3442909.309747453313321 -roundtrip 1 -accept 2 -1 -expect -1097074.948149705072865 3233611.728292400948703 +accept 0 0 +expect -7002162.71368800 4042700.52765732 roundtrip 1 -accept -2 1 -expect -1575486.353775786235929 3442168.342736063525081 +accept 45 0 +expect -3607318.74300090 8724734.47739539 roundtrip 1 -accept -2 -1 -expect -1575486.353772019501776 3234352.695310209877789 +accept 2 49 +expect -12216752.96471957 7182742.29453684 roundtrip 1 -operation +proj=isea +mode=hex +resolution=31 -accept 0 0 -expect failure +accept -110 54 +expect -16651351.97836581 4682660.28886354 +roundtrip 1 ------------------------------------------------------------------------------- -operation +proj=isea +R=6371007.18091875 +# +orient_lat / +orient_lon / +azi rotate the polyhedron on the sphere. +# Here the first vertex (apex) lands at the equator on the prime meridian. +operation +proj=tsea +R=6371007.18091875 +orient_lat=0 +orient_lon=0 +azi=0 ------------------------------------------------------------------------------- tolerance 0.2 mm -accept -168.75 58.282525588539 -expect -19186144.870842020958662 3323137.771944524254650 +accept 0 0 +expect -17160609.66302776 9907682.60840717 roundtrip 1 -accept 11.25 58.282525588539 -expect -15348915.896747918799520 9969413.315350906923413 +accept 45 0 +expect 10377221.85600363 8134897.79380681 roundtrip 1 -accept -110 54 -expect -15321401.505530973896384 3338358.859094056300819 +accept 2 49 +expect 13196457.75916024 3260942.81128537 roundtrip 1 -accept -75 45 -expect -12774358.709073608741164 4373188.646695702336729 -roundtrip 1 +------------------------------------------------------------------------------- +# +lat_0 / +lon_0 pan the projected output: (lat_0, lon_0) lands at (0, 0). +operation +proj=tsea +R=6371007.18091875 +lat_0=49 +lon_0=2 +------------------------------------------------------------------------------- +tolerance 0.2 mm accept 2 49 -expect -642252.939347098814324 8796229.009143760427833 +expect 0 0 roundtrip 1 accept 0 0 -expect -1331454.074623266700655 3323137.771634854841977 +expect 5214590.25103157 -3140041.76687951 roundtrip 1 -accept 90 0 -expect 8564460.639100870117545 593869.297485541785136 +------------------------------------------------------------------------------- +# +x_0 / +y_0 stack on top of the +lat_0 / +lon_0 pan. +operation +proj=tsea +R=6371007.18091875 +lat_0=49 +lon_0=2 +x_0=500000 +y_0=1000000 +------------------------------------------------------------------------------- +tolerance 0.2 mm + +accept 2 49 +expect 500000 1000000 roundtrip 1 -accept 0 45 -expect -837334.699958428042009 8323409.759132191538811 +accept 0 0 +expect 5714590.25103157 -2140041.76687951 roundtrip 1 +=============================================================================== +# Dodecahedral Snyder Equal Area +# Sph&Ell +=============================================================================== + ------------------------------------------------------------------------------- -operation +proj=isea +R=6371007.18091875 +orient=pole +operation +proj=dsea +R=6371007.18091875 ------------------------------------------------------------------------------- tolerance 0.2 mm -accept -168.75 58.282525588539 -expect -16702163.549901897087693 6386395.630649688653648 +accept 0 0 +expect 0 -3031675.05748085 roundtrip 1 -accept 11.25 58.282525588539 -expect 619648.646531744743697 6212947.536539182066917 +accept 45 0 +expect 5018607.99435854 -2644265.27628609 roundtrip 1 -accept -110 54 -expect -13285649.857057726010680 6149501.348902118392289 +accept 2 49 +expect 151412.11757341 2437388.28771295 roundtrip 1 -accept -75 45 -expect -7921366.529368571005762 4728387.055336073972285 +accept -110 54 +expect -3813952.97764993 8220035.42468884 roundtrip 1 +------------------------------------------------------------------------------- +# +lat_0 / +lon_0 pan the projected output: (lat_0, lon_0) lands at (0, 0). +operation +proj=dsea +R=6371007.18091875 +lat_0=49 +lon_0=2 +------------------------------------------------------------------------------- +tolerance 0.2 mm + accept 2 49 -expect 152616.434999307675753 5152048.791301283054054 +expect 0 0 roundtrip 1 accept 0 0 -expect 0 -195097.133640714135254 +expect -151412.11757342 -5469063.34519379 roundtrip 1 -accept 90 0 -expect 9593072.435467451811 0 +------------------------------------------------------------------------------- +# +net=a5 selects an alternate unfolding tree (and applies an A5 longitude +# offset of -93°). +operation +proj=dsea +R=6371007.18091875 +net=a5 +------------------------------------------------------------------------------- +tolerance 0.2 mm + +accept 0 0 +expect -6207281.59771581 -1472479.22190498 +roundtrip 1 + +accept 45 0 +expect -663460.15390426 -2932991.72938417 +roundtrip 1 + +=============================================================================== +# Icosahedral Snyder Equal Area +# Sph&Ell +=============================================================================== + +------------------------------------------------------------------------------- +# At default settings on a sphere, the projected origin is the geographic point +# at (lon=11.25°E, lat≈-30.0945°S) — the bbox centre of icosahedron face F8 +# (the central upper-band face), matching Snyder's published Figure 12. +# On an ellipsoid, +orient_lon defaults to 11.20° instead of 11.25°. +------------------------------------------------------------------------------- +operation +proj=isea +R=6371007.18091875 +------------------------------------------------------------------------------- +tolerance 0.2 mm + +accept 0 0 +expect -1331454.07462099 3323137.77178363 +roundtrip 1 + +accept 45 0 +expect 4039031.71269420 3206626.90638330 roundtrip 1 -accept 0 45 -expect 0 4726854.770339427515864 +accept 2 49 +expect -642252.939348975 8796229.00931178 +roundtrip 1 + +accept -110 54 +expect -15321401.5055844 3338358.85897048 +roundtrip 1 + +------------------------------------------------------------------------------- +# +lat_0 / +lon_0 pan the projected output: (lat_0, lon_0) lands at (0, 0). +operation +proj=isea +R=6371007.18091875 +lat_0=49 +lon_0=2 +------------------------------------------------------------------------------- +tolerance 0.2 mm + +accept 2 49 +expect 0 0 +roundtrip 1 + +accept 0 0 +expect -689201.135272012 -5473091.23752815 roundtrip 1 =============================================================================== diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 3fc36c9827..4b44143c74 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -256,6 +256,18 @@ add_test(NAME test_misc COMMAND test_misc) set_property(TEST test_misc PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) +add_executable(test_polyhedral + main.cpp + test_polyhedral.cpp) +target_link_libraries(test_polyhedral + PRIVATE GTest::gtest + PRIVATE ${PROJ_LIBRARIES}) +target_include_directories(test_polyhedral + PRIVATE ${PROJECT_SOURCE_DIR}/src/projections) +add_test(NAME test_polyhedral COMMAND test_polyhedral) +set_property(TEST test_polyhedral + PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) + if (Threads_FOUND AND CMAKE_USE_PTHREADS_INIT AND NOT APPLE) add_definitions(-DPROJ_HAS_PTHREADS) add_executable(test_fork diff --git a/test/unit/test_polyhedral.cpp b/test/unit/test_polyhedral.cpp new file mode 100644 index 0000000000..bdbef2afaa --- /dev/null +++ b/test/unit/test_polyhedral.cpp @@ -0,0 +1,135 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Tests for the polyhedral equal-area projections. + * Author: Felix Palmer + * + ****************************************************************************/ + +#include "gtest_include.h" + +// PROJ include order is sensitive +// clang-format off +#include "proj.h" +#include "proj_internal.h" +// clang-format on + +#include +#include + +namespace { + +struct RoundtripPoint { + double lat; + double lon; +}; + +static const RoundtripPoint test_points[] = { + {0.0, 0.0}, {45.0, 90.0}, {-45.0, -90.0}, {89.0, 0.0}, + {-89.0, 170.0}, {30.0, 120.0}, {-60.0, -150.0}, {15.0, 45.0}, + {-30.0, -60.0}, {70.0, -170.0}, +}; + +struct ProjConfig { + const char *str; + double tol; +}; + +// 1e-13 rad is ~0.6 nm at the equator on Earth — the precision the original +// (no-offset) implementation hits. Configs that compose extra rotations and +// projected-space offsets give up a few bits to cancellation near triangle +// edges; 1e-11 rad is still ~60 μm, well below any practical tolerance. +static const ProjConfig proj_strings[] = { + {"+proj=tsea +R=1", 1e-13}, + {"+proj=tsea +ellps=WGS84", 1e-13}, + {"+proj=tsea +R=1 +orient_lat=90 +orient_lon=0 +azi=0", 1e-13}, + {"+proj=tsea +R=1 +lat_0=10 +lon_0=20", 1e-11}, + {"+proj=tsea +R=1 +orient_lat=45 +orient_lon=30 +azi=60 +lat_0=15 " + "+lon_0=-25 +x_0=1000 +y_0=-2000", + 1e-11}, + {"+proj=dsea +R=1", 1e-13}, + {"+proj=dsea +ellps=WGS84", 1e-13}, + {"+proj=dsea +R=1 +lat_0=10 +lon_0=20", 1e-11}, + {"+proj=isea +R=1", 1e-13}, + {"+proj=isea +ellps=WGS84", 1e-13}, + {"+proj=isea +R=1 +lat_0=10 +lon_0=20", 1e-11}, +}; + +static void roundtrip_test(const char *proj_string, double tolerance) { + PJ *testProj = proj_create(PJ_DEFAULT_CTX, proj_string); + ASSERT_TRUE(testProj != nullptr) << "Failed to create: " << proj_string; + + for (const auto &point : test_points) { + PJ_COORD coord = proj_coord(0, 0, 0, 0); + coord.lp.lam = proj_torad(point.lon); + coord.lp.phi = proj_torad(point.lat); + + PJ_COORD flat = proj_trans(testProj, PJ_FWD, coord); + ASSERT_FALSE(flat.lp.lam == HUGE_VAL) + << "Forward failed for lon=" << point.lon << " lat=" << point.lat; + + PJ_COORD new_coord = proj_trans(testProj, PJ_INV, flat); + ASSERT_FALSE(new_coord.lp.lam == HUGE_VAL) + << "Inverse failed for lon=" << point.lon << " lat=" << point.lat; + + EXPECT_NEAR(coord.lp.lam, new_coord.lp.lam, tolerance) + << "lon roundtrip for lon=" << point.lon << " lat=" << point.lat; + EXPECT_NEAR(coord.lp.phi, new_coord.lp.phi, tolerance) + << "lat roundtrip for lon=" << point.lon << " lat=" << point.lat; + } + + proj_destroy(testProj); +} + +TEST(polyhedral, roundtrip) { + for (const auto &cfg : proj_strings) { + roundtrip_test(cfg.str, cfg.tol); + } +} + +// +lat_0 / +lon_0 must place the supplied geographic point at the projected +// origin (0, 0). +x_0 / +y_0 then stack on top via PROJ core. +TEST(polyhedral, lat0_lon0_at_origin) { + PJ *p = proj_create(PJ_DEFAULT_CTX, "+proj=tsea +R=1 +lat_0=10 +lon_0=20"); + ASSERT_TRUE(p != nullptr); + + PJ_COORD c = proj_coord(0, 0, 0, 0); + c.lp.lam = proj_torad(20.0); + c.lp.phi = proj_torad(10.0); + PJ_COORD f = proj_trans(p, PJ_FWD, c); + EXPECT_NEAR(f.xy.x, 0.0, 1e-13); + EXPECT_NEAR(f.xy.y, 0.0, 1e-13); + + proj_destroy(p); +} + +TEST(polyhedral, dsea_invalid_net_rejected) { + PJ *p = proj_create(PJ_DEFAULT_CTX, "+proj=dsea +R=1 +net=does_not_exist"); + EXPECT_TRUE(p == nullptr); + proj_destroy(p); +} + +TEST(polyhedral, dsea_known_nets_accepted) { + for (const char *net : {"dsea", "a5", "crescent", "flower"}) { + std::string s = std::string("+proj=dsea +R=1 +net=") + net; + PJ *p = proj_create(PJ_DEFAULT_CTX, s.c_str()); + ASSERT_TRUE(p != nullptr) << "Failed to create: " << s; + proj_destroy(p); + } +} + +TEST(polyhedral, x0_y0_stacks_on_lat0_lon0) { + PJ *p = proj_create(PJ_DEFAULT_CTX, + "+proj=tsea +R=1 +lat_0=10 +lon_0=20 +x_0=7 +y_0=-3"); + ASSERT_TRUE(p != nullptr); + + PJ_COORD c = proj_coord(0, 0, 0, 0); + c.lp.lam = proj_torad(20.0); + c.lp.phi = proj_torad(10.0); + PJ_COORD f = proj_trans(p, PJ_FWD, c); + EXPECT_NEAR(f.xy.x, 7.0, 1e-13); + EXPECT_NEAR(f.xy.y, -3.0, 1e-13); + + proj_destroy(p); +} + +} // namespace