From 2e104e092578347de11208e9a3a80a3bf711265d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 25 Sep 2020 22:53:38 +0200 Subject: Implement ellipsoidal formulation of +proj=ortho (fixes #397) - Map ESRI 'Local' to +proj=ortho when Scale_Factor = 1 and Azimuth = 0 - Map ESRI 'Orthographic' to a PROJ WKT2 'Orthographic (Spherical)' which maps to +proj=ortho +f=0 to froce spherical evaluation --- src/projections/ortho.cpp | 87 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 81 insertions(+), 6 deletions(-) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index 5f9366c3..ff9d9341 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -4,7 +4,7 @@ #include "proj_internal.h" #include -PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph"; +PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph&Ell"; namespace { // anonymous namespace enum Mode { @@ -19,6 +19,7 @@ namespace { // anonymous namespace struct pj_opaque { double sinph0; double cosph0; + double nu0; enum Mode mode; }; } // anonymous namespace @@ -121,6 +122,72 @@ static PJ_LP ortho_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, invers } +static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ + PJ_XY xy; + struct pj_opaque *Q = static_cast(P->opaque); + + // From EPSG guidance note 7.2 + const double cosphi = cos(lp.phi); + const double sinphi = sin(lp.phi); + const double coslam = cos(lp.lam); + const double sinlam = sin(lp.lam); + const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi); + xy.x = nu * cosphi * sinlam; + xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) + + P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0; + + return xy; +} + + + +static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ + PJ_LP lp; + struct pj_opaque *Q = static_cast(P->opaque); + + // From EPSG guidance note 7.2 + + // It suggests as initial guess: + // lp.lam = 0; + // lp.phi = P->phi0; + // But for poles, this will not converge well. Better use: + lp = ortho_s_inverse(xy, P); + + for( int i = 0; i < 20; i++ ) + { + const double cosphi = cos(lp.phi); + const double sinphi = sin(lp.phi); + const double coslam = cos(lp.lam); + const double sinlam = sin(lp.lam); + const double one_minus_es_sinphi2 = 1.0 - P->es * sinphi * sinphi; + const double nu = 1.0 / sqrt(one_minus_es_sinphi2); + PJ_XY xy_new; + xy_new.x = nu * cosphi * sinlam; + xy_new.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) + + P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0; + const double rho = (1.0 - P->es) * nu / one_minus_es_sinphi2; + const double J11 = -rho * sinphi * sinlam; + const double J12 = nu * cosphi * coslam; + const double J21 = rho * (cosphi * Q->cosph0 + sinphi * Q->sinph0 * coslam); + const double J22 = nu * Q->sinph0 * Q->cosph0 * sinlam; + const double D = J11 * J22 - J12 * J21; + const double dx = xy.x - xy_new.x; + const double dy = xy.y - xy_new.y; + const double dphi = (J22 * dx - J12 * dy) / D; + const double dlam = (-J21 * dx + J11 * dy) / D; + lp.phi += dphi; + if( lp.phi > M_PI_2) lp.phi = M_PI_2; + else if( lp.phi < -M_PI_2) lp.phi = -M_PI_2; + lp.lam += dlam; + if( fabs(dphi) < 1e-12 && fabs(dlam) < 1e-12 ) + { + return lp; + } + } + pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT); + return lp; +} + PJ *PROJECTION(ortho) { struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); @@ -128,17 +195,25 @@ PJ *PROJECTION(ortho) { return pj_default_destructor(P, ENOMEM); P->opaque = Q; + Q->sinph0 = sin(P->phi0); + Q->cosph0 = cos(P->phi0); if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10) Q->mode = P->phi0 < 0. ? S_POLE : N_POLE; else if (fabs(P->phi0) > EPS10) { Q->mode = OBLIQ; - Q->sinph0 = sin(P->phi0); - Q->cosph0 = cos(P->phi0); } else Q->mode = EQUIT; - P->inv = ortho_s_inverse; - P->fwd = ortho_s_forward; - P->es = 0.; + if( P->es == 0 ) + { + P->inv = ortho_s_inverse; + P->fwd = ortho_s_forward; + } + else + { + Q->nu0 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0); + P->inv = ortho_e_inverse; + P->fwd = ortho_e_forward; + } return P; } -- cgit v1.2.3 From 3d88b6fc89f95803bfd4d59c47eef1c05bda710c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 26 Sep 2020 14:40:29 +0200 Subject: Ortho: add visibility condition for ellipsoidal case. Credits to @cffk --- src/projections/ortho.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index ff9d9341..0594205b 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -50,6 +50,13 @@ static PJ_XY ortho_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forwar break; case OBLIQ: sinphi = sin(lp.phi); + + // Is the point visible from the projection plane ? + // From https://lists.osgeo.org/pipermail/proj/2020-September/009831.html + // this is the dot product of the normal of the ellipsoid at the center of + // the projection and at the point considered for projection. + // [cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)] + // Also from Snyder's Map Projection - A working manual, equation (5-3), page 149 if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) return forward_error(P, lp, xy); xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam; @@ -131,6 +138,14 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa const double sinphi = sin(lp.phi); const double coslam = cos(lp.lam); const double sinlam = sin(lp.lam); + + // Is the point visible from the projection plane ? + // Same condition as in spherical case + if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) { + xy.x = HUGE_VAL; xy.y = HUGE_VAL; + return forward_error(P, lp, xy); + } + const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi); xy.x = nu * cosphi * sinlam; xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) + -- cgit v1.2.3 From b04894819ea4e4d9d93e03015f0c7c9aa84642fe Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 26 Sep 2020 18:45:25 +0200 Subject: Ortho ellipsoidal inverse: add non iterative implementations for polar and equatorial --- src/projections/ortho.cpp | 60 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index 0594205b..ac54d88a 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -160,6 +160,66 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver PJ_LP lp; struct pj_opaque *Q = static_cast(P->opaque); + if( Q->mode == N_POLE || Q->mode == S_POLE ) + { + // Polar case. Forward case equations can be simplified as: + // xy.x = nu * cosphi * sinlam + // xy.y = nu * -cosphi * coslam * sign(phi0) + // ==> lam = atan2(xy.x, -xy.y * sign(phi0)) + // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2 + // rh^2 = cosphi^2 / (1 - es * sinphi^2) + // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) + + const double rh2 = xy.x * xy.x + xy.y * xy.y; + if (rh2 >= 1. - 1e-15) { + if ((rh2 - 1.) > EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + lp.phi = 0; + } + else { + lp.phi = asin(sqrt((1 - rh2) / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); + } + lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1)); + return lp; + } + + if( Q->mode == EQUIT ) + { + // Equatorial case. Forward case equations can be simplified as: + // xy.x = nu * cosphi * sinlam + // xy.y = nu * sinphi * (1 - P->es) + // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2 + // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2 + + const auto SQ = [](double x) { return x*x; }; + + // Equation of the ellipse + if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + + const double sinphi2 = xy.y == 0 ? 0 : 1.0 / (SQ((1 - P->es) / xy.y) + P->es); + if( sinphi2 > 1 - 1e-11 ) { + lp.phi = M_PI_2 * (xy.y > 0 ? 1 : -1); + lp.lam = 0; + return lp; + } + lp.phi = asin(sqrt(sinphi2)) * (xy.y > 0 ? 1 : -1); + const double sinlam = xy.x * sqrt((1 - P->es * sinphi2) / (1 - sinphi2)); + if( fabs(sinlam) - 1 > -1e-15 ) + lp.lam = M_PI_2 * (xy.x > 0 ? 1: -1); + else + lp.lam = asin(sinlam); + return lp; + } + // From EPSG guidance note 7.2 // It suggests as initial guess: -- cgit v1.2.3 From fedeeec68ff6a65126da35ae54ec75a719ff40ce Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 27 Sep 2020 13:27:01 +0200 Subject: Ortho ellipsoidal inverse: add domain check for oblique case, and slighly improve initial guessing --- src/projections/ortho.cpp | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index ac54d88a..b4ecff7f 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -20,6 +20,8 @@ struct pj_opaque { double sinph0; double cosph0; double nu0; + double y_shift; + double y_scale; enum Mode mode; }; } // anonymous namespace @@ -155,11 +157,12 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa } - static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ PJ_LP lp; struct pj_opaque *Q = static_cast(P->opaque); + const auto SQ = [](double x) { return x*x; }; + if( Q->mode == N_POLE || Q->mode == S_POLE ) { // Polar case. Forward case equations can be simplified as: @@ -170,7 +173,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // rh^2 = cosphi^2 / (1 - es * sinphi^2) // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) - const double rh2 = xy.x * xy.x + xy.y * xy.y; + const double rh2 = SQ(xy.x) + SQ(xy.y); if (rh2 >= 1. - 1e-15) { if ((rh2 - 1.) > EPS10) { proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); @@ -195,8 +198,6 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2 // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2 - const auto SQ = [](double x) { return x*x; }; - // Equation of the ellipse if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) { proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); @@ -220,13 +221,26 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver return lp; } + // Using Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam == 0 (visibity + // condition of the forward case) in the forward equations, and a lot of + // substition games... + PJ_XY xy_recentered; + xy_recentered.x = xy.x; + xy_recentered.y = (xy.y - Q->y_shift) / Q->y_scale; + if( SQ(xy.x) + SQ(xy_recentered.y) > 1 + 1e-11 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + // From EPSG guidance note 7.2 // It suggests as initial guess: // lp.lam = 0; // lp.phi = P->phi0; // But for poles, this will not converge well. Better use: - lp = ortho_s_inverse(xy, P); + lp = ortho_s_inverse(xy_recentered, P); for( int i = 0; i < 20; i++ ) { @@ -286,6 +300,8 @@ PJ *PROJECTION(ortho) { else { Q->nu0 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0); + Q->y_shift = P->es * Q->nu0 * Q->sinph0 * Q->cosph0; + Q->y_scale = 1.0 / sqrt(1.0 - P->es * Q->cosph0 * Q->cosph0); P->inv = ortho_e_inverse; P->fwd = ortho_e_forward; } -- cgit v1.2.3 From 5e7d58ec69764f7ffa8d09991aa0b0102d2f02ab Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 27 Sep 2020 15:58:54 +0200 Subject: Ortho ellipsoidal inverse: improve accuracy in polar case with (x,y) close to (0,0) --- src/projections/ortho.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index b4ecff7f..555f3f71 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -171,7 +171,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // ==> lam = atan2(xy.x, -xy.y * sign(phi0)) // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2 // rh^2 = cosphi^2 / (1 - es * sinphi^2) - // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) + // ==> cosphi^2 = rh^2 * (1 - es) / (1 - es * rh^2) const double rh2 = SQ(xy.x) + SQ(xy.y); if (rh2 >= 1. - 1e-15) { @@ -184,7 +184,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver lp.phi = 0; } else { - lp.phi = asin(sqrt((1 - rh2) / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); + lp.phi = acos(sqrt(rh2 * P->one_es / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); } lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1)); return lp; -- cgit v1.2.3 From 016e8e60747a2def2dfd7ffc7f6ad2e6aa8ba009 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Sep 2020 11:17:13 +0200 Subject: ortho.cpp: more precise reference to guidance note --- src/projections/ortho.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/projections') diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index 555f3f71..75b8199d 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -135,7 +135,7 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa PJ_XY xy; struct pj_opaque *Q = static_cast(P->opaque); - // From EPSG guidance note 7.2 + // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic const double cosphi = cos(lp.phi); const double sinphi = sin(lp.phi); const double coslam = cos(lp.lam); @@ -234,7 +234,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver return lp; } - // From EPSG guidance note 7.2 + // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic // It suggests as initial guess: // lp.lam = 0; -- cgit v1.2.3