Refactor grid shift functions.

This refactoring of the grid shift functions allows for easier access to
the actual grid values, as well as making it possible to perform
horizontal grid shift on a single coordinate without using the
pj_apply_gridshift* functions. The latter simplifies the execution path
of the forward and inverse functions in PJ_hgridshift.c.

This commit introduces proj_*grid_init, proj_*grid_value
and proj_hgrid_apply. The init functions initialises horizontal and
vertical grids respectivelive (wrappers for pj_gridlist_from_nadgrids
with simpler parameters). The proj_*grid_value functions returns the
specific grid value at coordinate lp. The proj_hgrid_apply function
applies the grid offset to the input coordinate and outputs the adjusted
coordinate.

1 files changed, 13 insertions, 39 deletions

diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 0adc9e00..659039ab 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -4,7 +4,6 @@
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
-
 static XYZ forward_3d(LPZ lpz, PJ *P) {
     PJ_TRIPLET point;
     point.lpz = lpz;
@@ -12,9 +11,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
     if (P->gridlist != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_gridshift_3( P->ctx, P->gridlist,
-                              P->gridlist_count, 1, 1, 0,
-                              &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD);
     }
 
     return point.xyz;
@@ -28,9 +25,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     if (P->gridlist != NULL) {
         /* Only try the gridshift if at least one grid is loaded,
          * otherwise just pass the coordinate through unchanged. */
-        pj_apply_gridshift_3( P->ctx, P->gridlist,
-                              P->gridlist_count, 0, 1, 0,
-                              &point.xyz.x, &point.xyz.y, &point.xyz.z );
+        point.lp = proj_hgrid_apply(P, point.lp, PJ_INV);
     }
 
     return point.lpz;
@@ -52,29 +47,8 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
 
 
 
-#if 0
-static XY forward_xy(LP lp, PJ *P) {
-    PJ_TRIPLET point;
-    point.lp = lp;
-    point.lpz.z = 0;
-    point.xyz = forward_3d (point.lpz, P);
-    return point.xy;
-}
-
-
-static LP reverse_lp(XY xy, PJ *P) {
-    PJ_TRIPLET point;
-    point.xy = xy;
-    point.xyz.z = 0;
-    point.lpz = reverse_3d (point.xyz, P);
-    return point.lp;
-}
-#endif
-
-
-
 PJ *PROJECTION(hgridshift) {
-    
+
     P->fwdobs = forward_obs;
     P->invobs = reverse_obs;
     P->fwd3d  = forward_3d;
@@ -90,14 +64,12 @@ PJ *PROJECTION(hgridshift) {
         return pj_default_destructor (P, PJD_ERR_NO_ARGS);
     }
 
-    /* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */
-    P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s,
-                                             &(P->gridlist_count) );
 
+    proj_hgrid_init(P, "grids");
     /* Was gridlist compiled properly? */
-    if ( pj_ctx_get_errno(pj_get_ctx(P)) ) {
+    if ( proj_errno(P) ) {
         proj_log_error(P, "hgridshift: could not find required grid(s).");
-        return pj_default_destructor (P, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
     }
 
     return P;
@@ -119,26 +91,28 @@ int pj_hgridshift_selftest (void) {
         proj_destroy (P);
         return 99;
     }
-    
+
     /* fail on purpose: open non-existing grid */
     P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb");
     if (0!=P) {
         proj_destroy (P);
         return 999;
     }
-        
+
     /* Failure most likely means the grid is missing */
     P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80");
     if (0==P)
         return 10;
-        
+
     a = proj_obs_null;
     a.coo.lpz.lam = PJ_TORAD(173);
     a.coo.lpz.phi = PJ_TORAD(-45);
-    
+
     dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
-    if (dist > 0.00000001)
+    if (dist > 0.00000001) {
+        printf("dist: %f\n",dist);
         return 1;
+    }
 
     expect.coo.lpz.lam = PJ_TORAD(172.999892181021551);
     expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613);