OSGeo · kbevers · Jul 13, 2017 · Jul 12, 2017 · Jul 12, 2017 · Jul 12, 2017
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
@@ -235,12 +235,14 @@ int pj_cart_selftest (void) {
     PJ_CONTEXT *ctx;
     PJ *P;
     PJ_OBS a, b, obs[2];
+    PJ_COORD coord[2];
     int err;
     size_t n, sz;
     double dist, h, t;
     char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
     char *arg = {" +proj=utm +zone=32 +ellps=GRS80"};
     char *arg_def;
+    char buf[40];
 
     /* An utm projection on the GRS80 ellipsoid */
     P = proj_create (0, arg);
@@ -419,9 +421,91 @@ int pj_cart_selftest (void) {
     if (50            != obs[1].coo.lpz.z)    return 27; /* NOTE: unchanged */
     if (50==h) return 28;
 
-    /* Clean up */
+    /* test proj_transform_obs() */
+
+    obs[0].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0);
+    obs[1].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0);
+
+    if (proj_transform_obs(P, PJ_FWD, 2, obs))
+        return 30;
+
+    if (a.coo.lpz.lam != obs[0].coo.lpz.lam)  return 31;
+    if (a.coo.lpz.phi != obs[0].coo.lpz.phi)  return 32;
+    if (a.coo.lpz.z   != obs[0].coo.lpz.z)    return 33;
+    if (b.coo.lpz.lam != obs[1].coo.lpz.lam)  return 34;
+    if (b.coo.lpz.phi != obs[1].coo.lpz.phi)  return 35;
+    if (b.coo.lpz.z   != obs[1].coo.lpz.z)    return 36;
+
+    /* test proj_transform_coord() */
+
+    coord[0] = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0);
+    coord[1] = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0);
+    if (proj_transform_coord(P, PJ_FWD, 2, coord))
+        return 40;
+
+    if (a.coo.lpz.lam != coord[0].lpz.lam)  return 41;
+    if (a.coo.lpz.phi != coord[0].lpz.phi)  return 42;
+    if (a.coo.lpz.z   != coord[0].lpz.z)    return 43;
+    if (b.coo.lpz.lam != coord[1].lpz.lam)  return 44;
+    if (b.coo.lpz.phi != coord[1].lpz.phi)  return 45;
+    if (b.coo.lpz.z   != coord[1].lpz.z)    return 46;
+
+    /* Clean up  after transform* tests */
     proj_destroy (P);
 
+    /* test proj_create_crs_to_crs() */
+    P = proj_create_crs_to_crs(0, "epsg:25832", "epsg:25833");
+    if (P==0)
+        return 50;
+
+    a.coo.xy.x =  700000.0;
+    a.coo.xy.y = 6000000.0;
+    b.coo.xy.x =  307788.8761171057;
+    b.coo.xy.y = 5999669.3036037628;
+
+    a = proj_trans_obs(P, PJ_FWD, a);
+    if (dist > 1e-7)
+        return 51;
+    proj_destroy(P);
+
+    /* let's make sure that only entries in init-files results in a usable PJ */
+    P = proj_create_crs_to_crs(0, "proj=utm +zone=32 +datum=WGS84", "proj=utm +zone=33 +datum=WGS84");
+    if (P != 0) {
+        proj_destroy(P);
+        return 52;
+    }
+    proj_destroy(P);
+
+    /* Test proj_has_inverse() */
+    P = proj_create(0, "+proj=august"); /* august has no inverse */
+    if (proj_has_inverse(P)) {
+        proj_destroy(P);
+        return 60;
+    }
+    proj_destroy(P);
+
+    P = proj_create(0, "+proj=merc"); /* merc has an inverse */
+    if (!proj_has_inverse(P)) {
+        proj_destroy(P);
+        return 61;
+    }
+    proj_destroy(P);
+
+    /* test proj_rtodms() and proj_dmstor() */
+    if (strcmp("180dN", proj_rtodms(buf, M_PI, 'N', 'S')))
+        return 70;
+
+    if (proj_dmstor(&buf[0], NULL) != M_PI)
+        return 71;
+
+    if (strcmp("114d35'29.612\"S", proj_rtodms(buf, -2.0, 'N', 'S')))
+        return 72;
+
+    /* we can't expect perfect numerical accuracy so testing with a tolerance */
+    if (fabs(-2.0 - proj_dmstor(&buf[0], NULL)) > 1e-7)
+        return 73;
+
+
     return 0;
 }
 #endif
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
@@ -42,6 +42,7 @@
 #include <errno.h>
 
 
+/* Initialize PJ_COORD struct */
 PJ_COORD proj_coord (double x, double y, double z, double t) {
     PJ_COORD res;
     res.v[0] = x;
@@ -51,6 +52,24 @@ PJ_COORD proj_coord (double x, double y, double z, double t) {
     return res;
 }
 
+/* Initialize PJ_OBS struct */
+PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags) {
+    PJ_OBS res;
+    res.coo.v[0] = x;
+    res.coo.v[1] = y;
+    res.coo.v[2] = z;
+    res.coo.v[3] = t;
+    res.anc.v[0] = o;
+    res.anc.v[1] = p;
+    res.anc.v[2] = k;
+    res.id       = id;
+    res.flags    = flags;
+
+    return res;
+}
+
+
+
 /* Geodesic distance between two points with angular 2D coordinates */
 double proj_lp_dist (PJ *P, LP a, LP b) {
     double s12, azi1, azi2;
@@ -197,7 +216,7 @@ size_t proj_transform (
     In most cases, the stride will be identical for x, y,z, and t, since they will
     typically be either individual arrays (stride = sizeof(double)), or strided
     views into an array of application specific data structures (stride = sizeof (...)).
-    
+
     But in order to support cases where x, y, z, and t come from heterogeneous
     sources, individual strides, sx, sy, sz, st, are used.
 
@@ -225,7 +244,7 @@ size_t proj_transform (
     if (0==ny) y = &null_broadcast;
     if (0==nz) z = &null_broadcast;
     if (0==nt) t = &null_broadcast;
-    
+
     /* nothing to do? */
     if (0==nx+ny+nz+nt)
         return 0;
@@ -296,10 +315,46 @@ size_t proj_transform (
         *z = coord.xyzt.z;
     if (nt==1)
         *t = coord.xyzt.t;
-    
+
     return i;
 }
 
+/*****************************************************************************/
+int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs) {
+/******************************************************************************
+    Batch transform an array of PJ_OBS.
+
+    Returns 0 if all observations are transformed without error, otherwise
+    returns error number.
+******************************************************************************/
+    size_t i;
+    for (i=0; i<n; i++) {
+        obs[i] = proj_trans_obs(P, direction, obs[i]);
+        if (proj_errno(P))
+            return proj_errno(P);
+    }
+
+    return 0;
+}
+
+/*****************************************************************************/
+int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord) {
+/******************************************************************************
+    Batch transform an array of PJ_COORD.
+
+    Returns 0 if all coordinates are transformed without error, otherwise
+    returns error number.
+******************************************************************************/
+    size_t i;
+    for (i=0; i<n; i++) {
+        coord[i] = proj_trans_coord(P, direction, coord[i]);
+        if (proj_errno(P))
+            return proj_errno(P);
+    }
+
+    return 0;
+}
+
 
 
 PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
@@ -314,6 +369,36 @@ PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv) {
     return pj_init_ctx (ctx, argc, argv);
 }
 
+/*****************************************************************************/
+PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *srid_from, const char *srid_to) {
+/******************************************************************************
+    Create a transformation pipeline between two known coordinate reference
+    systems.
+
+    srid_from and srid_to should be the value part of a +init=... parameter
+    set, i.e. "epsg:25833" or "IGNF:AMST63". Any projection definition that is
+    can be found in a init-file in PROJ_LIB is a valid input to this function.
+
+    For now the function mimics the cs2cs app: An input and an output CRS is
+    given and coordinates are transformed via a hub datum (WGS84). This
+    transformation strategy is referred to as "early-binding" by the EPSG. The
+    function can be extended to support "late-binding" transformations in the
+    future without affecting users of the function.
+
+    Example call:
+
+        PJ *P = proj_create_crs_to_crs(0, "epsg:25832", "epsg:25833");
+
+******************************************************************************/
+    PJ *P;
+    char buffer[256];
+
+    sprintf(buffer, "+proj=pipeline +step +init=%s +inv +step +init=%s", srid_from, srid_to);
+    P = proj_create(ctx, buffer);
+
+    return P;
+}
+
 PJ *proj_destroy (PJ *P) {
     pj_free (P);
     return 0;
@@ -374,13 +459,13 @@ int proj_errno_reset (PJ *P) {
 /******************************************************************************
     Clears errno in the PJ, and bubbles it up to the context and
     pj_errno levels through the low level pj_ctx interface.
-    
+
     Returns the previous value of the errno, for convenient reset/restore
     operations:
 
     void foo (PJ *P) {
         int last_errno = proj_errno_reset (P);
-        
+
         do_something_with_P (P);
 
         (* failure - keep latest error status *)
@@ -441,10 +526,14 @@ void *proj_release (void *buffer) {
 double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
 double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
 
+int proj_has_inverse(PJ *P) {
+    return (P->inv != 0 || P->inv3d != 0 || P->invobs != 0);
+}
 
-/* The shape of jazz to come! */
+double proj_dmstor(const char *is, char **rs) {
+    return dmstor(is, rs);
+}
 
-int proj_transform_obs   (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs);
-int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord);
-PJ_OBS   proj_obs   (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags);
-PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *def_from, const char *def_to);
+char*  proj_rtodms(char *s, double r, int pos, int neg) {
+    return rtodms(s, r, pos, neg);
+}
diff --git a/src/proj.def b/src/proj.def
@@ -92,42 +92,48 @@ EXPORTS
 
     proj_create                     @90
     proj_create_argv                @91
-    proj_destroy                    @92
+    proj_create_crs_to_crs          @92
+    proj_destroy                    @93
 
-    proj_trans_obs                  @93
-    proj_trans_coord                @94
-    proj_transform                  @95
-    proj_roundtrip                  @96
+    proj_trans_obs                  @94
+    proj_trans_coord                @95
+    proj_transform                  @96
+    proj_transform_obs              @97
+    proj_transform_coord            @98
+    proj_roundtrip                  @99
 
-    proj_coord                      @97
-    proj_coord_error                @98
-    proj_obs_error                  @99
+    proj_coord                     @100
+    proj_obs                       @101
+    proj_coord_error               @102
+    proj_obs_error                 @103
 
-    proj_errno                     @100
-    proj_errno_set                 @101
-    proj_errno_reset               @102
-    proj_errno_restore             @103
-    proj_context_errno_set         @104
+    proj_errno                     @104
+    proj_errno_set                 @105
+    proj_errno_reset               @106
+    proj_errno_restore             @107
+    proj_context_errno_set         @108
 
-    proj_context_create            @105
-    proj_context_set               @106
-    proj_context_inherit           @107
-    proj_context_destroy           @108
+    proj_context_create            @109
+    proj_context_set               @110
+    proj_context_inherit           @111
+    proj_context_destroy           @112
 
-    proj_lp_dist                   @109
-    proj_xy_dist                   @110
-    proj_xyz_dist                  @111
+    proj_lp_dist                   @113
+    proj_xy_dist                   @114
+    proj_xyz_dist                  @115
 
-    proj_log_level                 @112
-    proj_log_func                  @113
-    proj_log_error                 @114
-    proj_log_debug                 @115
-    proj_log_trace                 @116
+    proj_log_level                 @116
+    proj_log_func                  @117
+    proj_log_error                 @118
+    proj_log_debug                 @119
+    proj_log_trace                 @120
 
-    proj_definition_retrieve       @117
-    proj_release                   @118
-    proj_torad                     @119
-    proj_todeg                     @120
+    proj_definition_retrieve       @121
+    proj_release                   @122
+    proj_torad                     @123
+    proj_todeg                     @124
+    proj_has_inverse               @125
+    proj_rtodms                    @126
+    proj_dmstor                    @127
 
-
-    pj_find_file                   @121
+    pj_find_file                   @128