Fix segfault in deformation

When transforming coordinates outside the grid model the deformation
operation failed spectatularly. This is now fixed by checking that the
coordinate is inside the grid. If it isn't an error is returned and a
debugging log message is issued.

Closes #934

2 files changed, 13 insertions, 1 deletions

diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
index d6376e40..58c0e4df 100644
--- a/src/PJ_deformation.c
+++ b/src/PJ_deformation.c
@@ -82,6 +82,7 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
 ********************************************************************************/
     PJ_COORD geodetic, shift, temp;
     double sp, cp, sl, cl;
+    int previous_errno = proj_errno_reset(P);
 
     /* cartesian to geodetic */
     geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart);
@@ -90,6 +91,10 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
     shift.lp    = proj_hgrid_value(P, geodetic.lp);
     shift.enu.u = proj_vgrid_value(P, geodetic.lp);
 
+    if (proj_errno(P) == PJD_ERR_GRID_AREA)
+        proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model",
+                       proj_todeg(geodetic.lp.lam), proj_todeg(geodetic.lp.phi));
+
     /* grid values are stored as mm/yr, we need m/yr */
     shift.xyz.x /= 1000;
     shift.xyz.y /= 1000;
@@ -108,6 +113,8 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) {
 
     shift.xyz = temp.xyz;
 
+    proj_errno_restore(P, previous_errno);
+
     return shift.xyz;
 }
 
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index c503ec0b..1ef39b20 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -306,13 +306,18 @@ int proj_hgrid_init(PJ* P, const char *grids) {
 /********************************************/
 LP proj_hgrid_value(PJ *P, LP lp) {
     struct CTABLE *ct;
-    LP out;
+    LP out = proj_coord_error().lp;
 
     ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist);
+    if (ct == 0) {
+        pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA);
+        return out;
+    }
 
     /* normalize input to ll origin */
     lp.lam -= ct->ll.lam;
     lp.phi -= ct->ll.phi;
+
     lp.lam = adjlon(lp.lam - M_PI) + M_PI;
 
     out = nad_intr(lp, ct);