From 5cad2f68f2a68c8c6854d0bfc62700299598dfa1 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 27 Mar 2019 18:55:35 +0100 Subject: ocea: fix behaviour when +alpha is specified The one-point case was completely broken with lat_0 being ignored. I've fixed that, and also modified the alpha orientation, so that the angle corresponds to the one of the 2 point case when going from point 1 to point 2, similarly to what omerc does. This was found rather experimentaly with the added test cases that try to find equivalence between 1-point and 2-point cases. Fixes #1379 and adresses https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=13930 --- src/projections/ocea.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/projections/ocea.cpp b/src/projections/ocea.cpp index 4e28f727..3141dd11 100644 --- a/src/projections/ocea.cpp +++ b/src/projections/ocea.cpp @@ -50,7 +50,7 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ PJ *PROJECTION(ocea) { - double phi_0=0.0, phi_1, phi_2, lam_1, lam_2, lonz, alpha; + double phi_1, phi_2, lam_1, lam_2, lonz, alpha; struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); if (nullptr==Q) @@ -63,12 +63,17 @@ PJ *PROJECTION(ocea) { /*If the keyword "alpha" is found in the sentence then use 1point+1azimuth*/ if ( pj_param(P->ctx, P->params, "talpha").i) { /*Define Pole of oblique transformation from 1 point & 1 azimuth*/ - alpha = pj_param(P->ctx, P->params, "ralpha").f; + // ERO: I've added M_PI so that the alpha is the angle from point 1 to point 2 + // from the North in a clockwise direction + // (to be consistent with omerc behaviour) + alpha = M_PI + pj_param(P->ctx, P->params, "ralpha").f; lonz = pj_param(P->ctx, P->params, "rlonc").f; /*Equation 9-8 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - lam_p = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz; + // Actually slightliy modified to use atan2(), as it is suggested by + // Snyder for equation 9-1, but this is not mentionned here + lam_p = atan2(-cos(alpha) , -sin(P->phi0) * sin(alpha)) + lonz; /*Equation 9-7 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/ - phi_p = asin(cos(phi_0) * sin(alpha)); + phi_p = asin(cos(P->phi0) * sin(alpha)); /*If the keyword "alpha" is NOT found in the sentence then use 2points*/ } else { /*Define Pole of oblique transformation from 2 points*/ -- cgit v1.2.3