Merge remote-tracking branch 'osgeo/master' into docs-release-4.10.0

1 files changed, 389 insertions, 213 deletions

diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 5f4bf334..fb20978b 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -28,7 +28,6 @@
  *****************************************************************************/
 #include <stddef.h>
 #include <errno.h>
-#include <ctype.h>
 #include <proj.h>
 #include "proj_internal.h"
 #include "projects.h"
@@ -53,8 +52,8 @@ int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
     dir: {PJ_FWD, PJ_INV}
 ******************************************************************************/
     if (PJ_FWD==dir)
-        return pj_left (P)==PJ_IO_UNITS_RADIANS;
-    return pj_right (P)==PJ_IO_UNITS_RADIANS;
+        return pj_left (P)==PJ_IO_UNITS_ANGULAR;
+    return pj_right (P)==PJ_IO_UNITS_ANGULAR;
 }
 
 /*****************************************************************************/
@@ -68,30 +67,47 @@ int proj_angular_output (PJ *P, enum PJ_DIRECTION dir) {
 }
 
 
+/* Geodesic distance (in meter) + fwd and rev azimuth between two points on the ellipsoid */
+PJ_COORD proj_geod (const PJ *P, PJ_COORD a, PJ_COORD b) {
+    PJ_COORD c;
+    /* Note: the geodesic code takes arguments in degrees */
+    geod_inverse (P->geod,
+        PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
+        PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
+        c.v, c.v+1, c.v+2
+    );
+
+    return c;
+}
+
+
 /* Geodesic distance (in meter) between two points with angular 2D coordinates */
-double proj_lp_dist (const PJ *P, LP a, LP b) {
+double proj_lp_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
     double s12, azi1, azi2;
     /* Note: the geodesic code takes arguments in degrees */
-    geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2);
+    geod_inverse (P->geod,
+        PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
+        PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
+        &s12, &azi1, &azi2
+    );
     return s12;
 }
 
 /* The geodesic distance AND the vertical offset */
-double proj_lpz_dist (const PJ *P, LPZ a, LPZ b) {
-    PJ_COORD aa, bb;
-    aa.lpz = a;
-    bb.lpz = b;
-    return hypot (proj_lp_dist (P, aa.lp, bb.lp), a.z - b.z);
+double proj_lpz_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
+    if (HUGE_VAL==a.lpz.lam || HUGE_VAL==b.lpz.lam)
+        return HUGE_VAL;
+    return hypot (proj_lp_dist (P, a, b), a.lpz.z - b.lpz.z);
 }
 
 /* Euclidean distance between two points with linear 2D coordinates */
-double proj_xy_dist (XY a, XY b) {
-    return hypot (a.x - b.x, a.y - b.y);
+double proj_xy_dist (PJ_COORD a, PJ_COORD b) {
+    return hypot (a.xy.x - b.xy.x, a.xy.y - b.xy.y);
 }
 
 /* Euclidean distance between two points with linear 3D coordinates */
-double proj_xyz_dist (XYZ a, XYZ b) {
-    return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z);
+double proj_xyz_dist (PJ_COORD a, PJ_COORD b) {
+    return hypot (proj_xy_dist (a, b), a.xyz.z - b.xyz.z);
 }
 
 
@@ -99,7 +115,7 @@ double proj_xyz_dist (XYZ a, XYZ b) {
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
 double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     int i;
-    PJ_COORD o, u, org;
+    PJ_COORD t, org;
 
     if (0==P)
         return HUGE_VAL;
@@ -110,29 +126,36 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     }
 
     /* in the first half-step, we generate the output value */
-    u = org = *coo;
-    o = *coo = proj_trans (P, direction, u);
+    org  = *coo;
+    *coo = proj_trans (P, direction, org);
+    t = *coo;
 
-    /* now we take n-1 full steps */
-    for (i = 0;  i < n - 1;  i++) {
-        u = proj_trans (P, -direction, o);
-        o = proj_trans (P,  direction, u);
-    }
+    /* now we take n-1 full steps in inverse direction: We are */
+    /* out of phase due to the half step already taken */
+    for (i = 0;  i < n - 1;  i++)
+        t = proj_trans (P,  direction,  proj_trans (P, -direction, t) );
 
     /* finally, we take the last half-step */
-    u = proj_trans (P, -direction, o);
+    t = proj_trans (P, -direction, t);
 
     /* checking for angular *input* since we do a roundtrip, and end where we begin */
     if (proj_angular_input (P, direction))
-        return proj_lpz_dist (P, org.lpz, u.lpz);
+        return proj_lpz_dist (P, org, t);
 
-    return proj_xyz_dist (org.xyz, u.xyz);
+    return proj_xyz_dist (org, t);
 }
 
 
 
-/* Apply the transformation P to the coordinate coo */
+/**************************************************************************************/
 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/***************************************************************************************
+Apply the transformation P to the coordinate coo, preferring the 4D interfaces if
+available.
+
+See also pj_approx_2D_trans and pj_approx_3D_trans in pj_internal.c, which work
+similarly, but prefers the 2D resp. 3D interfaces if available.
+***************************************************************************************/
     if (0==P)
         return coo;
     if (P->inverted)
@@ -334,13 +357,13 @@ size_t proj_trans_generic (
 
 
 /*************************************************************************************/
-PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+PJ_COORD proj_geocentric_latitude (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
 /**************************************************************************************
     Convert geographical latitude to geocentric (or the other way round if
     direction = PJ_INV)
 
     The conversion involves a call to the tangent function, which goes through the
-    roof at the poles, so very close (the last few micrometers) to the poles no
+    roof at the poles, so very close (the last centimeter) to the poles no
     conversion takes place and the input latitude is copied directly to the output.
 
     Fortunately, the geocentric latitude converges to the geographical at the
@@ -349,7 +372,7 @@ PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
     For the spherical case, the geographical latitude equals the geocentric, and
     consequently, the input is copied directly to the output.
 **************************************************************************************/
-    const double limit = M_HALFPI - 1e-12;
+    const double limit = M_HALFPI - 1e-9;
     PJ_COORD res = coo;
     if ((coo.lp.phi > limit) || (coo.lp.phi < -limit) || (P->es==0))
         return res;
@@ -372,6 +395,131 @@ char*  proj_rtodms(char *s, double r, int pos, int neg) {
     return rtodms(s, r, pos, neg);
 }
 
+/*************************************************************************************/
+static PJ* skip_prep_fin(PJ *P) {
+/**************************************************************************************
+Skip prepare and finalize function for the various "helper operations" added to P when
+in cs2cs compatibility mode.
+**************************************************************************************/
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+    return P;
+}
+
+/*************************************************************************************/
+static int pj_cs2cs_emulation_setup (PJ *P) {
+/**************************************************************************************
+If any cs2cs style modifiers are given (axis=..., towgs84=..., ) create the 4D API
+equivalent operations, so the preparation and finalization steps in the pj_inv/pj_fwd
+invocators can emulate the behaviour of pj_transform and the cs2cs app.
+
+Returns 1 on success, 0 on failure
+**************************************************************************************/
+    PJ *Q;
+    paralist *p;
+    if (0==P)
+        return 0;
+
+    /* Don't recurse when calling proj_create (which calls us back) */
+    if (pj_param_exists (P->params, "break_cs2cs_recursion"))
+        return 1;
+
+    /* Swap axes? */
+    p = pj_param_exists (P->params, "axis");
+
+    /* Don't axisswap if data are already in "enu" order */
+    if (p && (0!=strcmp ("enu", p->param))) {
+        char *def = malloc (100+strlen(P->axis));
+        if (0==def)
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=axisswap  axis=%s", P->axis);
+        Q = proj_create (P->ctx, def);
+        free (def);
+        if (0==Q)
+            return 0;
+        P->axisswap = skip_prep_fin(Q);
+    }
+
+    /* Geoid grid(s) given? */
+    p = pj_param_exists (P->params, "geoidgrids");
+    if (p  &&  strlen (p->param) > strlen ("geoidgrids=")) {
+        char *gridnames = p->param + strlen ("geoidgrids=");
+        char *def = malloc (100+strlen(gridnames));
+        if (0==def)
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=vgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        free (def);
+        if (0==Q)
+            return 0;
+        P->vgridshift = skip_prep_fin(Q);
+    }
+
+    /* Datum shift grid(s) given? */
+    p = pj_param_exists (P->params, "nadgrids");
+    if (p  &&  strlen (p->param) > strlen ("nadgrids=")) {
+        char *gridnames = p->param + strlen ("nadgrids=");
+        char *def = malloc (100+strlen(gridnames));
+        if (0==def)
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=hgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        free (def);
+        if (0==Q)
+            return 0;
+        P->hgridshift = skip_prep_fin(Q);
+    }
+
+    /* We ignore helmert if we have grid shift */
+    p = P->hgridshift ? 0 : pj_param_exists (P->params, "towgs84");
+    while (p) {
+        char *def;
+        char *s = p->param;
+        double *d = P->datum_params;
+        size_t n = strlen (s);
+
+        /* We ignore null helmert shifts (common in auto-translated resource files, e.g. epsg) */
+        if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6])
+            break;
+
+        if (n <= 8) /* 8==strlen ("towgs84=") */
+            return 0;
+
+        def = malloc (100+n);
+        if (0==def)
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=helmert %s transpose", s);
+        Q = proj_create (P->ctx, def);
+        pj_inherit_ellipsoid_def (P, Q);
+        free (def);
+        if (0==Q)
+            return 0;
+        P->helmert = skip_prep_fin (Q);
+
+        break;
+    }
+
+    /* We also need cartesian/geographical transformations if we are working in */
+    /* geocentric/cartesian space or we need to do a Helmert transform.         */
+    if (P->is_geocent || P->helmert) {
+        char def[150];
+        sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  f=%40.20g", P->a, P->f);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->cart = skip_prep_fin (Q);
+
+        sprintf (def, "break_cs2cs_recursion     proj=cart  ellps=WGS84");
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->cart_wgs84 = skip_prep_fin (Q);
+    }
+
+    return 1;
+}
 
 
 
@@ -386,80 +534,40 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
     It may even use free formatting "proj  =  utm;  zone  =32  ellps= GRS80".
     Note that the semicolon separator is allowed, but not required.
 **************************************************************************************/
-    PJ   *P;
-    char *args, **argv;
-    int	  argc, i, j, last, n;
+    PJ    *P;
+    char  *args, **argv;
+    size_t argc, n;
+    int    ret;
 
     if (0==ctx)
         ctx = pj_get_default_ctx ();
 
     /* Make a copy that we can manipulate */
-    n = (int) strlen (definition);
+    n = strlen (definition);
     args = (char *) malloc (n + 1);
     if (0==args)
         return 0;
     strcpy (args, definition);
 
-    /* All-in-one: count args, eliminate superfluous whitespace, 0-terminate substrings */
-    for (i = j = argc = last = 0;  i < n;  ) {
-
-        /* Skip prefix whitespace */
-        while (isspace (args[i]))
-            i++;
-
-        /* Skip at most one prefix '+' */
-        if ('+'==args[i])
-            i++;
-
-        /* Whitespace after a '+' is a syntax error - but by Postel's prescription, we ignore and go on */
-        if (isspace (args[i]))
-            continue;
-
-        /* Move a whitespace delimited text string to the left, skipping over superfluous whitespace */
-        while ((0!=args[i]) && (!isspace (args[i])) && (';'!=args[i]))
-            args[j++] = args[i++];
-
-        /* Skip postfix whitespace */
-        while (isspace (args[i]) || ';'==args[i])
-            i++;
-
-        /* Greedy assignment operator: turn "a = b" into "a=b" */
-        if ('='==args[i]) {
-            args[j++] = '=';
-            i++;
-            while (isspace (args[i]))
-                i++;
-            while ((0!=args[i]) && (!isspace (args[i])) && (';'!=args[i]))
-                args[j++] = args[i++];
-            while (isspace (args[i]) || ';'==args[i])
-                i++;
-        }
-
-        /* terminate string - if that makes j pass i (often the case for first arg), let i catch up */
-        args[j++] = 0;
-        if (i < j)
-            i = j;
-
-        /* we finished another arg */
-        argc++;
+    argc = pj_trim_argc (args);
+    if (argc==0) {
+        pj_dealloc (args);
+        return 0;
     }
 
-    /* turn the massaged input into an array of strings */
-    argv = (char **) calloc (argc, sizeof (char *));
-    if (0==argv)
-        return pj_dealloc (args);
-    argv[0] = args;
-    for (i = 0, j = 1;  i < n;  i++) {
-        if (0==args[i])
-            argv[j++] = args + (i + 1);
-        if (j==argc)
-            break;
-    }
+    argv = pj_trim_argv (argc, args);
 
     /* ...and let pj_init_ctx do the hard work */
-    P = pj_init_ctx (ctx, argc, argv);
+    P = pj_init_ctx (ctx, (int) argc, argv);
+
     pj_dealloc (argv);
     pj_dealloc (args);
+
+    /* Support cs2cs-style modifiers */
+    ret = pj_cs2cs_emulation_setup  (P);
+    if (0==ret)
+        return proj_destroy (P);
+
     return P;
 }
 
@@ -474,11 +582,23 @@ a null-pointer is returned. The definition arguments may use '+' as argument sta
 indicator, as in {"+proj=utm", "+zone=32"}, or leave it out, as in {"proj=utm",
 "zone=32"}.
 **************************************************************************************/
+    PJ *P;
+    const char *c;
+
     if (0==argv)
         return 0;
     if (0==ctx)
         ctx = pj_get_default_ctx ();
-    return pj_init_ctx (ctx, argc, argv);
+
+    /* We assume that free format is used, and build a full proj_create compatible string */
+    c = pj_make_args (argc, argv);
+    if (0==c)
+        return 0;
+
+    P = proj_create (ctx, c);
+
+    pj_dealloc ((char *) c);
+    return P;
 }
 
 
@@ -531,11 +651,26 @@ PJ *proj_destroy (PJ *P) {
     return 0;
 }
 
+/*****************************************************************************/
 int proj_errno (const PJ *P) {
+/******************************************************************************
+    Read an error level from the context of a PJ.
+******************************************************************************/
     return pj_ctx_get_errno (pj_get_ctx ((PJ *) P));
 }
 
 /*****************************************************************************/
+int proj_context_errno (PJ_CONTEXT *ctx) {
+/******************************************************************************
+    Read an error directly from a context, without going through a PJ
+    belonging to that context.
+******************************************************************************/
+    if (0==ctx)
+        ctx = pj_get_default_ctx();
+    return pj_ctx_get_errno (ctx);
+}
+
+/*****************************************************************************/
 int proj_errno_set (const PJ *P, int err) {
 /******************************************************************************
     Set context-errno, bubble it up to the thread local errno, return err
@@ -599,6 +734,7 @@ int proj_errno_reset (const PJ *P) {
 
     pj_ctx_set_errno (pj_get_ctx ((PJ *) P), 0);
     errno = 0;
+    pj_errno = 0;
     return last_errno;
 }
 
@@ -622,70 +758,113 @@ PJ_CONTEXT *proj_context_destroy (PJ_CONTEXT *ctx) {
 }
 
 
+
+
+
+
 /*****************************************************************************/
-PJ_INFO proj_info(void) {
+static char *path_append (char *buf, const char *app, size_t *buf_size) {
 /******************************************************************************
-    Basic info about the current instance of the PROJ.4 library.
+    Helper for proj_info() below. Append app to buf, separated by a
+    semicolon. Also handle allocation of longer buffer if needed.
+
+    Returns buffer and adjusts *buf_size through provided pointer arg.
+******************************************************************************/
+    char *p;
+    size_t len, applen = 0, buflen = 0;
+#ifdef _WIN32
+    char *delim = ";";
+#else
+    char *delim = ":";
+#endif
+
+    /* Nothing to do? */
+    if (0 == app)
+        return buf;
+    applen = strlen (app);
+    if (0 == applen)
+        return buf;
+
+    /* Start checking whether buf is long enough */
+    if (0 != buf)
+        buflen = strlen (buf);
+    len = buflen+applen+strlen (delim) + 1;
+
+    /* "pj_realloc", so to speak */
+    if (*buf_size < len) {
+        p = pj_calloc (2 * len, sizeof (char));
+        if (0==p) {
+            pj_dealloc (buf);
+            return 0;
+        }
+        *buf_size = 2 * len;
+        if (buf != 0)
+            strcpy (p, buf);
+        pj_dealloc (buf);
+        buf = p;
+    }
 
-    Returns PJ_INFO struct. Searchpath member of the struct is truncated to 512
-    characters.
+    /* Only append a semicolon if something's already there */
+    if (0 != buflen)
+        strcat (buf, ";");
+    strcat (buf, app);
+    return buf;
+}
 
+static const char *empty = {""};
+static char version[64]  = {""};
+static PJ_INFO info = {0, 0, 0, 0, 0, 0, 0, 0};
+static volatile int info_initialized = 0;
+
+/*****************************************************************************/
+PJ_INFO proj_info (void) {
+/******************************************************************************
+    Basic info about the current instance of the PROJ.4 library.
+
+    Returns PJ_INFO struct.
 ******************************************************************************/
-    PJ_INFO info;
-    const char **paths;
-    char *tmpstr;
-    int i, n;
-    size_t len = 0;
+    const char * const *paths;
+    size_t i, n;
+
+    size_t  buf_size = 0;
+    char   *buf = 0;
 
-    memset(&info, 0, sizeof(PJ_INFO));
+    pj_acquire_lock ();
+
+    if (0!=info_initialized) {
+        pj_release_lock ();
+        return info;
+    }
 
     info.major = PROJ_VERSION_MAJOR;
     info.minor = PROJ_VERSION_MINOR;
     info.patch = PROJ_VERSION_PATCH;
 
     /* This is a controlled environment, so no risk of sprintf buffer
-       overflow. A normal version string is xx.yy.zz which is 8 characters
-       long and there is room for 64 bytes in the version string. */
-    sprintf(info.version, "%d.%d.%d", info.major, info.minor, info.patch);
-
-    pj_strlcpy(info.release, pj_get_release(), sizeof(info.release));
+    overflow. A normal version string is xx.yy.zz which is 8 characters
+    long and there is room for 64 bytes in the version string. */
+    sprintf (version, "%d.%d.%d", info.major, info.minor, info.patch);
 
+    info.searchpath = empty;
+    info.version    = version;
+    info.release    = pj_get_release ();
 
     /* build search path string */
-    tmpstr = getenv("HOME");
-    if (tmpstr != NULL) {
-        pj_strlcpy(info.searchpath, tmpstr, sizeof(info.searchpath));
-    }
-
-    tmpstr = getenv("PROJ_LIB");
-    if (tmpstr != NULL) {
-        if (strlen(info.searchpath) != 0) {
-            /* $HOME already in path */
-            strcat(info.searchpath, ";");
-            len = strlen(tmpstr);
-            strncat(info.searchpath, tmpstr, sizeof(info.searchpath)-len-1);
-        } else {
-            /* path is empty */
-            pj_strlcpy(info.searchpath, tmpstr, sizeof(info.searchpath));
-        }
-    }
+    buf = path_append (buf, getenv ("HOME"), &buf_size);
+    buf = path_append (buf, getenv ("PROJ_LIB"), &buf_size);
 
-    paths = proj_get_searchpath();
-    n = proj_get_path_count();
+    paths = proj_get_searchpath ();
+    n = (size_t) proj_get_path_count ();
 
-    for (i=0; i<n; i++) {
-        if (strlen(info.searchpath)+strlen(paths[i]) >= 511)
-            continue;
+    for (i = 0;  i < n;  i++)
+        buf = path_append (buf, paths[i], &buf_size);
+    info.searchpath = buf ? buf : empty;
 
-        if (strlen(info.searchpath) != 0) {
-            strcat(info.searchpath, ";");
-            len = strlen(paths[i]);
-            strncat(info.searchpath, paths[i], sizeof(info.searchpath)-len-1);
-        } else {
-            pj_strlcpy(info.searchpath, paths[i], sizeof(info.searchpath));
-        }
-    }
+    info.paths = paths;
+    info.path_count = n;
 
+    info_initialized = 1;
+    pj_release_lock ();
     return info;
 }
 
@@ -696,39 +875,41 @@ PJ_PROJ_INFO proj_pj_info(PJ *P) {
     Basic info about a particular instance of a projection object.
 
     Returns PJ_PROJ_INFO struct.
-
 ******************************************************************************/
-    PJ_PROJ_INFO info;
+    PJ_PROJ_INFO pjinfo;
     char *def;
 
-    memset(&info, 0, sizeof(PJ_PROJ_INFO));
+    memset(&pjinfo, 0, sizeof(PJ_PROJ_INFO));
 
     /* Expected accuracy of the transformation. Hardcoded for now, will be improved */
     /* later. Most likely to be used when a transformation is set up with           */
     /* proj_create_crs_to_crs in a future version that leverages the EPSG database. */
-    info.accuracy = -1.0;
+    pjinfo.accuracy = -1.0;
 
-    if (!P) {
-        return info;
-    }
+    if (0==P)
+        return pjinfo;
 
     /* projection id */
     if (pj_param(P->ctx, P->params, "tproj").i)
-        pj_strlcpy(info.id, pj_param(P->ctx, P->params, "sproj").s, sizeof(info.id));
+        pjinfo.id = pj_param(P->ctx, P->params, "sproj").s;
 
     /* projection description */
-    pj_strlcpy(info.description, P->descr, sizeof(info.description));
+    pjinfo.description = P->descr;
 
     /* projection definition */
-    def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */
-    pj_strlcpy(info.definition, &def[1], sizeof(info.definition)); /* def includes a leading space */
-    pj_dealloc(def);
-
-    /* this does not take into account that a pipeline potentially does not */
-    /* have an inverse.                                                     */
-    info.has_inverse = (P->inv != 0 || P->inv3d != 0 || P->inv4d != 0);
+    if (P->def_full)
+        def = P->def_full;
+    else
+        def = pj_get_def(P, 0); /* pj_get_def takes a non-const PJ pointer */
+    if (0==def)
+        pjinfo.definition = empty;
+    else
+        pjinfo.definition = pj_shrink (def);
+    /* Make pj_free clean this up eventually */
+    P->def_full = def;
 
-    return info;
+    pjinfo.has_inverse = pj_has_inverse(P);
+    return pjinfo;
 }
 
 
@@ -738,49 +919,53 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) {
     Information about a named datum grid.
 
     Returns PJ_GRID_INFO struct.
-
 ******************************************************************************/
-    PJ_GRID_INFO info;
+    PJ_GRID_INFO grinfo;
 
     /*PJ_CONTEXT *ctx = proj_context_create(); */
     PJ_CONTEXT *ctx = pj_get_default_ctx();
     PJ_GRIDINFO *gridinfo = pj_gridinfo_init(ctx, gridname);
-    memset(&info, 0, sizeof(PJ_GRID_INFO));
+    memset(&grinfo, 0, sizeof(PJ_GRID_INFO));
 
     /* in case the grid wasn't found */
     if (gridinfo->filename == NULL) {
         pj_gridinfo_free(ctx, gridinfo);
-        strcpy(info.format, "missing");
-        return info;
+        strcpy(grinfo.format, "missing");
+        return grinfo;
     }
 
+    /* The string copies below are automatically null-terminated due to */
+    /* the memset above, so strncpy is safe */
+
     /* name of grid */
-    pj_strlcpy(info.gridname, gridname, sizeof(info.gridname));
+    strncpy (grinfo.gridname, gridname, sizeof(grinfo.gridname) - 1);
 
     /* full path of grid */
-    pj_find_file(ctx, gridname, info.filename, sizeof(info.filename));
+    pj_find_file(ctx, gridname, grinfo.filename, sizeof(grinfo.filename) - 1);
 
     /* grid format */
-    pj_strlcpy(info.format, gridinfo->format, sizeof(info.format));
+    strncpy (grinfo.format, gridinfo->format, sizeof(grinfo.format) - 1);
 
     /* grid size */
-    info.n_lon = gridinfo->ct->lim.lam;
-    info.n_lat = gridinfo->ct->lim.phi;
+    grinfo.n_lon = gridinfo->ct->lim.lam;
+    grinfo.n_lat = gridinfo->ct->lim.phi;
 
     /* cell size */
-    info.cs_lon = gridinfo->ct->del.lam;
-    info.cs_lat = gridinfo->ct->del.phi;
+    grinfo.cs_lon = gridinfo->ct->del.lam;
+    grinfo.cs_lat = gridinfo->ct->del.phi;
 
     /* bounds of grid */
-    info.lowerleft  = gridinfo->ct->ll;
-    info.upperright.lam = info.lowerleft.lam + info.n_lon*info.cs_lon;
-    info.upperright.phi = info.lowerleft.phi + info.n_lat*info.cs_lat;
+    grinfo.lowerleft  = gridinfo->ct->ll;
+    grinfo.upperright.lam = grinfo.lowerleft.lam + grinfo.n_lon*grinfo.cs_lon;
+    grinfo.upperright.phi = grinfo.lowerleft.phi + grinfo.n_lat*grinfo.cs_lat;
 
     pj_gridinfo_free(ctx, gridinfo);
 
-    return info;
+    return grinfo;
 }
 
+
+
 /*****************************************************************************/
 PJ_INIT_INFO proj_init_info(const char *initname){
 /******************************************************************************
@@ -790,75 +975,77 @@ PJ_INIT_INFO proj_init_info(const char *initname){
 
     Returns PJ_INIT_INFO struct.
 
-    If the init file is not found all members of
-    the return struct are set to 0. If the init file is found, but it the
-    metadata is missing, the value is set to "Unknown".
+    If the init file is not found all members of the return struct are set
+    to the empty string.
 
+    If the init file is found, but the metadata is missing, the value is
+    set to "Unknown".
 ******************************************************************************/
-    int file_found, def_found=0;
+    int file_found;
     char param[80], key[74];
     paralist *start, *next;
-    PJ_INIT_INFO info;
+    PJ_INIT_INFO ininfo;
     PJ_CONTEXT *ctx = pj_get_default_ctx();
 
-    memset(&info, 0, sizeof(PJ_INIT_INFO));
+    memset(&ininfo, 0, sizeof(PJ_INIT_INFO));
 
-    file_found = pj_find_file(ctx, initname, info.filename, sizeof(info.filename));
+    file_found = pj_find_file(ctx, initname, ininfo.filename, sizeof(ininfo.filename));
     if (!file_found || strlen(initname) > 64) {
-        return info;
+        return ininfo;
     }
 
-    pj_strlcpy(info.name, initname, sizeof(info.name));
-    strcpy(info.origin, "Unknown");
-    strcpy(info.version, "Unknown");
-    strcpy(info.lastupdate, "Unknown");
+    /* The initial memset (0) makes strncpy safe here */
+    strncpy (ininfo.name, initname, sizeof(ininfo.name) - 1);
+    strcpy(ininfo.origin, "Unknown");
+    strcpy(ininfo.version, "Unknown");
+    strcpy(ininfo.lastupdate, "Unknown");
 
-    pj_strlcpy(key, initname, 64); /* make room for ":metadata\0" at the end */
+    strncpy (key, initname, 64); /* make room for ":metadata\0" at the end */
+    key[64] = 0;
     strncat(key, ":metadata", 9);
     strcpy(param, "+init=");
     strncat(param, key, 73);
 
     start = pj_mkparam(param);
-    next = pj_get_init(ctx, &start, start, key, &def_found);
+    pj_expand_init(ctx, start);
 
-    if (pj_param(ctx, start, "tversion").i) {
-        pj_strlcpy(info.version, pj_param(ctx, start, "sversion").s, sizeof(info.version));
-    }
+    if (pj_param(ctx, start, "tversion").i)
+        strncpy(ininfo.version, pj_param(ctx, start, "sversion").s, sizeof(ininfo.version) - 1);
 
-    if (pj_param(ctx, start, "torigin").i) {
-        pj_strlcpy(info.origin, pj_param(ctx, start, "sorigin").s, sizeof(info.origin));
-    }
+    if (pj_param(ctx, start, "torigin").i)
+        strncpy(ininfo.origin, pj_param(ctx, start, "sorigin").s, sizeof(ininfo.origin) - 1);
 
-    if (pj_param(ctx, start, "tlastupdate").i) {
-        pj_strlcpy(info.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(info.lastupdate));
-    }
+    if (pj_param(ctx, start, "tlastupdate").i)
+        strncpy(ininfo.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(ininfo.lastupdate) - 1);
 
     for ( ; start; start = next) {
         next = start->next;
         pj_dalloc(start);
     }
 
-   return info;
+   return ininfo;
 }
 
 
 
 /*****************************************************************************/
-PJ_FACTORS proj_factors(PJ *P, LP lp) {
+PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
 /******************************************************************************
     Cartographic characteristics at point lp.
 
     Characteristics include meridian, parallel and areal scales, angular
     distortion, meridian/parallel, meridian convergence and scale error.
 
-    returns PJ_FACTORS. If unsuccessfull error number is set and the returned
-    struct contains NULL data.
-
+    returns PJ_FACTORS. If unsuccessfull, error number is set and the
+    struct returned contains NULL data.
 ******************************************************************************/
-    PJ_FACTORS factors = {0,0,0, 0,0,0, 0,0};
+    PJ_FACTORS factors = {0,0,0,  0,0,0,  0,0,  0,0,0,0};
     struct FACTORS f;
 
-    if (pj_factors(lp, P, 0.0, &f))
+    if (0==P)
+        return factors;
+
+    if (pj_factors(lp.lp, P, 0.0, &f))
         return factors;
 
     factors.meridional_scale  =  f.h;
@@ -872,23 +1059,12 @@ PJ_FACTORS proj_factors(PJ *P, LP lp) {
     factors.tissot_semimajor  =  f.a;
     factors.tissot_semiminor  =  f.b;
 
-    return factors;
-}
+    /* Raw derivatives, for completeness's sake */
+    factors.dx_dlam = f.der.x_l;
+    factors.dx_dphi = f.der.x_p;
+    factors.dy_dlam = f.der.y_l;
+    factors.dy_dphi = f.der.y_p;
 
-
-const PJ_ELLPS *proj_list_ellps(void) {
-    return pj_get_ellps_ref();
-}
-
-const PJ_UNITS *proj_list_units(void) {
-    return pj_get_units_ref();
-}
-
-const PJ_OPERATIONS *proj_list_operations(void) {
-    return pj_get_list_ref();
-}
-
-const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) {
-    return pj_get_prime_meridians_ref();
+    return factors;
 }