From 299135153e33869157ecd432dd0194743e959170 Mon Sep 17 00:00:00 2001
From: Kristian Evers <kristianevers@gmail.com>
Date: Fri, 13 Apr 2018 22:42:54 +0200
Subject: 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
---
 src/PJ_deformation.c     | 7 +++++++
 src/pj_apply_gridshift.c | 7 ++++++-
 2 files changed, 13 insertions(+), 1 deletion(-)

(limited to 'src')

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);
-- 
cgit v1.2.3