obs_api revision and improvements: new namespace etc. (#530)

* obs_api revision and improvements: new namespace etc.

* Minor corrections: use unused functions, add missing prototype, added fwd/invcoord to PJ

* Minor correction: MSVC barfs on va_arg type specification. Trying once more with added parens

* Reverting paren correction, which appears to be a non-solution

* Significant improvements to the OBS_API, plus a number of corrections, mostly in test code, to reflect API changes

* Added two missing prototypes

* Adding the proj_transform function and some related checks

* Improvements to proj_transform etc. based on suggestions from Even Rouault

* Reducing the libc include footprint of proj.h - again based on suggestions from Even Rouault

* A few minor corrections

* Eliminate a MSVC warning about non-initialized usage. Not an actual issue, as another check has locked the path, but at least this should calm down MSVC

* Improved support for the errno reset/restore paradigm

* Introduced the internal header proj_internal.h; Moved most non-API stuff from pj_obs_api.c to pj_internal.c

* Adding proj_internal.h to HEADERS_LIBPROJ to solve build problems

* Correcting a few pj...proj blunders in PJ_pipeline.c

* Correcting a few additional blunders in PJ_pipeline.c

* Changed angle-brackets to quotation marks in includes of proj_internal.h

* Some minor build system repairs

* Some PJ_CONTEXT usage simplifications following suggestions by Kristian Evers @kbevers

* Added version numbering to proj.h - Fixes #529

* remove proj_errno_restore macro, provide function implementation

* Add proj_get_definition. Fixes #538

* Added library specific deallocator proj_buffer_free, eliminating a potential cross heap issues on Windows. Thx to Even Rouault for spotting this

* Got rid of a meaningless cast in proj_buffer_free

* Added some missing functions to proj.def (again spotted by @rouault); removed some not-yet-implemented material from proj.h

* Renamed proj_get_definition/proj_buffer_free to proj_definition_create/proj_definition_destroy, for symmetry and clarity

* Renaming the definition handlers to proj_definition_retrieve / proj_free

* Renaming proj_free to proj_release

25 files changed, 1301 insertions, 781 deletions

diff --git a/examples/pj_obs_api_mini_demo.c b/examples/pj_obs_api_mini_demo.c
new file mode 100644
index 00000000..855c37f8
--- /dev/null
+++ b/examples/pj_obs_api_mini_demo.c
@@ -0,0 +1,66 @@
+/*******************************************************************************
+    Tiny test of an evolving new API, demonstrating simple examples of
+    2D and 3D transformations.
+
+    The main transformation setup object is PJ, well known from the two
+    former proj APIs (projects.h and proj_api.h)
+
+    The main data objects PJ_COORD and PJ_OBS are new, but encapsulate
+    the older LP, XY etc. objects in a framework for storing a 3D
+    observation taken at a 4D point in space-time, and including some
+    additional metadata (e.g. a serial number or an epsg code).
+
+    PJ_OBS and PJ_COORD use unions to enforce explicit statement of what
+    kind of coordinates are expected at a given spot in the code, where
+    the old API uses type punning, implicitly assuming that "you know what
+    you do". For backward compatibility, the new API is not really type
+    safe in the sense that you cannot use a cartesian coordinate where a
+    geographic is expected - but it makes it possible to explicitly state
+    in the code whet the programmer expected and intended.
+
+    The proj thread contexts have not seen widespread use, so one of the
+    intentions with this new API is to make them less visible on the API
+    surface.
+    
+    A series of experiments have, however, shown that they, for (mostly)
+    historical reasons, are very hard to eliminate totally. But we have
+    reduced their API surface presence to a constructor and a destructor,
+    plus an extra argument to the PJ constructor, pj_create().
+    
+    For single threaded programs, the calls to the context constructor
+    and destructor may be left out, and the default context selected
+    by passing a null-pointer to pj_create.
+
+    Thomas Knudsen, 2016-10-30/2017-07-06
+*******************************************************************************/
+#include <stdio.h>
+#include <proj.h>
+
+int main (void) {
+    PJ_CONTEXT *C;
+    PJ *P;
+    PJ_COORD a, b;
+    
+    /* or you may set C=0 if you are sure you will use PJ objects from only one thread */
+    C = proj_context_create();
+
+    P = proj_create (C, "+proj=utm +zone=32 +ellps=GRS80");
+    if (0==P)
+        return puts ("Oops"), 0;
+
+    /* a coordinate union representing Copenhagen: 55d N, 12d E    */
+    /* note: PROJ.4 works in radians, hence the proj_torad() calls */
+    a = proj_coord (proj_torad(12), proj_torad(55), 0, 0);
+
+    /* transform to UTM zone 32, then back to geographical */
+    /* note the use of union selectors to indicate what kind of coordinates are expected */
+    b = proj_trans_coord (P, PJ_FWD, a);
+    printf ("easting: %g, northing: %g\n", b.en.e, b.en.n);
+    b = proj_trans_coord (P, PJ_INV, b);
+    printf ("longitude: %g, latitude: %g\n", b.lp.lam, b.lp.phi);
+
+    /* Clean up */
+    proj_destroy (P);
+    proj_context_destroy (C); /* may be omitted in the single threaded case */
+    return 0;
+}
diff --git a/examples/pj_obs_api_test.c b/examples/pj_obs_api_test.c
deleted file mode 100644
index bb792c62..00000000
--- a/examples/pj_obs_api_test.c
+++ /dev/null
@@ -1,243 +0,0 @@
-/*******************************************************************************
-    Tiny test of an evolving new API, demonstrating simple examples of
-    2D and 3D transformations.
-
-    The main transformation setup object is PJ, well known from the two
-    former proj APIs (projects.h and proj_api.h)
-
-    The main data object PJ_OBS is new, but encapsulates the older
-    LP, XY etc. objects in a framework for storing a 3D observation taken at
-    a 4D point in space-time, and including some additional metadata (e.g.
-    a serial number or an epsg code).
-
-    PJ_OBS uses unions to enforce explicit statement of what kind of
-    coordinates are expected at a given spot in the code.
-
-    The primary elements of the API are:
-
-    pj_create (char *desc):
-        Create a new PJ object from a text description.
-    pj_trans (PJ *P, int direction, PJ_OBS obs):
-        Forward or inverse transformation of obs.
-    pj_error (PJ *P):
-        Check error status of P.
-    pj_free (PJ *P):
-        Clean up when done with the transformation object.
-
-    These 4 functions should cover most of the common use cases. Additional
-    functionality for handling thread contexts are provided by the functions
-    pj_debug_set,  pj_error_set,  pj_log_set,  pj_context_renew,
-    pj_context_inherit,    pj_context_free,    pj_fileapi_set.
-
-    The proj thread contexts have not seen widespread use, so one of the
-    intentions with this new API is to make them less visible on the API
-    surface: Contexts do not have a life by themselves, they are only
-    visible through their associated PJs, and the number of functions
-    supporting them is limited.
-
-    Compilation example:
-    gcc -pedantic -W -Wall -Wextra -I../src -o pj_proj_test pj_proj_test.c -L../../../build/proj.4/lib -lproj_4_9
-
-    Thomas Knudsen, 2016-10-30/11-03
-*******************************************************************************/
-#include <proj.h>
-
-int pj_pipeline_selftest (void);
-
-int main (void) {
-    PJ *P;
-    PJ_OBS a, b;
-    char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
-    int err;
-    double dist;
-    XY cph_utm32;
-
-    /* Log everything libproj offers to log for you */
-    pj_log_level (0, PJ_LOG_TRACE);
-
-    /* An utm projection on the GRS80 ellipsoid */
-    P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
-    if (0==P)
-        return puts ("Oops"), 0;
-
-    /* Clean up */
-    pj_free (P);
-
-    /* Same projection, now using argc/argv style initialization */
-    P = pj_create_argv (3, args);
-    if (0==P)
-        return puts ("Oops"), 0;
-
-    /* zero initialize everything, then set (longitude, latitude) to (12, 55) */
-    a = pj_obs_null ();
-    /* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(55);
-
-    /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);
-    cph_utm32 = b.coo.xy;
-
-    /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
-    printf ("INV:   %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
-
-    /* Null projection */
-    a = pj_trans (P, PJ_IDENT, a);
-    printf ("IDENT: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
-
-    /* Forward again, to get two linear items for comparison */
-    a = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);
-
-    dist = pj_xy_dist (a.coo.xy, b.coo.xy);
-    printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);
-
-    /* should be identical - checking whether null-init is done */
-    dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz);
-    printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);
-
-    /* Invalid projection */
-    a = pj_trans (P, 42, a);
-    if (a.coo.lpz.lam!=DBL_MAX)
-        printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
-    err = pj_err_level (P, PJ_ERR_TELL);
-    printf ("pj_err_level: %d\n", err);
-
-    /* Clean up */
-    pj_free (P);
-
-    /* Now do some 3D transformations */
-    P = pj_create ("+proj=cart +ellps=GRS80");
-    if (0==P)
-        return puts ("Oops"), 0;
-
-    /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
-    a = b = pj_obs_null ();
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
-    a.coo.lpz.z   = 100;
-
-    /* Forward projection: 3D-Cartesian-to-Ellipsoidal */
-    b = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);
-
-    /* Check roundtrip precision for 10000 iterations each way */
-    dist = pj_roundtrip (P, PJ_FWD, 10000, a);
-    printf ("Roundtrip deviation, fwd (nm): %15.9f\n", dist*1e9*1e5);
-    dist = pj_roundtrip (P, PJ_INV, 10000, b);
-    printf ("Roundtrip deviation, inv (nm): %15.9f\n", dist*1e9);
-
-    /* Inverse projection: Ellipsoidal-to-3D-Cartesian */
-    b = pj_trans (P, PJ_INV, b);
-    printf ("INV:   %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);
-
-    /* Move p to another context */
-    pj_context_renew (P);
-    b = pj_trans (P, PJ_FWD, b);
-    printf ("CTX1:  %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);
-
-    /* Move it back to the default context */
-    pj_context_free (P);
-    b = pj_trans (P, PJ_INV, b);
-    printf ("CTX0:  %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);
-
-    pj_free (P);
-
-    /***************************************************************************
-
-                  P I P E L I N E   T E S T S
-
-    ***************************************************************************/
-
-    /* forward-reverse geo->utm->geo */
-    P = pj_create (
-            "+proj=pipeline +ellps=GRS80 +zone=32 +step "
-            "+proj=utm +step "
-            "+proj=utm +inv"
-        );
-    if (0==P)
-        return puts ("Oops"), 0;
-    /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
-    a = b = pj_obs_null ();
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
-    a.coo.lpz.z   = 100;
-    printf ("PRE:   %15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
-
-    /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi));
-
-    /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
-    printf ("INV:   %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
-
-    pj_free (P);
-
-
-    /* And now the back-to-back situation utm->geo->utm */
-    P = pj_create (
-        "+proj=pipeline +ellps=GRS80 +zone=32 +step "
-        "+proj=utm +inv +step "
-        "+proj=utm");
-    if (0==P)
-        return puts ("Oops"), 0;
-
-    /* zero initialize everything, then set (easting, northing) to utm(12, 55) */
-    a = b = pj_obs_null ();
-    a.coo.xy = cph_utm32;
-    printf ("PRE:   %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);
-
-    /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);
-
-    /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
-    printf ("INV:   %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);
-
-    pj_free (P);
-
-
-
-    /* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */
-    P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm");
-    if (0==P)
-        return puts ("Oops"), 0;
-
-    /* zero initialize everything, then set (easting, northing) to utm(12, 55) */
-    a = b = pj_obs_null ();
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
-    printf ("PRE:   %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
-    printf ("EXP:   %15.9f %15.9f\n", cph_utm32.x, cph_utm32.y);
-
-    /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    printf ("FWD:   %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);
-
-    /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
-    printf ("INV:   %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
-
-    /* Geodesic distance between two points with angular 2D coordinates */
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(60);
-    b.coo.lp.lam = TORAD(12);
-    b.coo.lp.phi = TORAD(61);
-    dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
-    printf ("1 deg at 60N:   %15.9f\n", dist);
-
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(0.);
-    b.coo.lp.lam = TORAD(12);
-    b.coo.lp.phi = TORAD(1.);
-    dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
-    printf ("1 deg at Equator:   %15.9f\n", dist);
-
-    pj_free (P);
-    pj_pipeline_selftest ();
-    return 0;
-}
diff --git a/src/Makefile.am b/src/Makefile.am
index 25a7eaf4..f2615412 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -32,7 +32,7 @@ lib_LTLIBRARIES = libproj.la
 libproj_la_LDFLAGS = -no-undefined -version-info 12:0:0
 
 libproj_la_SOURCES = \
-	pj_list.h \
+	pj_list.h proj_internal.h\
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
 	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
@@ -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_hgridshift.c PJ_unitconvert.c
+	PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c pj_internal.c
 
 install-exec-local:
 	rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index f3379a1b..26acac2d 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -41,7 +41,7 @@
  *****************************************************************************/
 
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 #include <assert.h>
 #include <stddef.h>
@@ -232,117 +232,196 @@ int pj_cart_selftest (void) {return 0;}
 #else
 /* Testing quite a bit of the pj_obs_api as a side effect (inspired by pj_obs_api_test.c) */
 int pj_cart_selftest (void) {
+    PJ_CONTEXT *ctx;
     PJ *P;
-    PJ_OBS a, b;
+    PJ_OBS a, b, obs[2];
     int err;
-    double dist;
+    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;
 
     /* An utm projection on the GRS80 ellipsoid */
-    P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
+    P = proj_create (0, arg);
     if (0==P)
         return 1;
 
+    /* note arg is handcrafted to go undisturbed through get def reconstruction */
+    arg_def = proj_definition_retrieve (P);
+    if (0!=strcmp(arg, arg_def))
+        return 44;
+    proj_release (arg_def);
+
     /* Clean up */
-    pj_free (P);
+    proj_destroy (P);
 
     /* Same projection, now using argc/argv style initialization */
-    P = pj_create_argv (3, args);
+    P = proj_create_argv (0, 3, args);
     if (0==P)
         return 2;
 
     /* zero initialize everything, then set (longitude, latitude) to (12, 55) */
-    a = pj_obs_null;
+    a = proj_obs_null;
     /* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(55);
+    a.coo.lp.lam = PJ_TORAD(12);
+    a.coo.lp.phi = PJ_TORAD(55);
 
     /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
+    b = proj_trans_obs (P, PJ_FWD, a);
 
     /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
+    a = proj_trans_obs (P, PJ_INV, b);
 
     /* Null projection */
-    a = pj_trans (P, PJ_IDENT, a);
+    a = proj_trans_obs (P, PJ_IDENT, a);
 
     /* Forward again, to get two linear items for comparison */
-    a = pj_trans (P, PJ_FWD, a);
+    a = proj_trans_obs (P, PJ_FWD, a);
 
-    dist = pj_xy_dist (a.coo.xy, b.coo.xy);
+    dist = proj_xy_dist (a.coo.xy, b.coo.xy);
     if (dist > 2e-9)
         return 3;
 
+    /* Clear any previous error */
+    proj_errno_set (P, 0);
+
     /* Invalid projection */
-    a = pj_trans (P, 42, a);
+    a = proj_trans_obs (P, 42, a);
     if (a.coo.lpz.lam!=HUGE_VAL)
         return 4;
-    err = pj_err_level (P, PJ_ERR_TELL);
+    err = proj_errno (P);
     if (0==err)
         return 5;
 
-    /* Clear error */
-    pj_err_level (P, 0);
+    /* Clear error again */
+    proj_errno_set (P, 0);
 
     /* Clean up */
-    pj_free (P);
+    proj_destroy (P);
 
     /* Now do some 3D transformations */
-    P = pj_create ("+proj=cart +ellps=GRS80");
+    P = proj_create (0, "+proj=cart +ellps=GRS80");
     if (0==P)
         return 6;
 
     /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
-    a = b = pj_obs_null;
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
+    a = b = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(12);
+    a.coo.lpz.phi = PJ_TORAD(55);
     a.coo.lpz.z   = 100;
 
     /* Forward projection: 3D-Cartesian-to-Ellipsoidal */
-    b = pj_trans (P, PJ_FWD, a);
+    b = proj_trans_obs (P, PJ_FWD, a);
 
     /* Check roundtrip precision for 10000 iterations each way */
-    dist = pj_roundtrip (P, PJ_FWD, 10000, a);
-    dist = pj_roundtrip (P, PJ_INV, 10000, b);
+    dist = proj_roundtrip (P, PJ_FWD, 10000, a);
+    dist = proj_roundtrip (P, PJ_INV, 10000, b);
     if (dist > 2e-9)
         return 7;
 
 
     /* Test at the North Pole */
-    a = b = pj_obs_null;
-    a.coo.lpz.lam = TORAD(0);
-    a.coo.lpz.phi = TORAD(90);
+    a = b = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(0);
+    a.coo.lpz.phi = PJ_TORAD(90);
     a.coo.lpz.z   = 100;
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 1e-12)
         return 8;
 
     /* Test at the South Pole */
-    a = b = pj_obs_null;
-    a.coo.lpz.lam = TORAD(0);
-    a.coo.lpz.phi = TORAD(-90);
+    a = b = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(0);
+    a.coo.lpz.phi = PJ_TORAD(-90);
     a.coo.lpz.z   = 100;
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 1e-12)
         return 9;
 
     /* Inverse projection: 3D-Cartesian-to-Ellipsoidal */
-    b = pj_trans (P, PJ_INV, b);
+    b = proj_trans_obs (P, PJ_INV, b);
 
     /* Move p to another context */
-    pj_context_renew (P);
-    b = pj_trans (P, PJ_FWD, b);
+    ctx = proj_context_create ();
+    if (ctx==pj_get_default_ctx())
+        return 10;
+    proj_context_set (P, ctx);
+    if (ctx != P->ctx)
+        return 11;
+    b = proj_trans_obs (P, PJ_FWD, b);
 
     /* Move it back to the default context */
-    pj_context_free (P);
-    b = pj_trans (P, PJ_INV, b);
+    proj_context_set (P, 0);
+    if (pj_get_default_ctx() != P->ctx)
+        return 12;
+    proj_context_destroy (ctx);
+
+    /* We go on with the work - now back on the default context */
+    b = proj_trans_obs (P, PJ_INV, b);
+    proj_destroy (P);
+
+
+    /* Testing the proj_transform nightmare */
+
+    /* An utm projection on the GRS80 ellipsoid */
+    P = proj_create (0, "+proj=utm +zone=32 +ellps=GRS80");
+    if (0==P)
+        return 13;
+
+    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);
+    sz = sizeof (PJ_OBS);
+
+    /* Forward projection */
+    a = proj_trans_obs (P, PJ_FWD, obs[0]);
+    b = proj_trans_obs (P, PJ_FWD, obs[1]);
+
+    n = proj_transform (
+        P, PJ_FWD, 
+        &(obs[0].coo.lpz.lam), sz, 2,
+        &(obs[0].coo.lpz.phi), sz, 2,
+        &(obs[0].coo.lpz.z),   sz, 2, 
+        0,                     sz, 0
+    );
+    if (2!=n)
+        return 14;
+    if (a.coo.lpz.lam != obs[0].coo.lpz.lam)  return 15;
+    if (a.coo.lpz.phi != obs[0].coo.lpz.phi)  return 16;
+    if (a.coo.lpz.z   != obs[0].coo.lpz.z)    return 17;
+    if (b.coo.lpz.lam != obs[1].coo.lpz.lam)  return 18;
+    if (b.coo.lpz.phi != obs[1].coo.lpz.phi)  return 19;
+    if (b.coo.lpz.z   != obs[1].coo.lpz.z)    return 20;
+
+    /* now test the case of constant z */
+    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);
+    h = 27;
+    t = 33;
+    n = proj_transform (
+        P, PJ_FWD, 
+        &(obs[0].coo.lpz.lam), sz, 2,
+        &(obs[0].coo.lpz.phi), sz, 2,
+        &h,                     0, 1, 
+        &t,                     0, 1
+    );
+    if (2!=n)
+        return 21;
+    if (a.coo.lpz.lam != obs[0].coo.lpz.lam)  return 22;
+    if (a.coo.lpz.phi != obs[0].coo.lpz.phi)  return 23;
+    if (45            != obs[0].coo.lpz.z)    return 24;
+    if (b.coo.lpz.lam != obs[1].coo.lpz.lam)  return 25;
+    if (b.coo.lpz.phi != obs[1].coo.lpz.phi)  return 26;
+    if (50            != obs[1].coo.lpz.z)    return 27; /* NOTE: unchanged */
+    if (50==h) return 28;
+
+    /* Clean up */
+    proj_destroy (P);
 
-    pj_free (P);
     return 0;
 }
 #endif
diff --git a/src/PJ_chamb.c b/src/PJ_chamb.c
index eb3ed479..e9e9848b 100644
--- a/src/PJ_chamb.c
+++ b/src/PJ_chamb.c
@@ -156,7 +156,7 @@ int pj_chamb_selftest (void) {
     double tolerance_lp = 1e-10;
     double tolerance_xy = 1e-7;
 
-    char s_args[] = {"+proj=chamb   +a=6400000    +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=chamb   +R=6400000    +lat_1=0.5 +lat_2=2"};
 
     LP fwd_in[] = {
         { 2, 1},
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index d5d5a683..e291c189 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -43,7 +43,7 @@ Last update: 2017-05-15
 ***********************************************************************/
 
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 #include <geocent.h>
 #include <assert.h>
@@ -154,16 +154,16 @@ static void update_parameters(PJ *P) {
     Q->theta = Q->theta_0 + Q->dtheta * dt;
 
     /* debugging output */
-    if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_TRACE) {
-        pj_log_trace(P, "Transformation parameters for observation epoch %g:", Q->t_obs);
-        pj_log_trace(P, "x: %g", Q->xyz.x);
-        pj_log_trace(P, "y: %g", Q->xyz.y);
-        pj_log_trace(P, "z: %g", Q->xyz.z);
-        pj_log_trace(P, "s: %g", Q->scale*1e-6);
-        pj_log_trace(P, "rx: %g", Q->opk.o);
-        pj_log_trace(P, "ry: %g", Q->opk.p);
-        pj_log_trace(P, "rz: %g", Q->opk.k);
-        pj_log_trace(P, "theta: %g", Q->theta);
+    if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) {
+        proj_log_trace(P, "Transformation parameters for observation epoch %g:", Q->t_obs);
+        proj_log_trace(P, "x: %g", Q->xyz.x);
+        proj_log_trace(P, "y: %g", Q->xyz.y);
+        proj_log_trace(P, "z: %g", Q->xyz.z);
+        proj_log_trace(P, "s: %g", Q->scale*1e-6);
+        proj_log_trace(P, "rx: %g", Q->opk.o);
+        proj_log_trace(P, "ry: %g", Q->opk.p);
+        proj_log_trace(P, "rz: %g", Q->opk.k);
+        proj_log_trace(P, "theta: %g", Q->theta);
     }
     return;
 }
@@ -312,11 +312,11 @@ static void build_rot_matrix(PJ *P) {
     }
 
     /* some debugging output */
-    if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_TRACE) {
-        pj_log_trace(P, "Rotation Matrix:");
-        pj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R00, R01, R02);
-        pj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R10, R11, R12);
-        pj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R20, R21, R22);
+    if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_TRACE) {
+        proj_log_trace(P, "Rotation Matrix:");
+        proj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R00, R01, R02);
+        proj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R10, R11, R12);
+        proj_log_trace(P, "  | % 6.6g  % 6.6g  % 6.6g |", R20, R21, R22);
     }
 
     return;
@@ -579,16 +579,16 @@ PJ *PROJECTION(helmert) {
     Q->theta  =  Q->theta_0;
 
     /* Let's help with debugging */
-    if (pj_err_level(P, PJ_ERR_TELL) >= PJ_LOG_DEBUG) {
-        pj_log_debug(P, "Helmert parameters:");
-        pj_log_debug(P, "x=  % 3.5f  y=  % 3.5f  z=  % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
-        pj_log_debug(P, "rx= % 3.5f  ry= % 3.5f  rz= % 3.5f",
+    if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
+        proj_log_debug(P, "Helmert parameters:");
+        proj_log_debug(P, "x=  % 3.5f  y=  % 3.5f  z=  % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
+        proj_log_debug(P, "rx= % 3.5f  ry= % 3.5f  rz= % 3.5f",
                 Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
-        pj_log_debug(P, "s=% 3.5f  approximate=% d  transpose=% d",
+        proj_log_debug(P, "s=% 3.5f  approximate=% d  transpose=% d",
                 Q->scale, Q->approximate, Q->transpose);
-        pj_log_debug(P, "dx= % 3.5f  dy= % 3.5f  dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
-        pj_log_debug(P, "drx=% 3.5f  dry=% 3.5f  drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
-        pj_log_debug(P, "ds=% 3.5f  epoch=% 5.5f  tobs=% 5.5f", Q->dscale, Q->epoch, Q->t_obs);
+        proj_log_debug(P, "dx= % 3.5f  dy= % 3.5f  dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
+        proj_log_debug(P, "drx=% 3.5f  dry=% 3.5f  drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
+        proj_log_debug(P, "ds=% 3.5f  epoch=% 5.5f  tobs=% 5.5f", Q->dscale, Q->epoch, Q->t_obs);
     }
 
     if ((Q->opk.o==0) && (Q->opk.p==0) && (Q->opk.k==0) && (Q->scale==0)) {
@@ -619,23 +619,26 @@ static int test (char *args, PJ_TRIPLET in, PJ_TRIPLET expect, double tol) {
         return 5;
 
     out.xyz = pj_fwd3d (in.lpz, P);
-    if (pj_xyz_dist (out.xyz, expect.xyz) > tol) {
-        pj_log_error(P, "Tolerance of forward calculation not met");
-        pj_log_error(P, "  Expect: %10.10f, %10.10f, %10.10f", expect.xyz.x, expect.xyz.y, expect.xyz.z);
-        pj_log_error(P, "  Out:    %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z);
-        pj_log_level(NULL, 0);
+    if (proj_xyz_dist (out.xyz, expect.xyz) > tol) {
+        proj_log_error(P, "Tolerance of forward calculation not met");
+        proj_log_error(P, "  Expect: %10.10f, %10.10f, %10.10f", expect.xyz.x, expect.xyz.y, expect.xyz.z);
+        proj_log_error(P, "  Out:    %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z);
+        proj_log_level(NULL, 0);
+        proj_destroy (P);
         return 1;
     }
 
     out.lpz = pj_inv3d (out.xyz, P);
-    if (pj_xyz_dist (out.xyz, in.xyz) > tol) {
-        pj_log_error(P, "Tolerance of inverse calculation not met");
-        pj_log_error(P, "  In:  %10.10f, %10.10f, %10.10f", in.xyz.x, in.xyz.y, in.xyz.z);
-        pj_log_error(P, "  Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z);
-        pj_log_level(NULL, 0);
+    if (proj_xyz_dist (out.xyz, in.xyz) > tol) {
+        proj_log_error(P, "Tolerance of inverse calculation not met");
+        proj_log_error(P, "  In:  %10.10f, %10.10f, %10.10f", in.xyz.x, in.xyz.y, in.xyz.z);
+        proj_log_error(P, "  Out: %10.10f, %10.10f, %10.10f", out.xyz.x, out.xyz.y, out.xyz.z);
+        proj_log_level(NULL, 0);
+        proj_destroy (P);
         return 2;
     }
 
+    proj_destroy (P);
     return 0;
 }
 
@@ -675,7 +678,7 @@ int pj_helmert_selftest (void) {
     PJ_TRIPLET in3 = {{3370658.378, 711877.314, 5349787.086}}; /* ITRF2000@2017.0 */
     PJ_TRIPLET expect3 = {{3370658.18890, 711877.42370, 5349787.12430}}; /* ITRF93@2017.0 */
     char args3[] = {
-        " +proj=helmert"
+        " +proj=helmert +ellps=GRS80"
         " +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195"
         " +rx=-0.00039 +ry=0.00080 +rz=-0.00114"
         " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001"
@@ -690,8 +693,9 @@ int pj_helmert_selftest (void) {
     PJ_OBS in4 = {{{3370658.378, 711877.314, 5349787.086, 2017.0}}, {{ 0, 0, 0}}, 0, 0};
     PJ_OBS out;
 
-    PJ *helmert = pj_create(
-        " +proj=helmert"
+    PJ *helmert = proj_create(
+        0, 
+        " +proj=helmert +ellps=GRS80"
         " +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195"
         " +rx=-0.00039 +ry=0.00080 +rz=-0.00114"
         " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001"
@@ -709,7 +713,7 @@ int pj_helmert_selftest (void) {
        to computed coordinates, hence the test tolerance is set accordingly. */
     PJ_TRIPLET in5 = {{2546506.957, 542256.609, 0}};
     PJ_TRIPLET expect5 = {{766563.675, 165282.277, 0}};
-    char args5[] = " +proj=helmert +x=-9597.3572 +y=.6112 +s=0.304794780637 +theta=-1.244048";
+    char args5[] = " +proj=helmert +ellps=GRS80 +x=-9597.3572 +y=.6112 +s=0.304794780637 +theta=-1.244048";
 
     /* Run tests 1-3 & 5 */
     ret = test (args1, in1, expect1, 1e-4);   if (ret)  return ret;
@@ -718,24 +722,24 @@ int pj_helmert_selftest (void) {
     ret = test (args5, in5, expect5, 0.001);  if (ret) return ret + 40;
 
     /* Run test 4 */
-    out = pj_trans(helmert, PJ_FWD, in4);
-    if (pj_xyz_dist (out.coo.xyz, expect4a) > 1e-4) {
-        pj_log_error(helmert, "Tolerance of test 4a not met!");
-        pj_log_error(helmert, "  In:  %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z);
-        pj_log_error(helmert, "  Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
+    out = proj_trans_obs (helmert, PJ_FWD, in4);
+    if (proj_xyz_dist (out.coo.xyz, expect4a) > 1e-4) {
+        proj_log_error(helmert, "Tolerance of test 4a not met!");
+        proj_log_error(helmert, "  In:  %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z);
+        proj_log_error(helmert, "  Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
         return 31;
     }
 
     in4.coo.xyzt.t = 2018.0;
-    out = pj_trans(helmert, PJ_FWD, in4);
-    if (pj_xyz_dist (out.coo.xyz, expect4b) > 1e-4) {
-        pj_log_error(helmert, "Tolerance of test 4b not met!");
-        pj_log_error(helmert, "  In:  %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z);
-        pj_log_error(helmert, "  Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
+    out = proj_trans_obs (helmert, PJ_FWD, in4);
+    if (proj_xyz_dist (out.coo.xyz, expect4b) > 1e-4) {
+        proj_log_error(helmert, "Tolerance of test 4b not met!");
+        proj_log_error(helmert, "  In:  %10.10f, %10.10f, %10.10f", in4.coo.xyz.x, in4.coo.xyz.y, in4.coo.xyz.z);
+        proj_log_error(helmert, "  Out: %10.10f, %10.10f, %10.10f", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
         return 32;
     }
-    pj_free(helmert);
 
+    proj_destroy(helmert);
     return 0;
 }
 
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 1b9d2727..41c2b629 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -1,5 +1,5 @@
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
@@ -70,7 +70,7 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
 PJ *PROJECTION(hgridshift) {
 
     if (!pj_param(P->ctx, P->params, "tgrids").i) {
-        pj_log_error(P, "hgridshift: +grids parameter missing.");
+        proj_log_error(P, "hgridshift: +grids parameter missing.");
         return freeup_msg(P, -1);
     }
 
@@ -80,7 +80,7 @@ PJ *PROJECTION(hgridshift) {
 
     /* Was gridlist compiled properly? */
     if ( pj_ctx_get_errno(P->ctx) ) {
-        pj_log_error(P, "hgridshift: could not find required grid(s).");
+        proj_log_error(P, "hgridshift: could not find required grid(s).");
         return freeup_msg(P, -38);
     }
 
@@ -108,38 +108,36 @@ int pj_hgridshift_selftest (void) {
     double dist;
 
     /* fail on purpose: +grids parameter is mandatory*/
-    P = pj_create("+proj=hgridshift");
+    P = proj_create(0, "+proj=hgridshift");
     if (0!=P)
         return 99;
 
     /* fail on purpose: open non-existing grid */
-    P = pj_create("+proj=hgridshift +grids=nonexistinggrid.gsb");
+    P = proj_create(0, "+proj=hgridshift +grids=nonexistinggrid.gsb");
     if (0!=P)
         return 999;
 
 
     /* Failure most likely means the grid is missing */
-    P = pj_create ("+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80");
+    P = proj_create (0, "+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);
+    a = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(173);
+    a.coo.lpz.phi = PJ_TORAD(-45);
 
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 0.00000001)
         return 1;
 
-    expect.coo.lpz.lam = TORAD(172.999892181021551);
-    expect.coo.lpz.phi = TORAD(-45.001620431954613);
-    b = pj_trans(P, PJ_FWD, a);
-    if (pj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4)
+    expect.coo.lpz.lam = PJ_TORAD(172.999892181021551);
+    expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613);
+    b = proj_trans_obs(P, PJ_FWD, a);
+    if (proj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4)
         return 2;
 
-
-    pj_free(P);
-
+    proj_destroy(P);
     return 0;
 }
 #endif
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
index 1b055a4e..56980f72 100644
--- a/src/PJ_horner.c
+++ b/src/PJ_horner.c
@@ -1,5 +1,5 @@
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 #include <assert.h>
 #include <stddef.h>
@@ -94,7 +94,7 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
 
 struct horner;
 typedef struct horner HORNER;
-static UV      horner (const HORNER *transformation, enum pj_direction, UV position);
+static UV      horner (const HORNER *transformation, enum proj_direction, UV position);
 static HORNER *horner_alloc (size_t order, int complex_polynomia);
 static void    horner_free (HORNER *h);
 
@@ -176,7 +176,7 @@ static HORNER *horner_alloc (size_t order, int complex_polynomia) {
 
 
 /**********************************************************************/
-static UV horner (const HORNER *transformation, enum pj_direction direction, UV position) {
+static UV horner (const HORNER *transformation, enum proj_direction direction, UV position) {
 /***********************************************************************
 
 A reimplementation of the classic Engsager/Poder 2D Horner polynomial
@@ -301,7 +301,7 @@ static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
 
 
 /**********************************************************************/
-static UV complex_horner (const HORNER *transformation, enum pj_direction direction, UV position) {
+static UV complex_horner (const HORNER *transformation, enum proj_direction direction, UV position) {
 /***********************************************************************
 
 A reimplementation of a classic Engsager/Poder Horner complex
@@ -401,7 +401,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
 
     buf = pj_calloc (strlen (param) + 2, sizeof(char));
     if (0==buf) {
-        pj_log_error (P, "Horner: Out of core");
+        proj_log_error (P, "Horner: Out of core");
         return 0;
     }
 
@@ -415,7 +415,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
     for (i = 0; i < ncoefs; i++) {
         if (i > 0) {
             if ( next == 0 || ','!=*next) {
-                pj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs);
+                proj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs);
                 return 0;
             }
             init = ++next;
@@ -443,7 +443,7 @@ PJ *PROJECTION(horner) {
     if (pj_param (P->ctx, P->params, "tdeg").i) /* degree specified? */
 		degree = pj_param(P->ctx, P->params, "ideg").i;
     else {
-        pj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)");
+        proj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)");
         return horner_freeup (P);
     }
 
@@ -525,48 +525,49 @@ int pj_horner_selftest (void) {
     double dist;
 
     /* Real polynonia relating the technical coordinate system TC32 to "System 45 Bornholm" */
-    P = pj_create (tc32_utm32);
+    P = proj_create (0, tc32_utm32);
     if (0==P)
         return 10;
 
-    a = b = pj_obs_null;
+    a = b = proj_obs_null;
     a.coo.uv.v = 6125305.4245;
     a.coo.uv.u =  878354.8539;
 
     /* Check roundtrip precision for 1 iteration each way, starting in forward direction */
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 0.01)
         return 1;
 
     /* The complex polynomial transformation between the "System Storebaelt" and utm32/ed50 */
-    P = pj_create (sb_utm32);
+    P = proj_create (0, sb_utm32);
     if (0==P)
         return 11;
 
     /* Test value: utm32_ed50(620000, 6130000) = sb_ed50(495136.8544, 6130821.2945) */
-    a = b = c = pj_obs_null;
+    a = b = c = proj_obs_null;
     a.coo.uv.v = 6130821.2945;
     a.coo.uv.u =  495136.8544;
     c.coo.uv.v = 6130000.0000;
     c.coo.uv.u =  620000.0000;
 
     /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    dist = pj_xy_dist (b.coo.xy, c.coo.xy);
+    b = proj_trans_obs (P, PJ_FWD, a);
+    dist = proj_xy_dist (b.coo.xy, c.coo.xy);
     if (dist > 0.001)
         return 2;
 
     /* Inverse projection */
-    b = pj_trans (P, PJ_INV, c);
-    dist = pj_xy_dist (b.coo.xy, a.coo.xy);
+    b = proj_trans_obs (P, PJ_INV, c);
+    dist = proj_xy_dist (b.coo.xy, a.coo.xy);
     if (dist > 0.001)
         return 3;
 
     /* Check roundtrip precision for 1 iteration each way */
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 0.01)
         return 4;
 
+    proj_destroy(P);
     return 0;
 }
 #endif
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index 387c0ca8..6e42c0cb 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -96,7 +96,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 ********************************************************************************/
 
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 
 #include <assert.h>
@@ -183,7 +183,7 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
     last_step  = P->opaque->steps + 1;
 
     for (i = first_step;  i != last_step;  i++) {
-        pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s",
+        proj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s",
             i-first_step, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
             P->opaque->pipeline[i]->descr
         );
@@ -191,13 +191,13 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
         if (P->opaque->omit_forward[i])
             continue;
         if (P->opaque->reverse_step[i])
-            point = pj_trans (P->opaque->pipeline[i], PJ_INV, point);
+            point = proj_trans_obs (P->opaque->pipeline[i], PJ_INV, point);
         else
-            point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point);
+            point = proj_trans_obs (P->opaque->pipeline[i], PJ_FWD, point);
         if (P->opaque->depth < PIPELINE_STACK_SIZE)
             P->opaque->stack[P->opaque->depth++] = point;
     }
-    pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
+    proj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
 
     P->opaque->depth = 0;    /* Clear the stack */
     return point;
@@ -210,20 +210,20 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
     first_step = P->opaque->steps;
     last_step  =  0;
     for (i = first_step;  i != last_step;  i--) {
-        pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f",
+        proj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f",
             i - 1, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
             P->opaque->pipeline[i]->a, P->opaque->pipeline[i]->rf
         );
         if (P->opaque->omit_inverse[i])
             continue;
         if (P->opaque->reverse_step[i])
-            point = pj_trans (P->opaque->pipeline[i], PJ_FWD, point);
+            point = proj_trans_obs (P->opaque->pipeline[i], PJ_FWD, point);
         else
-            point = pj_trans (P->opaque->pipeline[i], PJ_INV, point);
+            point = proj_trans_obs (P->opaque->pipeline[i], PJ_INV, point);
         if (P->opaque->depth < PIPELINE_STACK_SIZE)
             P->opaque->stack[P->opaque->depth++] = point;
     }
-    pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
+    proj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
 
     P->opaque->depth = 0;    /* Clear the stack */
     return point;
@@ -232,7 +232,7 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
 
 /* Delegate the work to pipeline_forward_obs() */
 static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) {
-    PJ_OBS point = pj_obs_null;
+    PJ_OBS point = proj_obs_null;
     point.coo.lpz = lpz;
     point = pipeline_forward_obs (point, P);
     return point.coo.xyz;
@@ -240,21 +240,21 @@ static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) {
 
 /* Delegate the work to pipeline_reverse_obs() */
 static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P) {
-    PJ_OBS point = pj_obs_null;
+    PJ_OBS point = proj_obs_null;
     point.coo.xyz = xyz;
     point = pipeline_reverse_obs (point, P);
     return point.coo.lpz;
 }
 
 static XY pipeline_forward (LP lp, PJ *P) {
-    PJ_OBS point = pj_obs_null;
+    PJ_OBS point = proj_obs_null;
     point.coo.lp = lp;
     point = pipeline_forward_obs (point, P);
     return point.coo.xy;
 }
 
 static LP pipeline_reverse (XY xy, PJ *P) {
-    PJ_OBS point = pj_obs_null;
+    PJ_OBS point = proj_obs_null;
     point.coo.xy = xy;
     point = pipeline_reverse_obs (point, P);
     return point.coo.lp;
@@ -275,7 +275,8 @@ static void *pipeline_freeup (PJ *P, int errlev) {         /* Destructor */
     if (0==P)
         return 0;
 
-    pj_err_level (P, errlev);
+    if (errlev)
+        proj_errno_set (P, errlev);
 
     if (0==P->opaque)
         return pj_dealloc (P);
@@ -382,7 +383,7 @@ PJ *PROJECTION(pipeline) {
     for (i = 0;  i < argc;  i++) {
         if (0==strcmp ("step", argv[i])) {
             if (-1==i_pipeline) {
-                pj_log_error (P, "Pipeline: +step before +proj=pipeline");
+                proj_log_error (P, "Pipeline: +step before +proj=pipeline");
                 return pipeline_freeup (P, -50);
             }
             if (0==nsteps)
@@ -393,7 +394,7 @@ PJ *PROJECTION(pipeline) {
 
         if (0==strcmp ("proj=pipeline", argv[i])) {
             if (-1 != i_pipeline) {
-                pj_log_error (P, "Pipeline: Nesting invalid");
+                proj_log_error (P, "Pipeline: Nesting invalid");
                 return pipeline_freeup (P, -50); /* ERROR: nested pipelines */
             }
             i_pipeline = i;
@@ -419,7 +420,7 @@ PJ *PROJECTION(pipeline) {
         PJ     *next_step = 0;
 
         /* Build a set of setup args for the current step */
-        pj_log_trace (P, "Pipeline: Building arg list for step no. %d", i);
+        proj_log_trace (P, "Pipeline: Building arg list for step no. %d", i);
 
         /* First add the step specific args */
         for (j = i_current_step + 1;  0 != strcmp ("step", argv[j]); j++)
@@ -445,21 +446,21 @@ PJ *PROJECTION(pipeline) {
                 P->opaque->reverse_step[i+1] = 1;
         }
 
-        pj_log_trace (P, "Pipeline: init - %s, %d", current_argv[0], current_argc);
+        proj_log_trace (P, "Pipeline: init - %s, %d", current_argv[0], current_argc);
         for (j = 1;  j < current_argc; j++)
-            pj_log_trace (P, "    %s", current_argv[j]);
+            proj_log_trace (P, "    %s", current_argv[j]);
 
         next_step = pj_init (current_argc, current_argv);
-        pj_log_trace (P, "Pipeline: Step %d at %p", i, next_step);
+        proj_log_trace (P, "Pipeline: Step %d at %p", i, next_step);
         if (0==next_step) {
-            pj_log_error (P, "Pipeline: Bad step definition: %s", current_argv[0]);
+            proj_log_error (P, "Pipeline: Bad step definition: %s", current_argv[0]);
             return pipeline_freeup (P, -50); /* ERROR: bad pipeline def */
         }
         P->opaque->pipeline[i+1] = next_step;
-        pj_log_trace (P, "Pipeline:    step done");
+        proj_log_trace (P, "Pipeline:    step done");
     }
 
-    pj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps);
+    proj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps);
 
     /* Determine forward input (= reverse output) data type */
 
@@ -468,7 +469,7 @@ PJ *PROJECTION(pipeline) {
         if (0==P->opaque->omit_forward[i+1])
             break;
     if (i==nsteps) {
-        pj_log_error (P, "Pipeline: No forward steps");
+        proj_log_error (P, "Pipeline: No forward steps");
         return pipeline_freeup (P, -50);
     }
 
@@ -491,7 +492,7 @@ PJ *PROJECTION(pipeline) {
         if (0==P->opaque->omit_inverse[i+1])
             break;
     if (i==-1) {
-        pj_log (pj_get_ctx (P), PJ_LOG_ERROR, "Pipeline: No reverse steps");
+        proj_log_error (P, "Pipeline: No reverse steps");
         return pipeline_freeup (P, -50);
     }
 
@@ -506,7 +507,7 @@ PJ *PROJECTION(pipeline) {
         else
             P->right = PJ_IO_UNITS_METERS;
     }
-    pj_log_trace (P, "Pipeline: Units - left: [%s], right: [%s]\n",
+    proj_log_trace (P, "Pipeline: Units - left: [%s], right: [%s]\n",
         P->left ==PJ_IO_UNITS_RADIANS? "angular": "linear",
         P->right==PJ_IO_UNITS_RADIANS? "angular": "linear");
 
@@ -525,94 +526,96 @@ int pj_pipeline_selftest (void) {
     double dist;
 
     /* forward-reverse geo->utm->geo */
-    P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +step +proj=utm +ellps=GRS80 +inv");
+    P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +step +proj=utm +ellps=GRS80 +inv");
     if (0==P)
         return 1000;
     /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 0) */
-    a = b = pj_obs_null;
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
+    a = b = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(12);
+    a.coo.lpz.phi = PJ_TORAD(55);
     a.coo.lpz.z   = 0;
 
     /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
+    b = proj_trans_obs (P, PJ_FWD, a);
+    if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
         return 1001;
 
     /* Inverse projection (still same result: pipeline is symmetrical) */
-    a = pj_trans (P, PJ_INV, b);
-    if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
+    a = proj_trans_obs (P, PJ_INV, b);
+    if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
         return 1002;
 
-    pj_free (P);
+    proj_destroy (P);
 
     /* And now the back-to-back situation utm->geo->utm */
-    P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +inv +step +proj=utm +ellps=GRS80");
+    P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 +inv +step +proj=utm +ellps=GRS80");
     if (0==P)
         return 2000;
 
     /* zero initialize everything, then set (easting, northing) to utm(12, 55) */
-    a = b = pj_obs_null;
+    a = b = proj_obs_null;
     a.coo.xy = cph_utm32;
 
     /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    if (pj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4)
+    b = proj_trans_obs (P, PJ_FWD, a);
+    if (proj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4)
         return 2001;
 
     /* Inverse projection */
-    a = pj_trans (P, PJ_INV, b);
-    if (pj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4)
+    a = proj_trans_obs (P, PJ_INV, b);
+    if (proj_xy_dist (a.coo.xy, b.coo.xy) > 1e-4)
         return 2001;
-    if (pj_xyz_dist (a.coo.xyz, b.coo.xyz) > 1e-4)
+    if (proj_xyz_dist (a.coo.xyz, b.coo.xyz) > 1e-4)
         return 2002;
 
-    pj_free (P);
+    proj_destroy (P);
 
 
     /* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */
-    P = pj_create ("+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 ");
+    P = proj_create (0, "+proj=pipeline +zone=32 +step +proj=utm +ellps=GRS80 ");
     if (0==P)
         return 3000;
 
 
-    a = b = pj_obs_null;
-    a.coo.lpz.lam = TORAD(12);
-    a.coo.lpz.phi = TORAD(55);
+    a = b = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(12);
+    a.coo.lpz.phi = PJ_TORAD(55);
 
     /* Forward projection */
-    b = pj_trans (P, PJ_FWD, a);
-    if (pj_xy_dist (cph_utm32, b.coo.xy) > 1e-4)
+    b = proj_trans_obs (P, PJ_FWD, a);
+    if (proj_xy_dist (cph_utm32, b.coo.xy) > 1e-4)
         return 3001;
 
     /* Inverse projection */
-    b = pj_trans (P, PJ_INV, b);
-    if (pj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
+    b = proj_trans_obs (P, PJ_INV, b);
+    if (proj_lp_dist (P, a.coo.lp, b.coo.lp) > 1e-4)
         return 3002;
 
 
     /* Since we use pj_lp_dist to determine success above, we should also test that it works */
 
     /* Geodesic distance between two points with angular 2D coordinates */
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(60);
-    b.coo.lp.lam = TORAD(12);
-    b.coo.lp.phi = TORAD(61);
-    dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
+    a.coo.lp.lam = PJ_TORAD(12);
+    a.coo.lp.phi = PJ_TORAD(60);
+    b.coo.lp.lam = PJ_TORAD(12);
+    b.coo.lp.phi = PJ_TORAD(61);
+    dist = proj_lp_dist (P, a.coo.lp, b.coo.lp);
     if (fabs (111420.727870234 - dist) > 1e-4)
         return 4001;
 
-    a.coo.lp.lam = TORAD(12);
-    a.coo.lp.phi = TORAD(0.);
-    b.coo.lp.lam = TORAD(12);
-    b.coo.lp.phi = TORAD(1.);
-    dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
+    a.coo.lp.lam = PJ_TORAD(12);
+    a.coo.lp.phi = PJ_TORAD(0.);
+    b.coo.lp.lam = PJ_TORAD(12);
+    b.coo.lp.phi = PJ_TORAD(1.);
+    dist = proj_lp_dist (P, a.coo.lp, b.coo.lp);
     if (fabs (110574.388554153 - dist) > 1e-4)
         return 4002;
 
+    proj_destroy (P);
 
     /* test a pipeline with several +init steps */
-    P = pj_create(
+    P = proj_create(
+        0,
         "+proj=pipeline "
         "+step +init=epsg:25832 +inv "
         "+step +init=epsg:25833  "
@@ -625,13 +628,13 @@ int pj_pipeline_selftest (void) {
     a.coo.xy.x =  700000.0;
     a.coo.xy.y = 6000000.0;
 
-    b = pj_trans(P, PJ_FWD, a);
-    dist = pj_xy_dist(a.coo.xy, b.coo.xy);
+    b = proj_trans_obs(P, PJ_FWD, a);
+    dist = proj_xy_dist(a.coo.xy, b.coo.xy);
     if (dist > 1e-7)
         return 5001;
 
 
-    pj_free (P);
+    proj_destroy (P);
     return 0;
 }
 #endif
diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index be9d2f70..27672c8d 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -65,8 +65,9 @@ Last update: 2017-05-16
 
 #define PJ_LIB__
 #include <time.h>
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
+#include <errno.h>
 
 PROJ_HEAD(unitconvert, "Unit conversion");
 
@@ -348,7 +349,7 @@ PJ *PROJECTION(unitconvert) {
         if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */
 
         Q->xy_in_id = i;
-        pj_log_debug(P, "xy_in unit: %s", pj_units[i].name);
+        proj_log_debug(P, "xy_in unit: %s", pj_units[i].name);
     }
 
     if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != NULL) {
@@ -357,7 +358,7 @@ PJ *PROJECTION(unitconvert) {
         if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */
 
         Q->xy_out_id = i;
-        pj_log_debug(P, "xy_out unit: %s", pj_units[i].name);
+        proj_log_debug(P, "xy_out unit: %s", pj_units[i].name);
     }
 
     if ((name = pj_param (P->ctx, P->params, "sz_in").s) != NULL) {
@@ -366,7 +367,7 @@ PJ *PROJECTION(unitconvert) {
         if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */
 
         Q->z_in_id = i;
-        pj_log_debug(P, "z_in unit: %s", pj_units[i].name);
+        proj_log_debug(P, "z_in unit: %s", pj_units[i].name);
     }
 
     if ((name = pj_param (P->ctx, P->params, "sz_out").s) != NULL) {
@@ -375,7 +376,7 @@ PJ *PROJECTION(unitconvert) {
         if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */
 
         Q->z_out_id = i;
-        pj_log_debug(P, "z_out unit: %s", pj_units[i].name);
+        proj_log_debug(P, "z_out unit: %s", pj_units[i].name);
     }
 
 
@@ -385,7 +386,7 @@ PJ *PROJECTION(unitconvert) {
         if (!s) return freeup_msg(P, -8); /* unknown unit conversion id */
 
         Q->t_in_id = i;
-        pj_log_debug(P, "t_in unit: %s", time_units[i].name);
+        proj_log_debug(P, "t_in unit: %s", time_units[i].name);
     }
 
     s = 0;
@@ -395,7 +396,7 @@ PJ *PROJECTION(unitconvert) {
             return freeup_msg(P, -8); /* unknown unit conversion id */
         }
         Q->t_out_id = i;
-        pj_log_debug(P, "t_out unit: %s", time_units[i].name);
+        proj_log_debug(P, "t_out unit: %s", time_units[i].name);
     }
 
     return P;
@@ -409,7 +410,7 @@ int pj_unitconvert_selftest (void) {return 0;}
 
 static int test_time(char* args, double tol, double t_in, double t_exp) {
     PJ_OBS in, out;
-    PJ *P = pj_create(args);
+    PJ *P = proj_create(0, args);
     int ret = 0;
 
     if (P == 0)
@@ -417,46 +418,46 @@ static int test_time(char* args, double tol, double t_in, double t_exp) {
 
     in.coo.xyzt.t = t_in;
 
-    out = pj_trans(P, PJ_FWD, in);
+    out = proj_trans_obs(P, PJ_FWD, in);
     if (fabs(out.coo.xyzt.t - t_exp) > tol) {
-        pj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_exp);
+        proj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_exp);
         ret = 1;
     }
-    out = pj_trans(P, PJ_INV, out);
+    out = proj_trans_obs(P, PJ_INV, out);
     if (fabs(out.coo.xyzt.t - t_in) > tol) {
-        pj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_in);
+        proj_log_error(P, "out: %10.10g, expect: %10.10g", out.coo.xyzt.t, t_in);
         ret = 2;
     }
     pj_free(P);
 
-    pj_log_level(NULL, 0);
+    proj_log_level(NULL, 0);
     return ret;
 }
 
 static int test_xyz(char* args, double tol, PJ_TRIPLET in, PJ_TRIPLET exp) {
     PJ_OBS out, obs_in;
-    PJ *P = pj_create(args);
+    PJ *P = proj_create(0, args);
     int ret = 0;
 
     if (P == 0)
         return 5;
 
     obs_in.coo.xyz = in.xyz;
-    out = pj_trans(P, PJ_FWD, obs_in);
-    if (pj_xyz_dist(out.coo.xyz, exp.xyz) > tol) {
+    out = proj_trans_obs(P, PJ_FWD, obs_in);
+    if (proj_xyz_dist(out.coo.xyz, exp.xyz) > tol) {
         printf("exp: %10.10g, %10.10g, %10.10g\n", exp.xyz.x, exp.xyz.y, exp.xyz.z);
         printf("out: %10.10g, %10.10g, %10.10g\n", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
         ret = 1;
     }
 
-    out = pj_trans(P, PJ_INV, out);
-    if (pj_xyz_dist(out.coo.xyz, in.xyz) > tol) {
+    out = proj_trans_obs(P, PJ_INV, out);
+    if (proj_xyz_dist(out.coo.xyz, in.xyz) > tol) {
         printf("exp: %g, %g, %g\n", in.xyz.x, in.xyz.y, in.xyz.z);
         printf("out: %g, %g, %g\n", out.coo.xyz.x, out.coo.xyz.y, out.coo.xyz.z);
         ret += 2;
     }
-    pj_free(P);
-    pj_log_level(NULL, 0);
+    proj_destroy(P);
+    proj_log_level(NULL, 0);
     return ret;
 }
 
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 72d97517..e790cd4d 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -1,5 +1,5 @@
 #define PJ_LIB__
-#include <proj.h>
+#include "proj_internal.h"
 #include <projects.h>
 
 PROJ_HEAD(vgridshift, "Vertical grid shift");
@@ -73,7 +73,7 @@ static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
 PJ *PROJECTION(vgridshift) {
 
    if (!pj_param(P->ctx, P->params, "tgrids").i) {
-        pj_log_error(P, "vgridshift: +grids parameter missing.");
+        proj_log_error(P, "vgridshift: +grids parameter missing.");
         return freeup_msg(P, -1);
     }
 
@@ -83,7 +83,7 @@ PJ *PROJECTION(vgridshift) {
 
     /* Was gridlist compiled properly? */
     if ( pj_ctx_get_errno(P->ctx) ) {
-        pj_log_error(P, "vgridshift: could not find required grid(s).");
+        proj_log_error(P, "vgridshift: could not find required grid(s).");
         pj_dalloc(P->gridlist);
         P->gridlist = NULL;
         return freeup_msg(P, -38);
@@ -113,32 +113,32 @@ int pj_vgridshift_selftest (void) {
     double dist;
 
     /* fail on purpose: +grids parameter it mandatory*/
-    P = pj_create("+proj=vgridshift");
+    P = proj_create(0, "+proj=vgridshift");
     if (0!=P)
         return 99;
 
     /* fail on purpose: open non-existing grid */
-    P = pj_create("+proj=vgridshift +grids=nonexistinggrid.gtx");
+    P = proj_create(0, "+proj=vgridshift +grids=nonexistinggrid.gtx");
     if (0!=P)
         return 999;
 
     /* Failure most likely means the grid is missing */
-    P = pj_create ("+proj=vgridshift +grids=egm96_15.gtx +ellps=GRS80");
+    P = proj_create (0, "+proj=vgridshift +grids=egm96_15.gtx +ellps=GRS80");
     if (0==P)
         return 10;
 
-    a = pj_obs_null;
-    a.coo.lpz.lam = TORAD(12.5);
-    a.coo.lpz.phi = TORAD(55.5);
+    a = proj_obs_null;
+    a.coo.lpz.lam = PJ_TORAD(12.5);
+    a.coo.lpz.phi = PJ_TORAD(55.5);
 
-    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 0.00000001)
         return 1;
 
     expect = a;
     expect.coo.lpz.z   = -36.021305084228515625;
-    b = pj_trans(P, PJ_FWD, a);
-    if (pj_xyz_dist(expect.coo.xyz, b.coo.xyz) > 1e-10)
+    b = proj_trans_obs(P, PJ_FWD, a);
+    if (proj_xyz_dist(expect.coo.xyz, b.coo.xyz) > 1e-10)
         return 2;
 
 
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 9ed1623c..39dab3a4 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -201,6 +201,8 @@ SET(SRC_LIBPROJ_CORE
         pj_msfn.c
         pj_mutex.c
         pj_obs_api.c
+        pj_internal.c
+        proj_internal.h
         pj_open_lib.c
         pj_param.c
         pj_phi2.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 8584d419..b37d3b7b 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -56,7 +56,8 @@ support = \
 	pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \
 	proj_mdist.obj pj_mutex.obj pj_initcache.obj \
 	pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \
-	pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj
+	pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj \
+	pj_internal.obj
 
 pipeline = \
 	pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \
diff --git a/src/pj_ctx.c b/src/pj_ctx.c
index 9ae3d4e0..89b2816f 100644
--- a/src/pj_ctx.c
+++ b/src/pj_ctx.c
@@ -60,8 +60,13 @@ void pj_set_ctx( projPJ pj, projCtx ctx )
 projCtx pj_get_default_ctx()
 
 {
+    /* If already initialized, don't bother locking */
+    if( default_context_initialized )
+        return &default_context;
+
     pj_acquire_lock();
 
+    /* Ask again, since it may have been initialized in another thread */
     if( !default_context_initialized )
     {
         default_context.last_errno = 0;
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 598c6dd2..02f18d1c 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -1,21 +1,21 @@
 /* general forward projection */
 #define PJ_LIB__
+#include <proj.h>
 #include <projects.h>
-#include <errno.h>
 # define EPS 1.0e-12
 	XY /* forward projection entry */
 pj_fwd(LP lp, PJ *P) {
 	XY xy;
     XY err;
     double t;
+    int last_errno;
 
     /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
     err.x = err.y = HUGE_VAL;
 
     if (0==P->fwd)
         return err;
-
-    P->ctx->last_errno = pj_errno = errno = 0;
+    last_errno = proj_errno_reset (P);
 
     /* Check input coordinates if angular */
     if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
@@ -39,7 +39,7 @@ pj_fwd(LP lp, PJ *P) {
 
     /* Do the transformation */
     xy = (*P->fwd)(lp, P);
-    if ( P->ctx->last_errno )
+    if ( proj_errno (P) )
         return err;
 
     /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
@@ -53,5 +53,6 @@ pj_fwd(LP lp, PJ *P) {
     xy.y = P->fr_meter * (xy.y + P->y0);
     /* z is not scaled since this is handled by vto_meter outside */
 
+    proj_errno_restore (P, last_errno);
     return xy;
 }
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
index 0fb2f633..be683745 100644
--- a/src/pj_fwd3d.c
+++ b/src/pj_fwd3d.c
@@ -1,4 +1,5 @@
 #define PJ_LIB__
+#include <proj.h>
 #include <projects.h>
 #include <errno.h>
 # define EPS 1.0e-12
@@ -9,6 +10,7 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
     XYZ xyz;
     XYZ err;
     double t;
+    int last_errno;
 
     /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
     err.x = err.y = err.z = HUGE_VAL;
@@ -16,8 +18,7 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
     if (0==P->fwd3d)
         return err;
 
-
-    P->ctx->last_errno = pj_errno = errno = 0;
+    last_errno = proj_errno_reset(P);
 
     /* Check input coordinates if angular */
     if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
@@ -55,5 +56,6 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
     xyz.y = P->fr_meter * (xyz.y + P->y0);
     /* z is not scaled since this is handled by vto_meter outside */
 
+    proj_errno_restore (P, last_errno);
     return xyz;
 }
diff --git a/src/pj_internal.c b/src/pj_internal.c
new file mode 100644
index 00000000..f3ca3566
--- /dev/null
+++ b/src/pj_internal.c
@@ -0,0 +1,236 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  This is primarily material originating from pj_obs_api.c,
+ *           that does not fit into the API category. Hence this pile of
+ *           tubings and fittings for PROJ.4 internal plumbing.
+ *
+ * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2017-07-05
+ *
+ ******************************************************************************
+ * 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"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+#define PJ_INTERNAL_C
+#include "proj_internal.h"
+#include <projects.h>
+#include <geodesic.h>
+
+#include <stddef.h>
+#include <stdarg.h>
+#include <errno.h>
+
+
+
+
+/* Used for zero-initializing new objects */
+const PJ_COORD proj_coord_null = {{0, 0, 0, 0}};
+const PJ_OBS   proj_obs_null = {
+    {{0, 0, 0, 0}},
+    {{0, 0, 0}},
+    0, 0
+};
+
+
+/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */
+PJ_COORD proj_coord_error (void) {
+    PJ_COORD c;
+    c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL;
+    return c;
+}
+
+PJ_OBS proj_obs_error (void) {
+    PJ_OBS obs;
+    obs.coo = proj_coord_error ();
+    obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL;
+    obs.id = obs.flags = 0;
+    return obs;
+}
+
+
+
+PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
+    if (0!=P->fwdobs) {
+        obs  =  P->fwdobs (obs, P);
+        return obs;
+    }
+    if (0!=P->fwd3d) {
+        obs.coo.xyz  =  pj_fwd3d (obs.coo.lpz, P);
+        return obs;
+    }
+    if (0!=P->fwd) {
+        obs.coo.xy  =  pj_fwd (obs.coo.lp, P);
+        return obs;
+    }
+    proj_errno_set (P, EINVAL);
+    return proj_obs_error ();
+}
+
+
+PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
+    if (0!=P->invobs) {
+        obs  =  P->invobs (obs, P);
+        return obs;
+    }
+    if (0!=P->inv3d) {
+        obs.coo.lpz  =  pj_inv3d (obs.coo.xyz, P);
+        return obs;
+    }
+    if (0!=P->inv) {
+        obs.coo.lp  =  pj_inv (obs.coo.xy, P);
+        return obs;
+    }
+    proj_errno_set (P, EINVAL);
+    return proj_obs_error ();
+}
+
+
+
+PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P) {
+    if (0!=P->fwdcoord)
+        return P->fwdcoord (coo, P);
+    if (0!=P->fwdobs) {
+        PJ_OBS obs = proj_obs_null;
+        obs.coo = coo;
+        obs  =  P->fwdobs (obs, P);
+        return obs.coo;
+    }
+    if (0!=P->fwd3d) {
+        coo.xyz  =  pj_fwd3d (coo.lpz, P);
+        return coo;
+    }
+    if (0!=P->fwd) {
+        coo.xy  =  pj_fwd (coo.lp, P);
+        return coo;
+    }
+    proj_errno_set (P, EINVAL);
+    return proj_coord_error ();
+}
+
+
+PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P) {
+    if (0!=P->invcoord)
+        return P->invcoord (coo, P);
+    if (0!=P->invobs) {
+        PJ_OBS obs = proj_obs_null;
+        obs.coo = coo;
+        obs  =  P->invobs (obs, P);
+        return obs.coo;
+    }
+    if (0!=P->inv3d) {
+        coo.lpz  =  pj_inv3d (coo.xyz, P);
+        return coo;
+    }
+    if (0!=P->inv) {
+        coo.lp  =  pj_inv (coo.xy, P);
+        return coo;
+    }
+    proj_errno_set (P, EINVAL);
+    return proj_coord_error ();
+}
+
+
+
+
+
+/* Move P to a new context - or to the default context if 0 is specified */
+void proj_context_set (PJ *P, PJ_CONTEXT *ctx) {
+    if (0==ctx)
+        ctx = pj_get_default_ctx ();
+    pj_set_ctx (P, ctx);
+    return;
+}
+
+void proj_context_inherit (PJ *parent, PJ *child) {
+    if (0==parent)
+        pj_set_ctx (child, pj_get_default_ctx());
+    else
+        pj_set_ctx (child, pj_get_ctx(parent));
+    return;
+}
+
+
+
+
+
+/* stuff below is *not* considered API, and will be moved to an "internal plumbing toolset" */
+
+
+
+void proj_context_errno_set (PJ_CONTEXT *ctx, int err) {
+    if (0==ctx)
+        ctx = pj_get_default_ctx();
+    pj_ctx_set_errno (ctx, err);
+    return;
+}
+
+
+
+/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
+enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level) {
+    enum proj_log_level previous;
+    if (0==ctx)
+        ctx = pj_get_default_ctx();
+    if (0==ctx)
+        return PJ_LOG_TELL;
+    previous = ctx->debug_level;
+    if (PJ_LOG_TELL==log_level)
+        return previous;
+    ctx->debug_level = log_level;
+    return previous;
+}
+
+
+
+/*  logging  */
+
+/* pj_vlog resides in pj_log.c and relates to pj_log as vsprintf relates to sprintf */
+void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args );
+
+void proj_log_error (PJ *P, const char *fmt, ...) {
+    va_list args;
+    va_start( args, fmt );
+    pj_vlog (pj_get_ctx (P), PJ_LOG_ERROR , fmt, args);
+    va_end( args );
+}
+
+void proj_log_debug (PJ *P, const char *fmt, ...) {
+    va_list args;
+    va_start( args, fmt );
+    pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MAJOR , fmt, args);
+    va_end( args );
+}
+
+void proj_log_trace (PJ *P, const char *fmt, ...) {
+    va_list args;
+    va_start( args, fmt );
+    pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MINOR , fmt, args);
+    va_end( args );
+}
+
+/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */
+void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log) {
+    if (0==ctx)
+        pj_get_default_ctx ();
+    if (0==ctx)
+        return;
+    ctx->app_data = app_data;
+    if (0!=log)
+        ctx->logger = log;
+}
diff --git a/src/pj_inv.c b/src/pj_inv.c
index d5393261..55fc917f 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -1,5 +1,6 @@
 /* general inverse projection */
 #define PJ_LIB__
+#include <proj.h>
 #include <projects.h>
 #include <errno.h>
 # define EPS 1.0e-12
@@ -8,6 +9,7 @@
 LP pj_inv(XY xy, PJ *P) {
     LP lp;
     LP err;
+    int last_errno;
 
     /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
     err.lam = err.phi = HUGE_VAL;
@@ -20,7 +22,8 @@ LP pj_inv(XY xy, PJ *P) {
         pj_ctx_set_errno( P->ctx, -15);
         return err;
     }
-    P->ctx->last_errno = errno = pj_errno = 0;
+
+    last_errno = proj_errno_reset (P);
 
     /* de-scale and de-offset */
     xy.x = (xy.x * P->to_meter - P->x0);
@@ -52,5 +55,6 @@ LP pj_inv(XY xy, PJ *P) {
             lp.phi = atan(P->one_es * tan(lp.phi));
     }
 
+    proj_errno_restore (P, last_errno);
     return lp;
-}
+}
\ No newline at end of file
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
index c6052ad4..a01cfa7e 100644
--- a/src/pj_inv3d.c
+++ b/src/pj_inv3d.c
@@ -1,4 +1,5 @@
 #define PJ_LIB__
+#include <proj.h>
 #include <projects.h>
 #include <errno.h>
 # define EPS 1.0e-12
@@ -8,6 +9,7 @@
 LPZ pj_inv3d (XYZ xyz, PJ *P) {
     LPZ lpz;
     LPZ err;
+    int last_errno;
 
     /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
     err.lam = err.phi = err.z = HUGE_VAL;
@@ -21,7 +23,7 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) {
         return err;
     }
 
-    P->ctx->last_errno = errno = pj_errno = 0;
+    last_errno = proj_errno_reset (P);
 
     /* de-scale and de-offset */
     /* z is not de-scaled since that is handled by vto_meter before we get here */
@@ -53,5 +55,6 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) {
             lpz.phi = atan(P->one_es * tan(lpz.phi));
     }
 
+    proj_errno_restore (P, last_errno);
     return lpz;
 }
diff --git a/src/pj_log.c b/src/pj_log.c
index 5b59e7a5..2525d050 100644
--- a/src/pj_log.c
+++ b/src/pj_log.c
@@ -45,9 +45,9 @@ void pj_stderr_logger( void *app_data, int level, const char *msg )
 /************************************************************************/
 /*                               pj_vlog()                              */
 /************************************************************************/
-
+void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args );
 /* Workhorse for the log functions - relates to pj_log as vsprintf relates to sprintf */
-static void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args )
+void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args )
 
 {
     char *msg_buf;
@@ -84,38 +84,3 @@ void pj_log( projCtx ctx, int level, const char *fmt, ... )
     pj_vlog( ctx, level, fmt, args );
     va_end( args );
 }
-
-
-
-/************************************************************************/
-/*             logging functions for the proj.h/obs_api                 */
-/************************************************************************/
-
-void pj_log_error (PJ *P, const char *fmt, ...)
-
-{
-    va_list args;
-    va_start( args, fmt );
-    pj_vlog (pj_get_ctx (P), PJ_LOG_ERROR , fmt, args);
-    va_end( args );
-}
-
-
-void pj_log_debug (PJ *P, const char *fmt, ...)
-
-{
-    va_list args;
-    va_start( args, fmt );
-    pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MAJOR , fmt, args);
-    va_end( args );
-}
-
-
-void pj_log_trace (PJ *P, const char *fmt, ...)
-
-{
-    va_list args;
-    va_start( args, fmt );
-    pj_vlog (pj_get_ctx (P), PJ_LOG_DEBUG_MINOR , fmt, args);
-    va_end( args );
-}
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index 37e55bc3..d8395b0f 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -34,118 +34,45 @@
  *****************************************************************************/
 #define PJ_OBS_API_C
 #include <proj.h>
-#include <projects.h>
+#include "proj_internal.h"
+#include "projects.h"
 #include <geodesic.h>
-#include <float.h>
-#include <math.h>
 
+#include <stddef.h>
+#include <errno.h>
 
-/* Used for zero-initializing new objects */
-const PJ_COORD pj_coo_null = {{0, 0, 0, 0}};
-const PJ_OBS   pj_obs_null = {
-    {{0, 0, 0, 0}},
-    {{0, 0, 0}},
-    0, 0
-};
 
-/* Magic object signaling proj system shutdown mode to routines taking a PJ * arg */
-const PJ *pj_shutdown = (PJ *) &pj_shutdown;
+PJ_COORD proj_coord (double x, double y, double z, double t) {
+    PJ_COORD res;
+    res.v[0] = x;
+    res.v[1] = y;
+    res.v[2] = z;
+    res.v[3] = t;
+    return res;
+}
 
 /* Geodesic distance between two points with angular 2D coordinates */
-double pj_lp_dist (PJ *P, LP a, LP b) {
+double proj_lp_dist (PJ *P, LP a, LP b) {
     double s12, azi1, azi2;
     /* Note: the geodesic code takes arguments in degrees */
-    geod_inverse (P->geod, TODEG(a.phi), TODEG(a.lam), TODEG(b.phi), TODEG(b.lam), &s12, &azi1, &azi2);
+    geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2);
     return s12;
 }
 
 /* Euclidean distance between two points with linear 2D coordinates */
-double pj_xy_dist (XY a, XY b) {
+double proj_xy_dist (XY a, XY b) {
     return hypot (a.x - b.x, a.y - b.y);
 }
 
 /* Euclidean distance between two points with linear 3D coordinates */
-double pj_xyz_dist (XYZ a, XYZ b) {
+double proj_xyz_dist (XYZ a, XYZ b) {
     return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z);
 }
 
 
-/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */
-static PJ_COORD pj_coo_error (void) {
-    PJ_COORD c;
-    c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL;
-    return c;
-}
-
-static PJ_OBS pj_obs_error (void) {
-    PJ_OBS obs;
-    obs.coo = pj_coo_error ();
-    obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL;
-    obs.id = obs.flags = 0;
-    return obs;
-}
-
-
-
-static PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
-    if (0!=P->fwdobs) {
-        obs  =  P->fwdobs (obs, P);
-        return obs;
-    }
-    if (0!=P->fwd3d) {
-        obs.coo.xyz  =  pj_fwd3d (obs.coo.lpz, P);
-        return obs;
-    }
-    if (0!=P->fwd) {
-        obs.coo.xy  =  pj_fwd (obs.coo.lp, P);
-        return obs;
-    }
-    pj_err_level (P, EINVAL);
-    return pj_obs_error ();
-}
-
-
-static PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
-    if (0!=P->invobs) {
-        obs  =  P->invobs (obs, P);
-        return obs;
-    }
-    if (0!=P->inv3d) {
-        obs.coo.lpz  =  pj_inv3d (obs.coo.xyz, P);
-        return obs;
-    }
-    if (0!=P->inv) {
-        obs.coo.lp  =  pj_inv (obs.coo.xy, P);
-        return obs;
-    }
-    pj_err_level (P, EINVAL);
-    return pj_obs_error ();
-}
-
-
-/* Apply the transformation P to the observation obs */
-PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs) {
-    if (0==P)
-        return obs;
-
-    switch (direction) {
-        case PJ_FWD:
-            return pj_fwdobs (obs, P);
-        case PJ_INV:
-            return  pj_invobs (obs, P);
-        case PJ_IDENT:
-            return obs;
-        default:
-            break;
-    }
-
-    pj_err_level (P, EINVAL);
-    return pj_obs_error ();
-}
-
 
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
-double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
+double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs) {
     int i;
     PJ_OBS o, u;
 
@@ -153,7 +80,7 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
         return HUGE_VAL;
 
     if (n < 1) {
-        pj_err_level (P, EINVAL);
+        proj_errno_set (P, EINVAL);
         return HUGE_VAL;
     }
 
@@ -170,116 +97,354 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
                 o  =  pj_fwdobs (u, P);
                 break;
             default:
-                pj_err_level (P, EINVAL);
+                proj_errno_set (P, EINVAL);
                 return HUGE_VAL;
         }
     }
 
-    return pj_xyz_dist (o.coo.xyz, obs.coo.xyz);
+    return proj_xyz_dist (o.coo.xyz, obs.coo.xyz);
 }
 
 
-PJ *pj_create (const char *definition) {
-    return pj_init_plus (definition);
-}
+/* Apply the transformation P to the coordinate coo */
+PJ_OBS proj_trans_obs (PJ *P, enum proj_direction direction, PJ_OBS obs) {
+    if (0==P)
+        return obs;
 
+    switch (direction) {
+        case PJ_FWD:
+            return pj_fwdobs (obs, P);
+        case PJ_INV:
+            return  pj_invobs (obs, P);
+        case PJ_IDENT:
+            return obs;
+        default:
+            break;
+    }
 
-PJ *pj_create_argv (int argc, char **argv) {
-    return pj_init (argc, argv);
+    proj_errno_set (P, EINVAL);
+    return proj_obs_error ();
 }
 
 
-/* Below: Minimum viable support for contexts. The first four functions   */
-/* relate to error reporting, debugging, and logging, hence being generically */
-/* useful. The remaining is a compact implementation of the more low level    */
-/* proj_api.h thread contexts, which may or may not be useful  */
 
-/* Set error level 0-3, or query current level */
-int pj_err_level (PJ *P, int err_level) {
-    int previous;
-    PJ_CONTEXT *ctx;
+/* Apply the transformation P to the coordinate coo */
+PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coo) {
     if (0==P)
-        ctx = pj_get_default_ctx();
-    else
-        ctx = pj_get_ctx (P);
-    previous = pj_ctx_get_errno (ctx);
-    if (PJ_ERR_TELL==err_level)
-        return previous;
-    pj_ctx_set_errno (ctx, err_level);
-    return previous;
-}
+        return coo;
 
+    switch (direction) {
+        case PJ_FWD:
+            return pj_fwdcoord (coo, P);
+        case PJ_INV:
+            return  pj_invcoord (coo, P);
+        case PJ_IDENT:
+            return coo;
+        default:
+            break;
+    }
 
-/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
-enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level) {
-    enum pj_log_level previous;
-    PJ_CONTEXT *ctx;
-    if (0==P)
-        ctx = pj_get_default_ctx();
-    else
-        ctx = pj_get_ctx (P);
-    previous = ctx->debug_level;
-    if (PJ_LOG_TELL==log_level)
-        return previous;
-    ctx->debug_level = log_level;
-    return previous;
+    proj_errno_set (P, EINVAL);
+    return proj_coord_error ();
 }
 
 
-/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */
-void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
-    PJ_CONTEXT *ctx = pj_get_ctx (P);
-    ctx->app_data = app_data;
-    if (0!=log)
-        ctx->logger = log;
-}
-
 
-/* Move P to a new context - initially a copy of the default context */
-int pj_context_renew (PJ *P) {
-    PJ_CONTEXT *ctx = pj_ctx_alloc ();
-    if (0==ctx) {
-        pj_err_level (P, ENOMEM);
-        return 1;
+/*************************************************************************************/
+size_t proj_transform (
+    PJ *P,
+    enum proj_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,
+    double *t, size_t st, size_t nt
+) {
+/**************************************************************************************
+
+    Transform a series of coordinates, where the individual coordinate dimension
+    may be represented by an array that is either
+
+        1. fully populated
+        2. a null pointer and/or a length of zero, which will be treated as a
+           fully populated array of zeroes
+        3. of length one, i.e. a constant, which will be treated as a fully
+           populated array of that constant value
+
+    The strides, sx, sy, sz, st, represent the step length, in bytes, between
+    consecutive elements of the corresponding array. This makes it possible for
+    proj_transform to handle transformation of a large class of application
+    specific data structures, without necessarily understanding the data structure
+    format, as in:
+
+        typedef struct {double x, y; int quality_level; char surveyor_name[134];} XYQS;
+        XYQS survey[345];
+        double height = 23.45;
+        PJ *P = {...};
+        size_t stride = sizeof (XYQS);
+        ...
+        proj_transform (
+            P, PJ_INV, sizeof(XYQS),
+            &(survey[0].x), stride, 345,  (*  We have 345 eastings  *)
+            &(survey[0].y), stride, 345,  (*  ...and 345 northings. *)
+            &height, 1,                   (*  The height is the constant  23.45 m *)
+            0, 0                          (*  and the time is the constant 0.00 s *)
+        );
+
+    This is similar to the inner workings of the pj_transform function, but the
+    stride functionality has been generalized to work for any size of basic unit,
+    not just a fixed number of doubles.
+
+    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.
+
+    Caveat: Since proj_transform does its work *in place*, this means that even the
+    supposedly constants (i.e. length 1 arrays) will return from the call in altered
+    state. Hence, remember to reinitialize between repeated calls.
+
+    Return value: Number of transformations completed.
+
+**************************************************************************************/
+    PJ_COORD coord = proj_coord_null;
+    size_t i, nmin;
+    double null_broadcast = 0;
+    if (0==P)
+        return 0;
+
+    /* ignore lengths of null arrays */
+    if (0==x) nx = 0;
+    if (0==y) ny = 0;
+    if (0==z) nz = 0;
+    if (0==t) nt = 0;
+
+    /* and make the nullities point to some real world memory for broadcasting nulls */
+    if (0==nx) x = &null_broadcast;
+    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;
+
+    /* arrays of length 1 are constants, which we broadcast along the longer arrays */
+    /* so we need to find the length of the shortest non-unity array to figure out  */
+    /* how many coordinate pairs we must transform */
+    nmin = (nx > 1)? nx: (ny > 1)? ny: (nz > 1)? nz: (nt > 1)? nt: 1;
+    if ((nx > 1) && (nx < nmin))  nmin = nx;
+    if ((ny > 1) && (ny < nmin))  nmin = ny;
+    if ((nz > 1) && (nz < nmin))  nmin = nz;
+    if ((nt > 1) && (nt < nmin))  nmin = nt;
+
+    /* Check validity of direction flag */
+    switch (direction) {
+        case PJ_FWD:
+        case PJ_INV:
+            break;
+        case PJ_IDENT:
+            return nmin;
+        default:
+            proj_errno_set (P, EINVAL);
+            return 0;
     }
 
-    pj_set_ctx (P, ctx);
-    return 0;
+    /* Arrays of length==0 are broadcast as the constant 0               */
+    /* Arrays of length==1 are broadcast as their single value           */
+    /* Arrays of length >1 are iterated over (for the first nmin values) */
+    /* The slightly convolved incremental indexing is used due           */
+    /* to the stride, which may be any size supported by the platform    */
+    for (i = 0;  i < nmin;  i++) {
+        coord.xyzt.x = *x;
+        coord.xyzt.y = *y;
+        coord.xyzt.z = *z;
+        coord.xyzt.t = *t;
+
+        if (PJ_FWD==direction)
+            coord = pj_fwdcoord (coord, P);
+        else
+            coord = pj_invcoord (coord, P);
+
+        /* in all full length cases, we overwrite the input with the output */
+        if (nx > 1)  {
+            *x = coord.xyzt.x;
+            x =  (double *) ( ((char *) x) + sx);
+        }
+        if (ny > 1)  {
+            *y = coord.xyzt.y;
+            y =  (double *) ( ((char *) y) + sy);
+        }
+        if (nz > 1)  {
+            *z = coord.xyzt.z;
+            z =  (double *) ( ((char *) z) + sz);
+        }
+        if (nt > 1)  {
+            *t = coord.xyzt.t;
+            t =  (double *) ( ((char *) t) + st);
+        }
+    }
+    /* Last time around, we update the length 1 cases with their transformed alter egos */
+    /* ... or would we rather not? Then what about the nmin==1 case?                    */
+    /* perhaps signalling the non-array case by setting all strides to 0?               */
+    if (nx==1)
+        *x = coord.xyzt.x;
+    if (ny==1)
+        *y = coord.xyzt.y;
+    if (nz==1)
+        *z = coord.xyzt.z;
+    if (nt==1)
+        *t = coord.xyzt.t;
+    
+    return i;
 }
 
 
-/* Move daughter to mother's context */
-void pj_context_inherit (PJ *mother, PJ *daughter) {
-    pj_set_ctx (daughter, pj_get_ctx (mother));
+
+PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
+    if (0==ctx)
+        ctx = pj_get_default_ctx ();
+    return pj_init_plus_ctx (ctx, definition);
 }
 
+PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv) {
+    if (0==ctx)
+        ctx = pj_get_default_ctx ();
+    return pj_init_ctx (ctx, argc, argv);
+}
 
-void pj_context_free (const PJ *P) {
-    PJ_CONTEXT *ctx;
+PJ *proj_destroy (PJ *P) {
+    pj_free (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.
+******************************************************************************/
+    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 */
+    proj_context_errno_set (pj_get_ctx (P), err);
+    errno = err;
+    return;
+}
+
+/*****************************************************************************/
+void proj_errno_restore (PJ *P, int err) {
+/******************************************************************************
+    Reduce some mental impedance in the canonical reset/restore use case:
+    Basically, proj_errno_restore() is a synonym for proj_errno_set(),
+    but the use cases are very different (_set: indicate an error to higher
+    level user code, _restore: pass previously set error indicators in case
+    of no errors at this level).
+
+    Hence, although the inner working is identical, we provide both options,
+    to avoid some rather confusing real world code.
+
+    See usage example under proj_errno_reset ()
+******************************************************************************/
+    proj_errno_set (P, err);
+}
 
-    /* During shutdown it should be OK to free the default context - but this needs more work */
-    if (pj_shutdown==P) {
-        /* The obvious solution "pj_ctx_free (pj_get_default_ctx ())" fails */
+/*****************************************************************************/
+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 *)
+        if (proj_errno(P))
+            return;
+        (* 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;
+}
 
-    ctx = pj_get_ctx ((PJ *) P);
 
-    /* Otherwise, trying to free the default context is a no-op */
+/* Create a new context - or provide a pointer to the default context */
+PJ_CONTEXT *proj_context_create (void) {
+    return pj_ctx_alloc ();
+}
+
+
+void proj_context_destroy (PJ_CONTEXT *ctx) {
+    if (0==ctx)
+        return;
+
+    /* Trying to free the default context is a no-op (since it is statically allocated) */
     if (pj_get_default_ctx ()==ctx)
         return;
 
-    /* The common (?) case: free the context and move the PJ to the default context */
     pj_ctx_free (ctx);
-    pj_set_ctx ((PJ *) P, pj_get_default_ctx ());
 }
 
 
-/* Minimum support for fileapi - which may never have been used... */
-void pj_fileapi_set (PJ *P, void *fileapi) {
-    PJ_CONTEXT *ctx = pj_get_ctx (P);
-    ctx->fileapi = (projFileAPI *) fileapi;
+/* Build a fully expanded proj_create() compatible representation of P */
+char *proj_definition_retrieve (PJ *P) {
+    if (0==P)
+        return 0;
+    return pj_get_def(P, 0);
 }
+
+/* ...and get rid of it safely */
+void *proj_release (void *buffer) {
+    return pj_dealloc (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);}
+
+
+/* The shape of jazz to come! */
+
+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);
diff --git a/src/proj.def b/src/proj.def
index 3792d3c0..314c4385 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -89,21 +89,45 @@ EXPORTS
     geod_polygon_testpoint          @87
     geod_polygon_clear              @88
     pj_run_selftests                @89
-    pj_create                       @90
-    pj_create_argv                  @91
-    pj_trans                        @92
-    pj_roundtrip                    @95
-    pj_log_level                    @96
-    pj_err_level                    @97
-    pj_log_func                     @98
-    pj_context_renew                @99
-    pj_context_inherit             @100
-    pj_context_free                @101
-    pj_fileapi_set                 @102
-    pj_log_error                   @103
-    pj_log_debug                   @104
-    pj_log_trace                   @105
-    pj_lp_dist                     @106
-    pj_xy_dist                     @107
-    pj_xyz_dist                    @108
-    pj_find_file                   @109
+
+    proj_create                     @90
+    proj_create_argv                @91
+    proj_destroy                    @92
+
+    proj_trans_obs                  @93
+    proj_trans_coord                @94
+    proj_transform                  @95
+    proj_roundtrip                  @96
+
+    proj_coord                      @97
+    proj_coord_error                @98
+    proj_obs_error                  @99
+
+    proj_errno                     @100
+    proj_errno_set                 @101
+    proj_errno_reset               @102
+    proj_errno_restore             @103
+    proj_context_errno_set         @104
+
+    proj_context_create            @105
+    proj_context_set               @106
+    proj_context_inherit           @107
+    proj_context_destroy           @108
+
+    proj_lp_dist                   @109
+    proj_xy_dist                   @110
+    proj_xyz_dist                  @111
+
+    proj_log_level                 @112
+    proj_log_func                  @113
+    proj_log_error                 @114
+    proj_log_debug                 @115
+    proj_log_trace                 @116
+
+    proj_definition_retrieve       @117
+    proj_release                   @118
+    proj_torad                     @119
+    proj_todeg                     @120
+
+
+    pj_find_file                   @121
diff --git a/src/proj.h b/src/proj.h
index f6b3bd2a..1977acde 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -14,7 +14,7 @@
  *
  *           This is an evolving attempt at creating a re-rationalized API
  *           with primary design goals focused on sanitizing the namespaces.
- *           Hence, all symbols exposed are being moved to the pj_ namespace,
+ *           Hence, all symbols exposed are being moved to the proj_ namespace,
  *           while all data types are being moved to the PJ_ namespace.
  *
  *           Please note that this API is *orthogonal* to  the previous APIs:
@@ -49,7 +49,8 @@
  *           without having any idea about its internals.
  *
  *           (obviously, in this example, struct projCtx_t should also be
- *           renamed struct pj_ctx some day...)
+ *           renamed struct pj_ctx some day, but for backwards compatibility
+ *           it remains as-is for now).
  *
  *           THIRD, I try to eliminate implicit type punning. Hence this API
  *           introduces the PJ_OBS ("observation") data type, for generic
@@ -73,20 +74,30 @@
  *           as large as the traditional XY and LP objects, timing results
  *           have shown the overhead to be very reasonable.
  *
- *           In its current incarnation, the API is focused on compactness:
- *           Currently it consists of only nine functions.
+ *           Contexts and thread safety
+ *           --------------------------
  *
- *           Hence, due to the proj_ctx subsystem being little used, it has
- *           been replaced by a tiny set of 3 functions, making it possible,
- *           but not very convenient, to do almost everything possible with
- *           the original ctx API.
+ *           After a year of experiments (and previous experience from the
+ *           trlib transformation library) it has become clear that the
+ *           context subsystem is unavoidable in a multi-threaded world.
+ *           Hence, instead of hiding it away, we move it into the limelight,
+ *           highly recommending (but not formally requiring) the bracketing
+ *           with calls to proj_context_create(...)/proj_context_destroy() of
+ *           any code block calling PROJ.4 functions.
  *
- *           See pj_proj_test.c for an example of how to use the API.
+ *           Legacy single threaded code need not do anything, but *may*
+ *           implement a bit of future compatibility by using the backward
+ *           compatible call proj_context_create(0), which will not create
+ *           a new context, but simply provide a pointer to the default one.
+ *
+ *           See pj_obs_api_test.c for an example of how to use the API.
  *
  * Author:   Thomas Knudsen, <thokn@sdfe.dk>
+ *           Benefitting from a large number of comments and suggestions
+ *           by (primarily) Kristian Evers and Even Rouault.
  *
  ******************************************************************************
- * 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"),
@@ -106,18 +117,15 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
+/*
 #ifdef _MSC_VER
+#ifndef _USE_MATH_DEFINES
 #define _USE_MATH_DEFINES
 #endif
-#include <math.h>
-
-#include <float.h>
-#include <limits.h>
+#endif
+#include <math.h>    For M_PI */
+#include <stddef.h>  /* For size_t */
 
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <errno.h>
 
 #ifdef PROJECTS_H
 #error proj.h must be included before projects.h
@@ -132,14 +140,33 @@
 extern "C" {
 #endif
 
+/************************************************************************
+    These version numbers should be updated with every release!
+    The format of PJ_VERSION is
+
+    * Before version 4.10.0: PJ_VERSION=MNP
+      where M, N, and P are the major, minor, and patch numbers;
+      e.g., PJ_VERSION=493 for version 4.9.3.
+
+    * Version 4.10.0 and later: PJ_VERSION=MMMNNNPP later
+      where MMM, NNN, PP are the major, minor, and patch numbers.
+      The minor and patch numbers are padded with leading zeros if
+      necessary);
+      e.g., PJ_VERSION=401000 for version 4.10.0.
+************************************************************************/
+#define PROJ_VERSION_MAJOR 10
+#define PROJ_VERSION_MINOR 0
+#define PROJ_VERSION_MICRO 0
+
+#ifndef PJ_VERSION
+#define PJ_VERSION 10##000##00
+#endif
+#ifndef PROJ_VERSION
+#define PROJ_VERSION PJ_VERSION
+#endif
 
-/* We need to access the PJ_VERSION symbol from proj_api.h */
-#define PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY
-#include <proj_api.h>
-#undef  PROJ_API_INCLUDED_FOR_PJ_VERSION_ONLY
 
 extern char const pj_release[]; /* global release id string */
-extern int        pj_errno;     /* global error return code */
 
 /* first forward declare everything needed */
 
@@ -182,6 +209,9 @@ typedef struct { double lam, phi,  z; }  LPZ;
 
 /* Ancillary pairs and triplets for geodetic computations */
 
+/* easting and northing */
+typedef struct { double   e, n; }  PJ_EN;
+
 /* Degrees, minutes, and seconds */
 typedef struct { double d, m,  s; }  PJ_DMS;
 
@@ -198,6 +228,7 @@ union PJ_COORD {
     PJ_ENHT enht;
     PJ_LPZT lpzt;
     PJ_ENH  enh;
+    PJ_EN   en;
     double v[4]; /* It's just a vector */
     XYZ  xyz;
     UVW  uvw;
@@ -227,6 +258,7 @@ union PJ_PAIR {
     LP     lp;
     UV     uv;
     PJ_AF  af;
+    PJ_EN  en;
     double v[2]; /* Yes - It's really just a vector! */
 };
 
@@ -240,91 +272,79 @@ struct PJ_OBS {
 /* The context type - properly namespaced synonym for projCtx */
 struct projCtx_t;
 typedef struct projCtx_t PJ_CONTEXT;
-typedef int *PJ_FILE;
+
+
+
+/* A P I */
+
+
+/* Functionality for handling thread contexts */
+PJ_CONTEXT *proj_context_create (void);
+void proj_context_destroy    (PJ_CONTEXT *ctx);
 
 
 /* Manage the transformation definition object PJ */
-PJ  *pj_create (const char *definition);
-PJ  *pj_create_argv (int argc, char **argv);
-void pj_free (PJ *P);
-int  pj_error (PJ *P);
+PJ  *proj_create (PJ_CONTEXT *ctx, const char *definition);
+PJ  *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv);
+PJ  *proj_destroy (PJ *P);
 
 
 /* Apply transformation to observation - in forward or inverse direction */
-enum pj_direction {
+enum proj_direction {
     PJ_FWD   =  1,   /* Forward    */
     PJ_IDENT =  0,   /* Do nothing */
     PJ_INV   = -1    /* Inverse    */
 };
-PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs);
 
-/* Measure internal consistency - in forward or inverse direction */
-double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs);
 
-/* Geodesic distance between two points with angular 2D coordinates */
-double pj_lp_dist (PJ *P, LP a, LP b);
+PJ_OBS   proj_trans_obs   (PJ *P, enum proj_direction direction, PJ_OBS obs);
+PJ_COORD proj_trans_coord (PJ *P, enum proj_direction direction, PJ_COORD coord);
 
-/* Euclidean distance between two points with linear 2D coordinates */
-double pj_xy_dist (XY a, XY b);
 
-/* Euclidean distance between two points with linear 3D coordinates */
-double pj_xyz_dist (XYZ a, XYZ b);
+size_t proj_transform (
+    PJ *P,
+    enum proj_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,
+    double *t, size_t st, size_t nt
+);
 
 
-#ifndef PJ_OBS_API_C
-/* Part of MSVC workaround: Make pj_*_null look function-like for symmetry with pj_*_error */
-#define pj_coo_null(x) pj_coo_null
-#define pj_obs_null(x) pj_obs_null
-extern const PJ_COORD pj_coo_null;
-extern const PJ_OBS   pj_obs_null;
-extern const PJ      *pj_shutdown;
-#endif
-
-#ifndef TODEG
-#define TODEG(rad)  ((rad)*180.0/M_PI)
-#endif
-#ifndef TORAD
-#define TORAD(deg)  ((deg)*M_PI/180.0)
-#endif
+/* not a constructor, but an initializer */
+PJ_COORD proj_coord (double x, double y, double z, double t);
 
 
+/* Measure internal consistency - in forward or inverse direction */
+double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs);
 
+/* Geodesic distance between two points with angular 2D coordinates */
+double proj_lp_dist (PJ *P, LP a, LP b);
 
+/* Euclidean distance between two points with linear 2D coordinates */
+double proj_xy_dist (XY a, XY b);
 
+/* Euclidean distance between two points with linear 3D coordinates */
+double proj_xyz_dist (XYZ a, XYZ b);
 
 
-/* High level functionality for handling thread contexts */
-enum pj_log_level {
-    PJ_LOG_NONE  = 0,
-    PJ_LOG_ERROR = 1,
-    PJ_LOG_DEBUG = 2,
-    PJ_LOG_TRACE = 3,
-    PJ_LOG_TELL  = 4,
-    PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */
-    PJ_LOG_DEBUG_MINOR = 3  /* for proj_api.h compatibility */
-};
-
 /* Set or read error level */
-#define PJ_ERR_TELL -56789
-int pj_err_level (PJ *P, int err_level);
-
-/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
-enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level);
-
-void pj_log_error (PJ *P, const char *fmt, ...);
-void pj_log_debug (PJ *P, const char *fmt, ...);
-void pj_log_trace (PJ *P, const char *fmt, ...);
-void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));
-
-
-/* Lower level functionality for handling thread contexts */
-int  pj_context_renew (PJ *P);
-void pj_context_inherit (PJ *mother, PJ *daughter);
-void pj_context_free    (const PJ *P);
-
-/* Lowest level: Minimum support for fileapi */
-void pj_fileapi_set (PJ *P, void *fileapi);
-
+int  proj_errno (PJ *P);
+void proj_errno_set (PJ *P, int err);
+int  proj_errno_reset (PJ *P);
+void proj_errno_restore (PJ *P, int err);
+
+/* Build a fully expanded proj_create() compatible representation of P */
+char *proj_definition_retrieve (PJ *P);
+/* ...and get rid of it safely */
+void *proj_release (void *buffer);
+
+
+/* These are trivial, and while occasionaly useful in real code, primarily here to       */
+/* simplify demo code, and in acknowledgement of the proj-internal discrepancy between   */
+/* angular units expected by classical proj, and by Charles Karney's geodesics subsystem */
+double proj_torad (double angle_in_degrees);
+double proj_todeg (double angle_in_radians);
 
 #ifdef __cplusplus
 }
diff --git a/src/proj_internal.h b/src/proj_internal.h
new file mode 100644
index 00000000..da457077
--- /dev/null
+++ b/src/proj_internal.h
@@ -0,0 +1,115 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Internal plumbing for the PROJ.4 library.
+ *
+ * Author:   Thomas Knudsen, <thokn@sdfe.dk>
+ *
+ ******************************************************************************
+ * 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"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO COORD SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifdef _MSC_VER
+#ifndef _USE_MATH_DEFINES
+#define _USE_MATH_DEFINES
+#endif
+#endif
+#include <math.h>   /* For M_PI */
+
+#include <proj.h>
+
+#ifdef PROJECTS_H
+#error proj_internal.h must be included before projects.h
+#endif
+#ifdef PROJ_API_H
+#error proj_internal.h must be included before proj_api.h
+#endif
+
+#ifndef PROJ_INTERNAL_H
+#define PROJ_INTERNAL_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+#ifndef PJ_TODEG
+#define PJ_TODEG(rad)  ((rad)*180.0/M_PI)
+#endif
+#ifndef PJ_TORAD
+#define PJ_TORAD(deg)  ((deg)*M_PI/180.0)
+#endif
+
+
+
+PJ_COORD proj_coord_error (void);
+PJ_OBS   proj_obs_error (void);
+#ifndef PJ_INTERNAL_C
+extern const PJ_COORD proj_coord_null;
+extern const PJ_OBS   proj_obs_null;
+#endif
+/* Part of MSVC workaround: Make proj_*_null look function-like for symmetry with proj_*_error */
+#define proj_coord_null(x) proj_coord_null
+#define proj_obs_null(x) proj_obs_null
+
+
+void proj_context_errno_set (PJ_CONTEXT *ctx, int err);
+void proj_context_set (PJ *P, PJ_CONTEXT *ctx);
+void proj_context_inherit (PJ *parent, PJ *child);
+
+
+PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P);
+PJ_OBS pj_invobs (PJ_OBS obs, PJ *P);
+PJ_COORD pj_fwdcoord (PJ_COORD coo, PJ *P);
+PJ_COORD pj_invcoord (PJ_COORD coo, PJ *P);
+
+
+
+/* High level functionality for handling thread contexts */
+enum proj_log_level {
+    PJ_LOG_NONE  = 0,
+    PJ_LOG_ERROR = 1,
+    PJ_LOG_DEBUG = 2,
+    PJ_LOG_TRACE = 3,
+    PJ_LOG_TELL  = 4,
+    PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */
+    PJ_LOG_DEBUG_MINOR = 3  /* for proj_api.h compatibility */
+};
+
+/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
+enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level);
+typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *);
+
+void proj_log_error (PJ *P, const char *fmt, ...);
+void proj_log_debug (PJ *P, const char *fmt, ...);
+void proj_log_trace (PJ *P, const char *fmt, ...);
+/*void proj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));*/
+void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log);
+
+
+/* Lowest level: Minimum support for fileapi */
+void proj_fileapi_set (PJ *P, void *fileapi);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* ndef PROJ_INTERNAL_H */
diff --git a/src/projects.h b/src/projects.h
index c1bf1b3a..1c99506f 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -175,6 +175,7 @@ typedef struct { double u, v, w; }     UVW;
 /* Forward declarations and typedefs for stuff needed inside the PJ object */
 struct PJconsts;
 struct PJ_OBS;
+union  PJ_COORD;
 struct geod_geodesic;
 struct pj_opaque;
 struct ARG_list;
@@ -190,6 +191,7 @@ enum pj_io_units {
 #ifndef PROJ_H
 typedef struct PJconsts PJ;         /* the PJ object herself */
 typedef struct PJ_OBS PJ_OBS;
+typedef union  PJ_COORD PJ_COORD;
 #endif
 
 struct PJ_REGION_S {
@@ -208,11 +210,12 @@ struct PJconsts {
 
     /*************************************************************************************
 
-                            G E N E R A L   C O N T E X T
+                         G E N E R A L   P A R A M E T E R   S T R U C T
 
     **************************************************************************************
 
         TODO: Need some description here - especially about the thread context...
+        This is the struct behind the PJ typedef
 
     **************************************************************************************/
 
@@ -245,11 +248,67 @@ struct PJconsts {
     LPZ (*inv3d)(XYZ, PJ *);
     PJ_OBS (*fwdobs)(PJ_OBS, PJ *);
     PJ_OBS (*invobs)(PJ_OBS, PJ *);
+    PJ_COORD (*fwdcoord)(PJ_COORD, PJ *);
+    PJ_COORD (*invcoord)(PJ_COORD, PJ *);
 
     void (*spc)(LP, PJ *, struct FACTORS *);
 
     void (*pfree)(PJ *);
 
+    /*************************************************************************************
+
+                              E R R O R   R E P O R T I N G
+
+    **************************************************************************************
+
+        Currently, we're doing error reporting through the context->last_errno indicator.
+
+        It is, however, not entirely sure this will be the right way to do it in all
+        cases: During allocation/initialization, it is certainly nice to have a higher
+        level error indicator, since we primarily signal "something went wrong", by
+        returning 0 from the pj_init family of functions - and with a null return we
+        cannot pass messages through internal state in the PJ object.
+
+        Historically, the errno variable has been used for that kind of messages, but
+        apparently, thread safety was added to PROJ.4 at a time where it was not clear
+        that errno is actually thread local.
+
+        Additionally, errno semantics has historically been misinterpreted in parts of
+        pj_init.c, a misinterpretation, that was mitigated by a hack in pj_malloc.c
+        some 15 years ago.
+
+        This PJ-local errno is a first step towards a more structured approach to
+        error reporting, being implemented in the OBS_API (cf. pj_obs_api.[ch]), and
+        related plumbing efforts.
+        
+        In due course this will let us get rid of the pj_malloc.c hack, and allow us
+        to introduce a more layered error reporting structure, where errors are
+        reported where they occur, and bubble up to the higher levels (context->errno,
+        then thread local errno), so the highest level indicate "something went wrong
+        somewhere", then the more localized ones can be used for pinpointing: 
+        
+        errno:               "something went wrong somewhere on this thread",
+        context->last_errno: "something went wrong in some PROJ.4 related code",
+        PJ->last_errno:      "It was in THIS PJ something went wrong",
+        pj_strerrno:         "This was what went wrong".
+
+        Which will be quite helpful, once fully implemented, especially for 
+        debugging complex transformation pipelines, while still maintaining backward
+        compatibility in the messaging system.
+
+        Note that there is even a global pj_errno, which is here and there accessed
+        without acquiring lock. This, and the practise of resetting the thread local
+        errno, should be given some consideration during the cleanup of the error
+        reporting system.
+
+        The name "last_errno", rather than "errno" is used partially for alignment
+        with the context->last_errno, partially because in a multithreaded environment,
+        errno is a macro, and spurious spaces turning "errno" into a separate token
+        will expose it to macro expansion, to the tune of much confusion and agony.
+
+    **************************************************************************************/
+    int last_errno;
+
 
     /*************************************************************************************
 
@@ -345,6 +404,11 @@ struct PJconsts {
 
                   D A T U M S   A N D   H E I G H T   S Y S T E M S
 
+    **************************************************************************************
+
+        It may be possible, and meaningful, to move the list parts of this up to the
+        PJ_CONTEXT level.
+
     **************************************************************************************/
 
     int     datum_type;                /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
@@ -374,9 +438,15 @@ struct PJconsts {
     struct _pj_gi *last_after_grid;     /* TODO: Description needed */
     PJ_Region     last_after_region;    /* TODO: Description needed */
     double        last_after_date;      /* TODO: Description needed */
-};
 
 
+    /*************************************************************************************
+
+                 E N D   O F    G E N E R A L   P A R A M E T E R   S T R U C T
+
+    **************************************************************************************/
+
+};
 
 
 
@@ -389,7 +459,6 @@ struct ARG_list {
 };
 
 
-
 typedef union { double  f; int  i; char *s; } PROJVALUE;
 
 
@@ -404,6 +473,7 @@ struct PJ_ELLPS {
     char    *ell;          /* elliptical parameter */
     char    *name;         /* comments */
 };
+
 struct PJ_UNITS {
     char    *id;           /* units keyword */
     char    *to_meter;     /* multiply by value to get meters */
@@ -490,8 +560,6 @@ extern struct PJ_LIST pj_list[];
 extern struct PJ_SELFTEST_LIST pj_selftest_list[];
 #endif
 
-
-
 #ifndef PJ_ELLPS__
 extern struct PJ_ELLPS pj_ellps[];
 #endif