1 files changed, 86 insertions, 68 deletions

diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index b2929a07..9f699fcf 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -12,7 +12,7 @@
  * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2016-06-09/2016-11-06
  *
  ******************************************************************************
- * Copyright (c) 2016, Thomas Knudsen/SDFE
+ * Copyright (c) 2016, 2017 Thomas Knudsen/SDFE 
  *
  * Permission is hereby granted, free of charge, to any person obtaining a
  * copy of this software and associated documentation files (the "Software"),
@@ -36,7 +36,7 @@
 #include <proj.h>
 #include "proj_internal.h"
 #include "projects.h"
-#include <geodesic.h>
+#include "geodesic.h"
 #include <stddef.h>
 #include <errno.h>
 
@@ -69,8 +69,8 @@ PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, dou
 
 
 
-/* Geodesic distance between two points with angular 2D coordinates */
-double proj_lp_dist (PJ *P, LP a, LP b) {
+/* Geodesic distance (in meter) between two points with angular 2D coordinates */
+double proj_lp_dist (const PJ *P, LP a, LP b) {
     double s12, azi1, azi2;
     /* Note: the geodesic code takes arguments in degrees */
     geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2);
@@ -90,9 +90,10 @@ double proj_xyz_dist (XYZ a, XYZ b) {
 
 
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
-double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs) {
+double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) {
     int i;
-    PJ_OBS o, u;
+    PJ_COORD o, u;
+    enum pj_io_units unit;
 
     if (0==P)
         return HUGE_VAL;
@@ -102,30 +103,48 @@ double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs)
         return HUGE_VAL;
     }
 
-    o.coo = obs.coo;
-
-    for (i = 0;  i < n;  i++) {
-        switch (direction) {
-            case PJ_FWD:
-                u  =  pj_fwdobs (o, P);
-                o  =  pj_invobs (u, P);
-                break;
-            case PJ_INV:
-                u  =  pj_invobs (o, P);
-                o  =  pj_fwdobs (u, P);
-                break;
-            default:
-                proj_errno_set (P, EINVAL);
-                return HUGE_VAL;
-        }
+    o = coo;
+
+    switch (direction) {
+        case PJ_FWD:
+            for (i = 0;  i < n;  i++) {
+                u  =  pj_fwdcoord (o, P);
+                o  =  pj_invcoord (u, P);
+            }
+            break;
+        case PJ_INV:
+            for (i = 0;  i < n;  i++) {
+                u  =  pj_invcoord (o, P);
+                o  =  pj_fwdcoord (u, P);
+            }
+            break;
+        default:
+            proj_errno_set (P, EINVAL);
+            return HUGE_VAL;
     }
 
-    return proj_xyz_dist (o.coo.xyz, obs.coo.xyz);
+    /* left when forward, because we do a roundtrip, and end where we begin */
+    unit = direction==PJ_FWD?  P->left:  P->right;
+    if (unit==PJ_IO_UNITS_RADIANS)
+        return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z);
+    
+    return proj_xyz_dist (coo.xyz, coo.xyz);
 }
 
 
+
+
+
+
+
+
+
+
+
+
+
 /* Apply the transformation P to the coordinate coo */
-PJ_OBS proj_trans_obs (PJ *P, enum proj_direction direction, PJ_OBS obs) {
+PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs) {
     if (0==P)
         return obs;
 
@@ -147,7 +166,7 @@ PJ_OBS proj_trans_obs (PJ *P, enum proj_direction direction, PJ_OBS obs) {
 
 
 /* Apply the transformation P to the coordinate coo */
-PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coo) {
+PJ_COORD proj_trans_coord (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
     if (0==P)
         return coo;
 
@@ -171,7 +190,7 @@ PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coo) {
 /*************************************************************************************/
 size_t proj_transform (
     PJ *P,
-    enum proj_direction direction,
+    PJ_DIRECTION direction,
     double *x, size_t sx, size_t nx,
     double *y, size_t sy, size_t ny,
     double *z, size_t sz, size_t nz,
@@ -319,7 +338,7 @@ size_t proj_transform (
 }
 
 /*****************************************************************************/
-int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs) {
+int proj_transform_obs (PJ *P, PJ_DIRECTION direction, size_t n, PJ_OBS *obs) {
 /******************************************************************************
     Batch transform an array of PJ_OBS.
 
@@ -337,7 +356,7 @@ int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *
 }
 
 /*****************************************************************************/
-int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord) {
+int proj_transform_coord (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord) {
 /******************************************************************************
     Batch transform an array of PJ_COORD.
 
@@ -407,34 +426,20 @@ PJ *proj_destroy (PJ *P) {
     return 0;
 }
 
-/* For now, if PJ itself is clean, we return the thread local error level. */
-/* This may change as OBS_API error reporting matures */
 int proj_errno (PJ *P) {
-    if (0==P)
-        return pj_ctx_get_errno (pj_get_default_ctx ());
-    if (0 != P->last_errno)
-        return P->last_errno;
     return pj_ctx_get_errno (pj_get_ctx (P));
 }
 
 /*****************************************************************************/
 void proj_errno_set (PJ *P, int err) {
 /******************************************************************************
-    Sets errno in the PJ, and bubbles it up to the context and pj_errno levels
-    through the low level pj_ctx interface.
+    Sets errno at the context and bubble it up to the thread local errno
 ******************************************************************************/
-    if (0==P) {
-        errno = EINVAL;
-        return;
-    }
-
     /* Use proj_errno_reset to explicitly clear the error status */
     if (0==err)
         return;
 
-    /* set local error level */
-    P->last_errno = err;
-    /* and let it bubble up */
+    /* For P==0 err goes to the default context */
     proj_context_errno_set (pj_get_ctx (P), err);
     errno = err;
     return;
@@ -454,14 +459,16 @@ void proj_errno_restore (PJ *P, int err) {
 
     See usage example under proj_errno_reset ()
 ******************************************************************************/
+    if (0==err)
+        return;
     proj_errno_set (P, err);
 }
 
 /*****************************************************************************/
 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.
+    Clears errno in the context and thread local levels
+    through the low level pj_ctx interface.
 
     Returns the previous value of the errno, for convenient reset/restore
     operations:
@@ -471,24 +478,17 @@ int proj_errno_reset (PJ *P) {
 
         do_something_with_P (P);
 
-        (* failure - keep latest error status *)
+        // failure - keep latest error status
         if (proj_errno(P))
             return;
-        (* success - restore previous error status *)
+        // success - restore previous error status
         proj_errno_restore (P, last_errno);
         return;
     }
 ******************************************************************************/
     int last_errno;
-    if (0==P) {
-        errno = EINVAL;
-        return EINVAL;
-    }
     last_errno = proj_errno (P);
 
-    /* set local error level */
-    P->last_errno = 0;
-    /* and let it bubble up */
     pj_ctx_set_errno (pj_get_ctx (P), 0);
     errno = 0;
     return last_errno;
@@ -501,15 +501,16 @@ PJ_CONTEXT *proj_context_create (void) {
 }
 
 
-void proj_context_destroy (PJ_CONTEXT *ctx) {
+PJ_CONTEXT *proj_context_destroy (PJ_CONTEXT *ctx) {
     if (0==ctx)
-        return;
+        return 0;
 
     /* Trying to free the default context is a no-op (since it is statically allocated) */
     if (pj_get_default_ctx ()==ctx)
-        return;
+        return 0;
 
     pj_ctx_free (ctx);
+    return 0;
 }
 
 
@@ -582,7 +583,7 @@ PJ_INFO proj_info(void) {
 
 
 /*****************************************************************************/
-PJ_PROJ_INFO proj_pj_info(const PJ *P) {
+PJ_PROJ_INFO proj_pj_info(PJ *P) {
 /******************************************************************************
     Basic info about a particular instance of a projection object.
 
@@ -594,6 +595,11 @@ PJ_PROJ_INFO proj_pj_info(const PJ *P) {
 
     memset(&info, 0, sizeof(PJ_PROJ_INFO));
 
+    /* Expected accuracy of the transformation. Hardcoded for now, will be improved */
+    /* later. Most likely to be used when a transformation is set up with           */
+    /* proj_create_crs_to_crs in a future version that leverages the EPSG database. */
+    info.accuracy = -1.0;
+
     if (!P) {
         return info;
     }
@@ -606,7 +612,7 @@ PJ_PROJ_INFO proj_pj_info(const PJ *P) {
     pj_strlcpy(info.description, P->descr, sizeof(info.description));
 
     /* projection definition */
-    def = pj_get_def((PJ *)P, 0); /* pj_get_def takes a non-const PJ pointer */
+    def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */
     pj_strlcpy(info.definition, &def[1], sizeof(info.definition)); /* def includes a leading space */
     pj_dealloc(def);
 
@@ -729,7 +735,7 @@ PJ_INIT_INFO proj_init_info(const char *initname){
 
 
 /*****************************************************************************/
-PJ_DERIVS proj_derivatives(const PJ *P, const LP lp) {
+PJ_DERIVS proj_derivatives(PJ *P, const LP lp) {
 /******************************************************************************
     Derivatives of coordinates.
 
@@ -739,8 +745,7 @@ PJ_DERIVS proj_derivatives(const PJ *P, const LP lp) {
 ******************************************************************************/
     PJ_DERIVS derivs;
 
-    /* casting to struct DERIVS for compatibility reasons */
-    if (pj_deriv(lp, 1e-5, (PJ *)P, (struct DERIVS *)&derivs)) {
+    if (pj_deriv(lp, 1e-5, P, &derivs)) {
         /* errno set in pj_derivs */
         memset(&derivs, 0, sizeof(PJ_DERIVS));
     }
@@ -750,7 +755,7 @@ PJ_DERIVS proj_derivatives(const PJ *P, const LP lp) {
 
 
 /*****************************************************************************/
-PJ_FACTORS proj_factors(const PJ *P, const LP lp) {
+PJ_FACTORS proj_factors(PJ *P, const LP lp) {
 /******************************************************************************
     Cartographic characteristics at point lp.
 
@@ -766,8 +771,7 @@ PJ_FACTORS proj_factors(const PJ *P, const LP lp) {
     /* pj_factors rely code being zero */
     factors.code = 0;
 
-    /* casting to struct FACTORS for compatibility reasons */
-    if (pj_factors(lp, (PJ *)P, 0.0, (struct FACTORS *)&factors)) {
+    if (pj_factors(lp, P, 0.0, &factors)) {
         /* errno set in pj_factors */
         memset(&factors, 0, sizeof(PJ_FACTORS));
     }
@@ -776,6 +780,22 @@ PJ_FACTORS proj_factors(const PJ *P, const LP lp) {
 }
 
 
+const PJ_ELLPS *proj_list_ellps(void) {
+    return pj_get_ellps_ref();
+}
+
+const PJ_UNITS *proj_list_units(void) {
+    return pj_get_units_ref();
+}
+
+const PJ_OPERATIONS *proj_list_operations(void) {
+    return pj_get_list_ref();
+}
+
+const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) {
+    return pj_get_prime_meridians_ref();
+}
+
 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);}
 
@@ -787,5 +807,3 @@ double proj_dmstor(const char *is, char **rs) {
 char*  proj_rtodms(char *s, double r, int pos, int neg) {
     return rtodms(s, r, pos, neg);
 }
-
-