Merge remote-tracking branch 'origin/master'

7 files changed, 33 insertions, 16 deletions

diff --git a/src/transformations/affine.cpp b/src/transformations/affine.cpp
index bda54f1e..28f73b9a 100644
--- a/src/transformations/affine.cpp
+++ b/src/transformations/affine.cpp
@@ -25,7 +25,6 @@
 #include <errno.h>
 #include <math.h>
 
-#include "proj_internal.h"
 #include "proj.h"
 #include "proj_internal.h"
 
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index c28e1489..5bb86909 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -56,7 +56,6 @@ grid-values in units of mm/year in ENU-space.
 #include "proj.h"
 #include "proj_internal.h"
 #include "proj_math.h"
-#include "proj_internal.h"
 
 PROJ_HEAD(deformation, "Kinematic grid shift");
 
@@ -92,7 +91,7 @@ static PJ_XYZ get_grid_shift(PJ* P, PJ_XYZ cartesian) {
 
     /* look up correction values in grids */
     shift.lp    = proj_hgrid_value(P, geodetic.lp);
-    shift.enu.u = proj_vgrid_value(P, geodetic.lp);
+    shift.enu.u = proj_vgrid_value(P, geodetic.lp, 1.0);
 
     if (proj_errno(P) == PJD_ERR_GRID_AREA)
         proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp
index 034f76f4..63785ea5 100644
--- a/src/transformations/helmert.cpp
+++ b/src/transformations/helmert.cpp
@@ -53,7 +53,6 @@ Last update: 2018-10-26
 #include <math.h>
 
 #include "proj_internal.h"
-#include "proj_internal.h"
 #include "geocent.h"
 
 PROJ_HEAD(helmert, "3(6)-, 4(8)- and 7(14)-parameter Helmert shift");
@@ -567,14 +566,14 @@ PJ *TRANSFORMATION(helmert, 0) {
     if (pj_param_exists (P->params, "theta")) {
         P->left  = PJ_IO_UNITS_PROJECTED;
         P->right = PJ_IO_UNITS_PROJECTED;
+        P->fwd    = helmert_forward;
+        P->inv    = helmert_reverse;
     }
 
     P->fwd4d  = helmert_forward_4d;
     P->inv4d  = helmert_reverse_4d;
     P->fwd3d  = helmert_forward_3d;
     P->inv3d  = helmert_reverse_3d;
-    P->fwd    = helmert_forward;
-    P->inv    = helmert_reverse;
 
     Q = (struct pj_opaque_helmert *)P->opaque;
 
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 2e2294cb..90633939 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -6,7 +6,6 @@
 #include <time.h>
 
 #include "proj_internal.h"
-#include "proj_internal.h"
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp
index f5d749c4..a6638773 100644
--- a/src/transformations/horner.cpp
+++ b/src/transformations/horner.cpp
@@ -85,7 +85,6 @@
 
 #include "proj.h"
 #include "proj_internal.h"
-#include "proj_internal.h"
 
 PROJ_HEAD(horner, "Horner polynomial evaluation");
 
diff --git a/src/transformations/molodensky.cpp b/src/transformations/molodensky.cpp
index c389fd32..7d17f64c 100644
--- a/src/transformations/molodensky.cpp
+++ b/src/transformations/molodensky.cpp
@@ -49,7 +49,6 @@
 
 #include "proj.h"
 #include "proj_internal.h"
-#include "proj_internal.h"
 
 PROJ_HEAD(molodensky, "Molodensky transform");
 
@@ -140,16 +139,26 @@ static PJ_LPZ calc_standard_params(PJ_LPZ lpz, PJ *P) {
 
     /* ellipsoid radii of curvature */
     double rho = RM(a, P->es, lpz.phi);
-    double nu  = RN(a, P->e2s, lpz.phi);
+    double nu  = RN(a, P->es, lpz.phi);
 
     /* delta phi */
     dphi  = (-dx*sphi*clam) - (dy*sphi*slam) + (dz*cphi)
             + ((nu * P->es * sphi * cphi * da) / a)
             + (sphi*cphi * ( rho/(1-f) + nu*(1-f))*df);
-    dphi /= (rho + lpz.z);
+    const double dphi_denom = rho + lpz.z;
+    if( dphi_denom == 0.0 ) {
+        lpz.lam = HUGE_VAL;
+        return lpz;
+    }
+    dphi /= dphi_denom;
 
     /* delta lambda */
-    dlam = (-dx*slam + dy*clam) / ((nu+lpz.z)*cphi);
+    const double dlam_denom = (nu+lpz.z)*cphi;
+    if( dlam_denom == 0.0 ) {
+        lpz.lam = HUGE_VAL;
+        return lpz;
+    }
+    dlam = (-dx*slam + dy*clam) / dlam_denom;
 
     /* delta h */
     dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - (a/nu)*da + nu*(1-f)*sphi*sphi*df;
@@ -183,7 +192,12 @@ static PJ_LPZ calc_abridged_params(PJ_LPZ lpz, PJ *P) {
 
     /* delta lambda */
     dlam = -dx*slam + dy*clam;
-    dlam /= RN(P->a, P->e2s, lpz.phi)*cphi;
+    const double dlam_denom = RN(P->a, P->es, lpz.phi)*cphi;
+    if( dlam_denom == 0.0 ) {
+        lpz.lam = HUGE_VAL;
+        return lpz;
+    }
+    dlam /= dlam_denom;
 
     /* delta h */
     dh = dx*cphi*clam + dy*cphi*slam + dz*sphi - da + adffda*sphi*sphi;
@@ -230,6 +244,10 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     } else {
         lpz = calc_standard_params(lpz, P);
     }
+    if( lpz.lam == HUGE_VAL ) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return proj_coord_error().xyz;
+    }
 
     /* offset coordinate */
     point.lpz.phi += lpz.phi;
@@ -258,6 +276,11 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     else
         lpz = calc_standard_params(point.lpz, P);
 
+    if( lpz.lam == HUGE_VAL ) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return proj_coord_error().lpz;
+    }
+
     /* offset coordinate */
     point.lpz.phi -= lpz.phi;
     point.lpz.lam -= lpz.lam;
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index fda38ec3..de0cdd8c 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -6,7 +6,6 @@
 #include <time.h>
 
 #include "proj_internal.h"
-#include "proj_internal.h"
 
 PROJ_HEAD(vgridshift, "Vertical grid shift");
 
@@ -26,7 +25,7 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     if (P->vgridlist_geoid != nullptr) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z += Q->forward_multiplier * proj_vgrid_value(P, point.lp);
+        point.xyz.z += proj_vgrid_value(P, point.lp, Q->forward_multiplier);
     }
 
     return point.xyz;
@@ -41,7 +40,7 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     if (P->vgridlist_geoid != nullptr) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        point.xyz.z -= Q->forward_multiplier * proj_vgrid_value(P, point.lp);
+        point.xyz.z -= proj_vgrid_value(P, point.lp, Q->forward_multiplier);
     }
 
     return point.lpz;