Added horizontal gridshift driver.

Until now gridshifts has not been working with the new API in
proj.h since parsing of +nadgrids and +geoidgrids is build
into pj_transform(). This commit introduces the possibility to
do horizontal gridshift with the pipeline API. The hgridshift
driver is a simple wrapper for the funtions in pj_apply_gridshift.c
that is also used by pj_transform().

src/Makefile.am

@@ -75,7 +75,7 @@ libproj_la_SOURCES = \
 	pj_strtod.c \
 	\
 	pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \
-	PJ_vgridshift.c
+	PJ_vgridshift.c PJ_hgridshift.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)

src/PJ_hgridshift.c

@@ -0,0 +1,131 @@
+#define PJ_LIB__
+#include <proj.h>
+#include <projects.h>
+
+PROJ_HEAD(hgridshift, "Horizontal grid shift");
+
+static void *freeup_msg (PJ *P, int errlev) {
+    if (0==P)
+        return 0;
+
+    if (0!=P->ctx)
+        pj_ctx_set_errno (P->ctx, errlev);
+
+    return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+    freeup_msg (P, 0);
+    return;
+}
+
+
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+    PJ_TRIPLET point;
+    point.lpz = lpz;
+
+    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 );
+    }
+
+    return point.xyz;
+}
+
+
+static LPZ reverse_3d(XYZ xyz, PJ *P) {
+    PJ_TRIPLET point;
+    point.xyz = xyz;
+
+    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 );
+    }
+
+    return point.lpz;
+}
+
+
+static PJ_OBS forward_obs(PJ_OBS obs, PJ *P) {
+    PJ_OBS point;
+    point.coo.xyz = forward_3d (obs.coo.lpz, P);
+    return point;
+}
+
+
+static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
+    PJ_OBS point;
+    point.coo.lpz = reverse_3d (obs.coo.xyz, P);
+    return point;
+}
+
+
+PJ *PROJECTION(hgridshift) {
+
+    if (!pj_param(P->ctx, P->params, "tgrids").i) {
+        pj_log_error(P, "hgridshift: +grids parameter missing.");
+        return freeup_msg(P, -1);
+    }
+
+    /* 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) );
+
+    /* Was gridlist compiled properly? */
+    if ( pj_ctx_get_errno(P->ctx) ) {
+        pj_log_error(P, "hgridshift: could not find required grid(s).");
+        return freeup_msg(P, -38);
+    }
+
+    P->fwdobs = forward_obs;
+    P->invobs = reverse_obs;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = 0;
+    P->inv    = 0;
+
+    P->left  = PJ_IO_UNITS_RADIANS;
+    P->right = PJ_IO_UNITS_RADIANS;
+
+    return P;
+}
+
+
+#ifndef PJ_SELFTEST
+/* selftest stub */
+int pj_hgridshift_selftest (void) {return 0;}
+#else
+int pj_hgridshift_selftest (void) {
+    PJ *P;
+    PJ_OBS a;
+    double dist;
+
+    /* Since the nzgd2kgrid0005.gsb grid is not part of the
+     * basic PROJ.4 distribution we only include it as an
+     * optional grid in the test. The test is fairly useless
+     * if the grid isn't found, but at least it won't fail
+     * because of a missing file */
+    P = pj_create ("+proj=hgridshift +grids=@nzgd2kgrid0005.gsb +ellps=GRS80");
+    if (0==P)
+        return 10;
+
+    a = pj_obs_null;
+    a.coo.lpz.lam = TORAD(173);
+    a.coo.lpz.phi = TORAD(-45);
+
+    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    if (dist > 0.00000001)
+        return 1;
+
+    pj_free(P);
+
+    return 0;
+}
+#endif

src/lib_proj.cmake

@@ -78,6 +78,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_hammer.c
         PJ_hatano.c
         PJ_helmert.c
+        PJ_hgridshift.c
         PJ_horner.c
         PJ_igh.c
         PJ_isea.c

src/makefile.vc

@@ -60,7 +60,7 @@ support = \
 
 pipeline = \
 	pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
-	PJ_vgridshift.obj
+	PJ_vgridshift.obj PJ_hgridshift.obj
 
 geodesic = geodesic.obj
 

src/pj_list.h

@@ -53,6 +53,7 @@ PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area")
 PROJ_HEAD(healpix, "HEALPix")
 PROJ_HEAD(rhealpix, "rHEALPix")
 PROJ_HEAD(helmert,   "3- and 7-parameter Helmert shift")
+PROJ_HEAD(hgridshift, "Horizontal grid shift")
 PROJ_HEAD(horner,    "Horner polynomial evaluation")
 PROJ_HEAD(igh,  "Interrupted Goode Homolosine")
 PROJ_HEAD(imw_p, "International Map of the World Polyconic")