From e734459c28ed7914562c41e6a8062dc66c4d8a1b Mon Sep 17 00:00:00 2001 From: Jerome St-Louis Date: Wed, 31 Jul 2024 11:34:41 -0400 Subject: [PATCH] ISEA Inverse projection (#3047) - Integrating code ported from Java to eC, then eC to C++ implementing the inverse ISEA projection - Originally from 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 (MIT License) - Draft implementation (for now only implemented for +proj=isea +R=6371007.18091875 with no other parameters) --- src/projections/isea.cpp | 399 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 399 insertions(+) diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 8ff5009edc..2e3abe002e 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -1011,6 +1011,8 @@ static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ 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( @@ -1022,6 +1024,7 @@ PJ *PJ_PROJECTION(isea) { // 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->dgg); Q->dgg.output = ISEA_PLANE; @@ -1107,3 +1110,399 @@ PJ *PJ_PROJECTION(isea) { #undef TABLE_H #undef ISEA_STD_LAT #undef ISEA_STD_LONG + +/* + Icosahedron Snyder equal-area (ISEA) Inverse Projection + -------------------------------------------------------------------------- + This 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 + */ + +struct GeoPoint { double lat, lon; }; // In radians + +#define Pi 3.1415926535897932384626433832795028841971 + +#define Degrees(x) ((x) * Pi / 180) + +#define wgs84InvFlattening 298.257223563 +#define wgs84Major 6378137.0 // Meters +#define wgs84Minor (wgs84Major - (wgs84Major / wgs84InvFlattening)) // 6356752.3142451792955399 + +// #define phi ((1 + sqrt(5)) / 2) +// radius +// #define a2 (wgs84Major * wgs84Major) +// #define c2 (wgs84Minor * wgs84Minor) + +// #define eccentricity 0.0818191908426 // sqrt(1.0 - (c2/a2)); +// #define log1pe_1me 0.1640050079196 // log((1+e)/(1-e))) +// #define S (Pi * (2 * a2 + c2/eccentricity * log1pe_1me)) +#define earthAuthalicRadius 6371007.18091875 //sqrt(S / (4 * Pi)); // R -- authalic sphere radius for WGS84 [6371007.1809184747 m] +#define R2 (earthAuthalicRadius * earthAuthalicRadius) // R^2 +#define RprimeOverR 0.9103832815095032 // (1 / (2 * sqrt(5)) + 1 / 6.0) * sqrt(Pi * sqrt(3)); // R' / R +#define Rprime (RprimeOverR * earthAuthalicRadius) // R' + +#define Min std::min +#define Max std::max + +#define inf std::numeric_limits::infinity() + +// distortion +// static double maximumAngularDistortion = 17.27; +// static double maximumScaleVariation = 1.163; +// static double minimumScaleVariation = .860; + +// Vertices of dodecahedron centered in icosahedron faces +#define E Degrees(52.6226318593487) //(Degrees)atan((3 + sqrt(5)) / 4); // Latitude of center of top icosahedron faces +#define F Degrees(10.8123169635739) //(Degrees)atan((3 - sqrt(5)) / 4); // Latitude of center of faces mirroring top icosahedron face +#define numIcosahedronFaces 20 +static const GeoPoint facesCenterDodecahedronVertices[numIcosahedronFaces] = +{ + { E, Degrees(-144) }, { E, Degrees(-72) }, { E, Degrees( 0) }, { E, Degrees( 72) }, { E, Degrees(144) }, + { F, Degrees(-144) }, { F, Degrees(-72) }, { F, Degrees( 0) }, { F, Degrees( 72) }, { F, Degrees(144) }, + { -F, Degrees(-108) }, { -F, Degrees(-36) }, { -F, Degrees(36) }, { -F, Degrees(108) }, { -F, Degrees(180) }, + { -E, Degrees(-108) }, { -E, Degrees(-36) }, { -E, Degrees(36) }, { -E, Degrees(108) }, { -E, Degrees(180) } +}; +// static define precision = Degrees(1e-9); +#define precision Degrees(1e-11) +#define precisionPerDefinition Degrees(1e-5) + +#define Rprime2X (2 * Rprime) +#define AzMax Degrees(120) +#define sdc2vos (F + 2 * Degrees(58.2825255885418) /*(Degrees)atan(phi)*/ - Degrees(90)) // g -- sphericalDistanceFromCenterToVerticesOnSphere +#define tang 0.763932022500419 //tan(sdc2vos) +#define cotTheta (1.0 / tan(Degrees(30))) +#define RprimeTang (Rprime * tang) // twice the center-to-base distance +#define centerToBase (RprimeTang / 2) +#define triWidth (RprimeTang * 1.73205080756887729352744634150587236694280525381038) //sqrt(3) +#define Rprime2Tan2g (RprimeTang * RprimeTang) +#define cosG cos(Degrees(36)) +#define sinGcosSDC2VoS sin(Degrees(36)) * cos(sdc2vos) // sin G cos g +#define westVertexLon Degrees(-144) + +static const double yOffsets[4] = { -2*centerToBase, -4*centerToBase, -5*centerToBase, -7*centerToBase }; + +static inline double latGeocentricToGeodetic(double theta) +{ + static const double a2ob2 = (wgs84Major * wgs84Major) / (wgs84Minor * wgs84Minor); + return atan(tan(theta) * a2ob2); +} + +struct ISEAFacePoint +{ + int face; + double x, y; +}; + +class ISEAPlanarProjection +{ +public: + ISEAPlanarProjection(const GeoPoint & value) + { + orientation = value; + sinOrientationLat = sin(value.lat); + cosOrientationLat = cos(value.lat); + } + + bool cartesianToGeo(const PJ_XY & position, GeoPoint & result) + { + bool r = false; + static const double epsilon = 1E-11; //0.000000001; + int face = 0; + // Rotate and shear to determine face if not stored in position.z + static const double sr = -0.86602540378443864676372317075293618347140262690519, cr = 0.5; // sin(-60), cos(-60) + static const double shearX = 1.0 / 1.73205080756887729352744634150587236694280525381038; //sqrt(3); // 0.5773502691896257 -- 2*centerToBase / triWidth + static const double sx = 1.0 / triWidth, sy = 1.0 / (3*centerToBase); + double yp = -(position.x * sr + position.y * cr); + double x = (position.x * cr - position.y * sr + yp * shearX) * sx; + double y = yp * sy; + + 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) + { + // Degrees faceLon = facesCenterDodecahedronVertices[face].lon; + int fy = (face-1) / 5, fx = (face-1) - 5*fy; + double x = position.x - (2*fx + fy/2 + 1) * triWidth/2; + double y = position.y - (yOffsets[fy] + 3 * centerToBase); + GeoPoint dst; + + r = icosahedronToSphere({ face - 1, x, y }, dst); + + // REVIEW: It seems like the forward transformation +proj=isea +R=6371007.18091875 does not apply this geodetic to geocentric conversion + // In FMSDI phase 3, this conversion was necessary to align data from PYXIS server properly. Should this be applied or not? + // dst.lat = latGeocentricToGeodetic(dst.lat); + + if(dst.lon < -Pi - epsilon) dst.lon += 2*Pi; + else if(dst.lon > Pi + epsilon) dst.lon -= 2*Pi; + + result = { dst.lat, dst.lon }; + } + return r; + } + + // Converts coordinates on the icosahedron to geographic coordinates (inverse projection) + bool icosahedronToSphere(const ISEAFacePoint & c, 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 = Rprime2Tan2g / (2 * (cotAz + cotTheta)); // A_G or A_{ABD} + double deltaAz = 10 * precision; + double degAreaOverR2Plus180Minus36 = area / 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 = RprimeTang / (cosAz + sinAz * cotTheta); // d' + double f = d / (Rprime2X * sin(q / 2)); // f + double z = 2 * asin(rho / (Rprime2X * f)); + + Az_earth -= AzAdjustment; + { + double lat0 = facesCenterDodecahedronVertices[c.face].lat; + double sinLat0 = sin(lat0), cosLat0 = cos(lat0); + double sinZ = sin(z), cosZ = cos(z); + double cosLat0SinZ = cosLat0 * sinZ; + double lat = asin(sinLat0 * cosZ + cosLat0SinZ * cos(Az_earth)); + 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 < Degrees(-90) + precisionPerDefinition || c.lat > Degrees(90) - precisionPerDefinition) ? 0 : c.lon; + if(orientation.lat || orientation.lon) + { + 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 void guessFace(const GeoPoint & c, int result[2]) + { + double lat = c.lat, lon = c.lon; + if(lat > E - F) + { + if (lon < Degrees( -108 )) result[0] = 0, result[1] = 5; + else if (lon < Degrees( -36 )) result[0] = 1, result[1] = 6; + else if (lon < Degrees( 36 )) result[0] = 2, result[1] = 7; + else if (lon < Degrees( 108 )) result[0] = 3, result[1] = 8; + else result[0] = 4, result[1] = 9; + } + else if(lat < F - E) + { + if (lon < Degrees( -144 )) result[0] = 19, result[1] = 14; + else if (lon < Degrees( -72 )) result[0] = 15, result[1] = 10; + else if (lon < Degrees( 0 )) result[0] = 16, result[1] = 11; + else if (lon < Degrees( 72 )) result[0] = 17, result[1] = 12; + else if (lon < Degrees( 144 )) result[0] = 18, result[1] = 13; + else result[0] = 19, result[1] = 14; + } + else + { + if (lon < Degrees( -144 )) result[0] = 5, result[1] = 14; + else if (lon < Degrees( -108 )) result[0] = 5, result[1] = 10; + else if (lon < Degrees( -72 )) result[0] = 6, result[1] = 10; + else if (lon < Degrees( -36 )) result[0] = 6, result[1] = 11; + else if (lon < Degrees( 0 )) result[0] = 7, result[1] = 11; + else if (lon < Degrees( 36 )) result[0] = 7, result[1] = 12; + else if (lon < Degrees( 72 )) result[0] = 8, result[1] = 12; + else if (lon < Degrees( 108 )) result[0] = 8, result[1] = 13; + else if (lon < Degrees( 144 )) result[0] = 9, result[1] = 13; + else result[0] = 9, result[1] = 14; + } + } + + static inline double faceOrientation(int face) + { + return (face <= 4 || (10 <= face && face <= 14)) ? 0 : Degrees(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({ (E + F) / 2 /* Degrees(90 - 58.282525590) = 31.7174744114613 */, Degrees(-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 }); + +// NOTE: This currently assumes +proj=isea +R=6371007.18091875 (authalic sphere with surface area of WGS84) +static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P) +{ + PJ_XY input { xy.x * wgs84Major, xy.y * wgs84Major }; + GeoPoint result; + + /* + Alaska Left Vertex + OGC:1534 x = 0 y = 0 + WGS84 58.454376167463 N 168.75 W + +proj=isea x = -19175223.42 y = 3342054.28 + + Top polar vertices triWidth/2 triWidth*sqrt(3)/2 + OGC:1534 x = 3837228.9741880306974053 y = 6646275.5435691606253386 (five points) + x = 11511686.9225640930235386 y = 6646275.5435691606253386 + x = 19186144.8709401525557041 y = 6646275.5435691606253386 + x = 26860602.8193162158131599 y = 6646275.5435691606253386 + x = 34535060.7676922753453255 y = 6646275.5435691606253386 + WGS84 58.454376167463 N 11.25 W + +proj=isea x = -8262744.20 y = 8638708.65 + + Base of left tri triWidth/2 + OGC:1534 x = 3837228.9741880306974053 y = 0 + WGS84 54.1828156920591 N 110.4674744114474 W + +proj=isea x = -15360252.72 y = 3339436.98 + + Ottawa + OGC:1534 x = 6431690.5052938889712 y = 1039748.3954550651833 + WGS84 45 N 75 W + +proj=isea x = -12774358.71 y = 4373188.65 + + Paris + OGC:1534 x = 18541657.1754211001098 y = 5450744.7171042216942 + WGS84 49 N 2 E + +proj=isea x = -642252.94 y = 8796229.01 + + 0 N, 0 E + OGC:1534 x = 17854690.7963187731802 y = 0 + WGS84 0 N 0 E + +proj=isea x = -1331454.07 y = 3323137.77 + */ + // triWidth: 7674457.9483760613948107 + // centerToBase: 2215425.1811896911822259 + + // Default origin of +proj=isea is different + static const double xo = 2.5 * triWidth; + static const double yo = -1.5 * centerToBase; + + input.x = xy.x * earthAuthalicRadius /*wgs84Major*/ + xo; + input.y = xy.y * earthAuthalicRadius /*wgs84Major*/ + yo; + + if(standardISEA.cartesianToGeo(input, result)) + // NOTE: Use polarISEA instead for +orient=pole + return { .lam = result.lon, .phi = result.lat }; + else + return { .lam = inf, .phi = inf }; +}