diff --git a/docs/source/operations/projections/isea.rst b/docs/source/operations/projections/isea.rst index 73f0af7aca..e01655b095 100644 --- a/docs/source/operations/projections/isea.rst +++ b/docs/source/operations/projections/isea.rst @@ -4,18 +4,24 @@ Icosahedral Snyder Equal Area ******************************************************************************** -Snyder's Icosahedral Equal Area map projections on polyhedral globes for the -dodecahedron and truncated icosahedron offer relatively low scale and -angular distortion. The equations involved are relatively straight-forward, -and for certain instructional tools and databases, the projections are useful -for world maps. The interruptions remain a disadvantage, as with any low-error -projection system applied to the entire globe :cite:`Snyder1992`. - +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. + +.. 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. +---------------------+----------------------------------------------------------+ | **Classification** | Polyhedral, equal area | +---------------------+----------------------------------------------------------+ -| **Available forms** | Forward, spherical | +| **Available forms** | Forward and inverse, spherical | +---------------------+----------------------------------------------------------+ | **Defined area** | Global | +---------------------+----------------------------------------------------------+ @@ -44,27 +50,35 @@ Parameters .. option:: +orient= Can be set to either ``isea`` or ``pole``. See Snyder's Figure 12 for pole orientation :cite:`Snyder1992`. - + *Defaults to isea* .. option:: +azi= Azimuth. + Not supported by the inverse. + *Defaults to 0.0* .. option:: +aperture= + Not supported by the inverse. + *Defaults to 3.0* .. option:: +resolution= + Not supported by the inverse. + *Defaults to 4.0* .. option:: +mode= Can be either ``plane``, ``di``, ``dd`` or ``hex``. - + + Only ``plane`` supported by the inverse. + *Defaults to plane* .. include:: ../options/lon_0.rst diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 8ff5009edc..789bf528af 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -1,6 +1,57 @@ /* - * This code was entirely written by Nathan Wagner - * and is in the public domain. + 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 @@ -10,6 +61,7 @@ #include #include +#include #include #include "proj.h" @@ -30,17 +82,52 @@ /* 26.565051177 degrees */ #define V_LAT 0.46364760899944494524 -/* 52.62263186 */ -#define E_RAD 0.91843818702186776133 +// 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) -/* 10.81231696 */ -#define F_RAD 0.18871053072122403508 +// tang * sin(60 deg) +#define TABLE_G (tang * sin60) -/* R tan(g) sin(60) */ -#define TABLE_G 0.6615845383 +// (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.1909830056 +#define TABLE_H (0.25 * tang) /* in radians */ #define ISEA_STD_LAT 1.01722196792335072101 @@ -135,9 +222,9 @@ static void hexbin2(double width, double x, double y, long *i, long *j) { *j = h.y; } +#define numIcosahedronFaces 20 + namespace { // anonymous namespace -enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 }; -enum isea_topology { ISEA_HEXAGON = 6, ISEA_TRIANGLE = 3, ISEA_DIAMOND = 4 }; enum isea_address_form { ISEA_GEO, ISEA_Q2DI, @@ -149,9 +236,7 @@ enum isea_address_form { ISEA_VERTEX2DD, ISEA_HEX }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_dgg { int polyhedron; /* ignored, icosahedron */ double o_lat, o_lon, o_az; /* orientation, radians */ @@ -164,23 +249,22 @@ struct isea_dgg { int quad; /* quad of last transformed point */ unsigned long serial; }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_pt { double x, y; }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_geo { double longitude, lat; }; + } // anonymous namespace -/* ENDINC */ +/* No longer used + +enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = numIcosahedronFaces }; +enum isea_topology { ISEA_HEXAGON = 6, ISEA_TRIANGLE = 3, ISEA_DIAMOND = 4 }; -namespace { // anonymous namespace enum snyder_polyhedron { SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON, @@ -190,17 +274,14 @@ enum snyder_polyhedron { SNYDER_POLY_DODECAHEDRON, SNYDER_POLY_ICOSAHEDRON }; -} // anonymous namespace -namespace { // anonymous namespace struct snyder_constants { double g, G, theta; - /* cppcheck-suppress unusedStructMember */ + // cppcheck-suppress unusedStructMember double ea_w, ea_a, ea_b, g_w, g_a, g_b; }; -} // anonymous namespace -/* TODO put these in radians to avoid a later conversion */ +// TODO put these in radians to avoid a later conversion static const struct snyder_constants constants[] = { {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0}, {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027}, @@ -210,6 +291,7 @@ static const struct snyder_constants constants[] = { {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0}, }; +*/ static struct isea_geo vertex[] = { {0.0, DEG90}, {DEG180, V_LAT}, {-DEG108, V_LAT}, {-DEG36, V_LAT}, @@ -250,9 +332,8 @@ static double az_adjustment(int triangle) { static struct isea_pt isea_triangle_xy(int triangle) { struct isea_pt c; - const double Rprime = 0.91038328153090290025; - triangle = (triangle - 1) % 20; + triangle = (triangle - 1) % numIcosahedronFaces; c.x = TABLE_G * ((triangle % 5) - 2) * 2.0; if (triangle > 9) { @@ -275,8 +356,8 @@ static struct isea_pt isea_triangle_xy(int triangle) { /* should be impossible */ exit(EXIT_FAILURE); } - c.x *= Rprime; - c.y *= Rprime; + c.x *= RprimeOverR; + c.y *= RprimeOverR; return c; } @@ -302,47 +383,20 @@ static double sph_azimuth(double f_lon, double f_lat, double t_lon, static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { int i; - /* - * spherical distance from center of polygon face to any of its - * vertices on the globe - */ - double g; - - /* - * spherical angle between radius vector to center and adjacent edge - * of spherical polygon on the globe - */ - double G; - - /* - * plane angle between radius vector to center and adjacent edge of - * plane polygon - */ - double theta; - /* additional variables from snyder */ double q, H, Ag, Azprime, Az, dprime, f, rho, x, y; /* variables used to store intermediate results */ - double cot_theta, tan_g, az_offset; + double az_offset; /* how many multiples of 60 degrees we adjust the azimuth */ int Az_adjust_multiples; - struct snyder_constants c; - /* * TODO by locality of reference, start by trying the same triangle * as last time */ - - /* TODO put these constants in as radians to begin with */ - c = constants[SNYDER_POLY_ICOSAHEDRON]; - theta = PJ_TORAD(c.theta); - g = PJ_TORAD(c.g); - G = PJ_TORAD(c.G); - - for (i = 1; i <= 20; i++) { + for (i = 1; i <= numIcosahedronFaces; i++) { double z; struct isea_geo center; @@ -353,7 +407,7 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { cos(center.lat) * cos(ll->lat) * cos(ll->longitude - center.longitude)); /* not on this triangle */ - if (z > g + 0.000005) { /* TODO DBL_EPSILON */ + if (z > sdc2vos /*g*/ + 0.000005) { /* TODO DBL_EPSILON */ continue; } @@ -391,12 +445,9 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { } /* step 3 */ - cot_theta = 1.0 / tan(theta); - tan_g = tan(g); /* TODO this is a constant */ /* Calculate q from eq 9. */ - /* TODO cot_theta is cot(30) */ - q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta); + q = atan2(tang, cos(Az) + sin(Az) * cotTheta); /* not in this triangle */ if (z > q + 0.000005) { @@ -407,30 +458,29 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { /* Apply equations 5-8 and 10-12 in order */ /* eq 5 */ - /* Rprime = 0.9449322893 * R; */ - /* R' in the paper is for the truncated */ - const double Rprime = 0.91038328153090290025; + /* R' in the paper is for the truncated (icosahedron?) */ /* eq 6 */ - H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G)); + H = acos(sin(Az) * sinGcosSDC2VoS /* sin(G) * cos(g) */ - + cos(Az) * cosG); /* eq 7 */ /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */ - Ag = Az + G + H - DEG180; + Ag = Az + DEG_TO_RAD * 36 /* G */ + H - DEG180; /* eq 8 */ - Azprime = atan2(2.0 * Ag, - Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta); + Azprime = atan2(2.0 * Ag, RprimeOverR * RprimeOverR * tang * tang - + 2.0 * Ag * cotTheta); /* eq 10 */ /* cot(theta) = 1.73205080756887729355 */ - dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta); + dprime = RprimeOverR * tang / (cos(Azprime) + sin(Azprime) * cotTheta); /* eq 11 */ - f = dprime / (2.0 * Rprime * sin(q / 2.0)); + f = dprime / (2.0 * RprimeOverR * sin(q / 2.0)); /* eq 12 */ - rho = 2.0 * Rprime * f * sin(z / 2.0); + rho = 2.0 * RprimeOverR * f * sin(z / 2.0); /* * add back the same 60 degree multiple adjustment from step @@ -561,7 +611,7 @@ static int isea_grid_init(struct isea_dgg *g) { if (!g) return 0; - g->polyhedron = 20; + g->polyhedron = numIcosahedronFaces; g->o_lat = ISEA_STD_LAT; g->o_lon = ISEA_STD_LONG; g->o_az = 0.0; @@ -984,8 +1034,33 @@ static struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in) { PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph"; namespace { // anonymous namespace + +struct GeoPoint { + double lat, lon; +}; // In radians + +class ISEAPlanarProjection; + +struct isea_inverse_params { + 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); +}; + struct pj_isea_data { struct isea_dgg dgg; + struct isea_inverse_params inv; }; } // anonymous namespace @@ -999,6 +1074,9 @@ static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ in.lat = lp.phi; try { + // TODO: Convert geodetic latitude to authalic latitude if not + // spherical + // as in eqearth, healpix, laea, etc. out = isea_forward(&Q->dgg, &in); } catch (const char *) { proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); @@ -1011,6 +1089,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,9 +1102,10 @@ 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; + /* P->dgg.radius = P->a; / * otherwise defaults to 1 */ /* calling library will scale, I think */ @@ -1089,9 +1170,377 @@ PJ *PJ_PROJECTION(isea) { Q->dgg.aperture = 3; } + Q->inv.initialize(P); + return P; } +#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 +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}}; + +// 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 isea_inverse_params ¶ms, GeoPoint &result) { + bool r = false; + static const double epsilon = + 2E-8; // NOTE: 1E-11 seems too small for forward projection + // precision at boundaries + 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 isea_inverse_params ¶ms, 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; + { + 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 latSin = + sinLat0 * cosZ + cosLat0SinZ * cos(Az_earth); + // REVIEW: Is there a way we could use atan2() + // for latitude as well for better precision? + double lat = fabs(latSin - 1.0) < 1E-15 ? M_PI / 2 + : fabs(latSin + 1.0) < 1E-15 + ? -M_PI / 2 + : asin(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 isea_inverse_params::initialize(const PJ *P) { + struct pj_isea_data *Q = static_cast(P->opaque); + // Only supporting default planar options for now + if (Q->dgg.output == ISEA_PLANE && Q->dgg.o_az == 0.0 && + Q->dgg.aperture == 3.0 && Q->dgg.resolution == 4.) { + // Only supporting +orient=isea and +orient=pole for now + if (Q->dgg.o_lat == ISEA_STD_LAT && Q->dgg.o_lon == ISEA_STD_LONG) + p = &standardISEA; + else if (Q->dgg.o_lat == M_PI / 2.0 && Q->dgg.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) { + struct pj_isea_data *Q = static_cast(P->opaque); + ISEAPlanarProjection *p = Q->inv.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->inv.xo, xy.y * P->a + Q->inv.yo}; + GeoPoint result; + + if (p->cartesianToGeo(input, Q->inv, 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 @@ -1101,9 +1550,5 @@ PJ *PJ_PROJECTION(isea) { #undef DEG180 #undef ISEA_SCALE #undef V_LAT -#undef E_RAD -#undef F_RAD #undef TABLE_G #undef TABLE_H -#undef ISEA_STD_LAT -#undef ISEA_STD_LONG diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 98ea577449..866f2f4d52 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -2762,20 +2762,64 @@ expect 0.000898315284 0.000000000000 ------------------------------------------------------------------------------- operation +proj=isea +a=6400000 ------------------------------------------------------------------------------- -tolerance 0.1 mm +tolerance 1 mm accept 2 1 -expect -1097074.948022474 3442909.309037183 +expect -1097074.948153483448550 3442909.309747451916337 +roundtrip 1 + accept 2 -1 -expect -1097074.948264795 3233611.728585708 +expect -1097074.948149696690962 3233611.728292400948703 +roundtrip 1 + accept -2 1 -expect -1575486.353641554 3442168.342028188 +expect -1575486.353775785770267 3442168.342736063525081 +roundtrip 1 + accept -2 -1 -expect -1575486.353880283 3234352.695594706 +expect -1575486.353772019501776 3234352.695310209877789 +roundtrip 1 operation +proj=isea +mode=hex +resolution=31 accept 0 0 expect failure +------------------------------------------------------------------------------- +operation +proj=isea +R=6371007.18091875 +------------------------------------------------------------------------------- +tolerance 1 mm + +accept -168.75 58.282525588539 +expect -19186144.870842020958662 3323137.771944524254650 +roundtrip 1 + +accept 11.25 58.282525588539 +expect -15348915.896747924387455 9969413.315350895747542 +roundtrip 1 + +accept -110 54 +expect -15321401.505530973896384 3338358.859094054903835 +roundtrip 1 + +accept -75 45 +expect -12774358.709073614329100 4373188.646695707924664 +roundtrip 1 + +accept 2 49 +expect -642252.939347098581493 8796229.009143756702542 +roundtrip 1 + +accept 0 0 +expect -1331454.074623258318752 3323137.771634854841977 +roundtrip 1 + +accept 90 0 +expect 8564460.639100877568126 593869.297485543764196 +roundtrip 1 + +accept 0 45 +expect -837334.699958425830118 8323409.759132194332778 +roundtrip 1 + =============================================================================== # Kavrayskiy V # PCyl., Sph. @@ -5988,7 +6032,7 @@ accept 1 0.5 expect 45 0 accept 0 0 expect -45 -35.446011426401625 -accept 1 0 +accept 1 0 expect 45 -35.446011426401625 accept 1 1 expect 45 35.446011426401625