cpp conversion: move source files in apps/ iso19111/ conversions/ projections/ transformations/ tests/ subdirectories

5 files changed, 1194 insertions, 0 deletions

diff --git a/src/conversions/PJ_axisswap.cpp b/src/conversions/PJ_axisswap.cpp
new file mode 100644
index 00000000..8714ec85
--- /dev/null
+++ b/src/conversions/PJ_axisswap.cpp
@@ -0,0 +1,302 @@
+/***********************************************************************
+
+        Axis order operation for use with transformation pipelines.
+
+                Kristian Evers, kreve@sdfe.dk, 2017-10-31
+
+************************************************************************
+
+Change the order and sign of 2,3 or 4 axes. Each of the possible four
+axes are numbered with 1-4, such that the first input axis is 1, the
+second is 2 and so on. The output ordering is controlled by a list of the
+input axes re-ordered to the new mapping. Examples:
+
+Reversing the order of the axes:
+
+    +proj=axisswap +order=4,3,2,1
+
+Swapping the first two axes (x and y):
+
+    +proj=axisswap +order=2,1,3,4
+
+The direction, or sign, of an axis can be changed by adding a minus in
+front of the axis-number:
+
+    +proj=axisswap +order=1,-2,3,4
+
+It is only necessary to specify the axes that are affected by the swap
+operation:
+
+    +proj=axisswap +order=2,1
+
+************************************************************************
+* Copyright (c) 2017, Kristian Evers / 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_LIB__
+#include <errno.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "proj_internal.h"
+#include "projects.h"
+
+PROJ_HEAD(axisswap, "Axis ordering");
+
+namespace { // anonymous namespace
+struct pj_opaque {
+    unsigned int axis[4];
+    int sign[4];
+};
+} // anonymous namespace
+
+static int sign(int x) {
+    return (x > 0) - (x < 0);
+}
+
+static XY forward_2d(LP lp, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD out, in;
+
+    in.lp = lp;
+    out = proj_coord_error();
+
+    for (i=0; i<2; i++)
+        out.v[i] = in.v[Q->axis[i]] * Q->sign[i];
+
+    return out.xy;
+}
+
+
+static LP reverse_2d(XY xy, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD out, in;
+
+    in.xy = xy;
+    out = proj_coord_error();
+
+    for (i=0; i<2; i++)
+        out.v[Q->axis[i]] = in.v[i] * Q->sign[i];
+
+    return out.lp;
+}
+
+
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD out, in;
+
+    in.lpz = lpz;
+    out = proj_coord_error();
+
+    for (i=0; i<3; i++)
+        out.v[i] = in.v[Q->axis[i]] * Q->sign[i];
+
+    return out.xyz;
+}
+
+static LPZ reverse_3d(XYZ xyz, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD in, out;
+
+    out = proj_coord_error();
+    in.xyz = xyz;
+
+    for (i=0; i<3; i++)
+        out.v[Q->axis[i]] = in.v[i] * Q->sign[i];
+
+    return out.lpz;
+}
+
+
+static PJ_COORD forward_4d(PJ_COORD coo, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD out;
+
+    out = proj_coord_error();
+
+    for (i=0; i<4; i++)
+        out.v[i] = coo.v[Q->axis[i]] * Q->sign[i];
+
+    return out;
+}
+
+
+static PJ_COORD reverse_4d(PJ_COORD coo, PJ *P) {
+    struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
+    unsigned int i;
+    PJ_COORD out;
+
+    out = proj_coord_error();
+
+    for (i=0; i<4; i++)
+        out.v[Q->axis[i]] = coo.v[i] * Q->sign[i];
+
+    return out;
+}
+
+
+/***********************************************************************/
+PJ *CONVERSION(axisswap,0) {
+/***********************************************************************/
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+    char *s;
+    unsigned int i, j, n = 0;
+
+    if (nullptr==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = (void *) Q;
+
+
+    /* +order and +axis are mutually exclusive */
+    if ( !pj_param_exists(P->params, "order") == !pj_param_exists(P->params, "axis") )
+        return pj_default_destructor(P, PJD_ERR_AXIS);
+
+    /* fill axis list with indices from 4-7 to simplify duplicate search further down */
+    for (i=0; i<4; i++) {
+        Q->axis[i] = i+4;
+        Q->sign[i] = 1;
+    }
+
+    /* if the "order" parameter is used */
+    if ( pj_param_exists(P->params, "order") ) {
+        /* read axis order */
+        char *order = pj_param(P->ctx, P->params, "sorder").s;
+
+        /* check that all characters are valid */
+        for (i=0; i<strlen(order); i++)
+            if (strchr("1234-,", order[i]) == nullptr) {
+                proj_log_error(P, "axisswap: unknown axis '%c'", order[i]);
+                return pj_default_destructor(P, PJD_ERR_AXIS);
+            }
+
+        /* read axes numbers and signs */
+        s = order;
+        n = 0;
+        while ( *s != '\0' && n < 4 ) {
+            Q->axis[n] = abs(atoi(s))-1;
+            if (Q->axis[n] > 3) {
+                proj_log_error(P, "axisswap: invalid axis '%d'", Q->axis[n]);
+                return pj_default_destructor(P, PJD_ERR_AXIS);
+            }
+            Q->sign[n++] = sign(atoi(s));
+            while ( *s != '\0' && *s != ',' )
+                s++;
+            if ( *s == ',' )
+                s++;
+        }
+    }
+
+    /* if the "axis" parameter is used */
+    if ( pj_param_exists(P->params, "axis") ) {
+        /* parse the classic PROJ.4 enu axis specification */
+        for (i=0; i < 3; i++) {
+            switch(P->axis[i]) {
+                case 'w':
+                    Q->sign[i] = -1;
+                    Q->axis[i] = 0;
+                    break;
+                case 'e':
+                    Q->sign[i] = 1;
+                    Q->axis[i] = 0;
+                    break;
+                case 's':
+                    Q->sign[i] = -1;
+                    Q->axis[i] = 1;
+                    break;
+                case 'n':
+                    Q->sign[i] = 1;
+                    Q->axis[i] = 1;
+                    break;
+                case 'd':
+                    Q->sign[i] = -1;
+                    Q->axis[i] = 2;
+                    break;
+                case 'u':
+                    Q->sign[i] = 1;
+                    Q->axis[i] = 2;
+                    break;
+                default:
+                    proj_log_error(P, "axisswap: unknown axis '%c'", P->axis[i]);
+                    return pj_default_destructor(P, PJD_ERR_AXIS);
+            }
+        }
+        n = 3;
+    }
+
+    /* check for duplicate axes */
+    for (i=0; i<4; i++)
+        for (j=0; j<4; j++) {
+            if (i==j)
+                continue;
+            if (Q->axis[i] == Q->axis[j]) {
+                proj_log_error(P, "swapaxis: duplicate axes specified");
+                return pj_default_destructor(P, PJD_ERR_AXIS);
+            }
+        }
+
+
+    /* only map fwd/inv functions that are possible with the given axis setup */
+    if (n == 4) {
+        P->fwd4d = forward_4d;
+        P->inv4d = reverse_4d;
+    }
+    if (n == 3 && Q->axis[0] < 3 && Q->axis[1] < 3 && Q->axis[2] < 3) {
+        P->fwd3d  = forward_3d;
+        P->inv3d  = reverse_3d;
+    }
+    if (n == 2 && Q->axis[0] < 2 && Q->axis[1] < 2) {
+        P->fwd    = forward_2d;
+        P->inv    = reverse_2d;
+    }
+
+
+    if (P->fwd4d == nullptr && P->fwd3d == nullptr && P->fwd == nullptr) {
+        proj_log_error(P, "swapaxis: bad axis order");
+        return pj_default_destructor(P, PJD_ERR_AXIS);
+    }
+
+    if (pj_param(P->ctx, P->params, "tangularunits").i) {
+        P->left  = PJ_IO_UNITS_ANGULAR;
+        P->right = PJ_IO_UNITS_ANGULAR;
+    } else {
+        P->left  = PJ_IO_UNITS_WHATEVER;
+        P->right = PJ_IO_UNITS_WHATEVER;
+    }
+
+
+    /* Preparation and finalization steps are skipped, since the raison   */
+    /* d'etre of axisswap is to bring input coordinates in line with the  */
+    /* the internally expected order (ENU), such that handling of offsets */
+    /* etc. can be done correctly in a later step of a pipeline */
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+
+    return P;
+}
diff --git a/src/conversions/PJ_cart.cpp b/src/conversions/PJ_cart.cpp
new file mode 100644
index 00000000..6fed9985
--- /dev/null
+++ b/src/conversions/PJ_cart.cpp
@@ -0,0 +1,219 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Convert between ellipsoidal, geodetic coordinates and
+ *           cartesian, geocentric coordinates.
+ *
+ *           Formally, this functionality is also found in the PJ_geocent.c
+ *           code.
+ *
+ *           Actually, however, the PJ_geocent transformations are carried
+ *           out in concert between 2D stubs in PJ_geocent.c and 3D code
+ *           placed in pj_transform.c.
+ *
+ *           For pipeline-style datum shifts, we do need direct access
+ *           to the full 3D interface for this functionality.
+ *
+ *           Hence this code, which may look like "just another PJ_geocent"
+ *           but really is something substantially different.
+ *
+ * Author:   Thomas Knudsen, thokn@sdfe.dk
+ *
+ ******************************************************************************
+ * Copyright (c) 2016, 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_LIB__
+
+#include "proj_internal.h"
+#include "projects.h"
+#include "proj_math.h"
+
+PROJ_HEAD(cart,    "Geodetic/cartesian conversions");
+
+
+/**************************************************************
+                CARTESIAN / GEODETIC CONVERSIONS
+***************************************************************
+    This material follows:
+
+    Bernhard Hofmann-Wellenhof & Helmut Moritz:
+    Physical Geodesy, 2nd edition.
+    Springer, 2005.
+
+    chapter 5.6: Coordinate transformations
+    (HM, below),
+
+    and
+
+    Wikipedia: Geographic Coordinate Conversion,
+    https://en.wikipedia.org/wiki/Geographic_coordinate_conversion
+
+    (WP, below).
+
+    The cartesian-to-geodetic conversion is based on Bowring's
+    celebrated method:
+
+    B. R. Bowring:
+    Transformation from spatial to geographical coordinates
+    Survey Review 23(181), pp. 323-327, 1976
+
+    (BB, below),
+
+    but could probably use some TLC from a newer and faster
+    algorithm:
+
+    Toshio Fukushima:
+    Transformation from Cartesian to Geodetic Coordinates
+    Accelerated by Halley’s Method
+    Journal of Geodesy, February 2006
+
+    (TF, below).
+
+    Close to the poles, we avoid singularities by switching to an
+    approximation requiring knowledge of the geocentric radius
+    at the given latitude. For this, we use an adaptation of the
+    formula given in:
+
+    Wikipedia: Earth Radius
+    https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude
+    (Derivation and commentary at https://gis.stackexchange.com/q/20200)
+
+    (WP2, below)
+
+    These routines are probably not as robust at those in
+    geocent.c, at least thay haven't been through as heavy
+    use as their geocent sisters. Some care has been taken
+    to avoid singularities, but extreme cases (e.g. setting
+    es, the squared eccentricity, to 1), will cause havoc.
+
+**************************************************************/
+
+
+/*********************************************************************/
+static double normal_radius_of_curvature (double a, double es, double phi) {
+/*********************************************************************/
+    double s = sin(phi);
+    if (es==0)
+        return a;
+    /* This is from WP.  HM formula 2-149 gives an a,b version */
+    return a / sqrt (1 - es*s*s);
+}
+
+/*********************************************************************/
+static double geocentric_radius (double a, double b, double phi) {
+/*********************************************************************
+    Return the geocentric radius at latitude phi, of an ellipsoid
+    with semimajor axis a and semiminor axis b.
+
+    This is from WP2, but uses hypot() for potentially better
+    numerical robustness
+***********************************************************************/
+    return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi));
+}
+
+
+/*********************************************************************/
+static XYZ cartesian (LPZ geod,  PJ *P) {
+/*********************************************************************/
+    double N, cosphi = cos(geod.phi);
+    XYZ xyz;
+
+    N   =  normal_radius_of_curvature(P->a, P->es, geod.phi);
+
+    /* HM formula 5-27 (z formula follows WP) */
+    xyz.x = (N + geod.z) * cosphi      * cos(geod.lam);
+    xyz.y = (N + geod.z) * cosphi      * sin(geod.lam);
+    xyz.z = (N * (1 - P->es) + geod.z) * sin(geod.phi);
+
+    return xyz;
+}
+
+
+/*********************************************************************/
+static LPZ geodetic (XYZ cart,  PJ *P) {
+/*********************************************************************/
+    double N, p, theta, c, s;
+    LPZ lpz;
+
+    /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */
+    p = hypot (cart.x, cart.y);
+
+    /* HM eq. (5-37) */
+    theta  =  atan2 (cart.z * P->a,  p * P->b);
+
+    /* HM eq. (5-36) (from BB, 1976) */
+    c  =  cos(theta);
+    s  =  sin(theta);
+    lpz.phi  =  atan2 (cart.z + P->e2s*P->b*s*s*s,  p - P->es*P->a*c*c*c);
+    lpz.lam  =  atan2 (cart.y, cart.x);
+    N        =  normal_radius_of_curvature (P->a, P->es, lpz.phi);
+
+
+    c  =  cos(lpz.phi);
+    if (fabs(c) < 1e-6) {
+        /* poleward of 89.99994 deg, we avoid division by zero   */
+        /* by computing the height as the cartesian z value      */
+        /* minus the geocentric radius of the Earth at the given */
+        /* latitude                                              */
+        double r = geocentric_radius (P->a, P->b, lpz.phi);
+        lpz.z = fabs (cart.z) - r;
+    }
+    else
+        lpz.z =  p / c  -  N;
+
+    return lpz;
+}
+
+
+
+/* In effect, 2 cartesian coordinates of a point on the ellipsoid. Rather pointless, but... */
+static XY cart_forward (LP lp, PJ *P) {
+    PJ_COORD point;
+    point.lp = lp;
+    point.lpz.z = 0;
+
+    point.xyz = cartesian (point.lpz, P);
+    return point.xy;
+}
+
+/* And the other way round. Still rather pointless, but... */
+static LP cart_reverse (XY xy, PJ *P) {
+    PJ_COORD point;
+    point.xy = xy;
+    point.xyz.z = 0;
+
+    point.lpz = geodetic (point.xyz, P);
+    return point.lp;
+}
+
+
+
+/*********************************************************************/
+PJ *CONVERSION(cart,1) {
+/*********************************************************************/
+    P->fwd3d  =  cartesian;
+    P->inv3d  =  geodetic;
+    P->fwd    =  cart_forward;
+    P->inv    =  cart_reverse;
+    P->left   =  PJ_IO_UNITS_ANGULAR;
+    P->right  =  PJ_IO_UNITS_CARTESIAN;
+    return P;
+}
diff --git a/src/conversions/PJ_geoc.cpp b/src/conversions/PJ_geoc.cpp
new file mode 100644
index 00000000..0455fada
--- /dev/null
+++ b/src/conversions/PJ_geoc.cpp
@@ -0,0 +1,59 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Conversion from geographic to geocentric latitude and back.
+ * Author:   Thomas Knudsen (2017)
+ *
+ ******************************************************************************
+ * Copyright (c) 2017, SDFE, http://www.sdfe.dk
+ * Copyright (c) 2017, Thomas Knudsen
+ *
+ * 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_LIB__
+
+#include <math.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+#include "projects.h"
+
+PROJ_HEAD(geoc, "Geocentric Latitude");
+
+/* Geographical to geocentric */
+static PJ_COORD forward(PJ_COORD coo, PJ *P) {
+    return pj_geocentric_latitude (P, PJ_FWD, coo);
+}
+
+/* Geocentric to geographical */
+static PJ_COORD inverse(PJ_COORD coo, PJ *P) {
+    return pj_geocentric_latitude (P, PJ_INV, coo);
+}
+
+
+static PJ *CONVERSION(geoc, 1) {
+    P->inv4d = inverse;
+    P->fwd4d = forward;
+
+    P->left   =  PJ_IO_UNITS_ANGULAR;
+    P->right  =  PJ_IO_UNITS_ANGULAR;
+
+    P->is_latlong = 1;
+    return P;
+}
diff --git a/src/conversions/PJ_unitconvert.cpp b/src/conversions/PJ_unitconvert.cpp
new file mode 100644
index 00000000..b25fd5d2
--- /dev/null
+++ b/src/conversions/PJ_unitconvert.cpp
@@ -0,0 +1,552 @@
+/***********************************************************************
+
+            Unit conversion pseudo-projection for use with
+                       transformation pipelines.
+
+                      Kristian Evers, 2017-05-16
+
+************************************************************************
+
+A pseudo-projection that can be used to convert units of input and
+output data. Primarily useful in pipelines.
+
+Unit conversion is performed by means of a pivot unit. The pivot unit
+for distance units are the meter and for time we use the modified julian
+date. A time unit conversion is performed like
+
+    Unit A -> Modified Julian date -> Unit B
+
+distance units are converted in the same manner, with meter being the
+central unit.
+
+The modified Julian date is chosen as the pivot unit since it has a
+fairly high precision, goes sufficiently long backwards in time, has no
+danger of hitting the upper limit in the near future and it is a fairly
+common time unit in astronomy and geodesy. Note that we are using the
+Julian date and not day. The difference being that the latter is defined
+as an integer and is thus limited to days in resolution. This approach
+has been extended wherever it makes sense, e.g. the GPS week unit also
+has a fractional part that makes it possible to determine the day, hour
+and minute of an observation.
+
+In- and output units are controlled with the parameters
+
+    +xy_in, +xy_out, +z_in, +z_out, +t_in and +t_out
+
+where xy denotes horizontal units, z vertical units and t time units.
+
+************************************************************************
+
+Kristian Evers, kreve@sdfe.dk, 2017-05-09
+Last update: 2017-05-16
+
+************************************************************************
+* Copyright (c) 2017, Kristian Evers / 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_LIB__
+
+#include <errno.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+
+#include "proj_internal.h"
+#include "proj_math.h"
+#include "projects.h"
+
+PROJ_HEAD(unitconvert, "Unit conversion");
+
+typedef double (*tconvert)(double);
+
+namespace { // anonymous namespace
+struct TIME_UNITS {
+    const char  *id;        /* units keyword */
+    tconvert     t_in;      /* unit -> mod. julian date function pointer */
+    tconvert     t_out;     /* mod. julian date > unit function pointer */
+    const char  *name;      /* comments */
+};
+} // anonymous namespace
+
+namespace { // anonymous namespace
+struct pj_opaque_unitconvert {
+    int     t_in_id;        /* time unit id for the time input unit   */
+    int     t_out_id;       /* time unit id for the time output unit  */
+    double  xy_factor;      /* unit conversion factor for horizontal components */
+    double  z_factor;       /* unit conversion factor for vertical components */
+};
+} // anonymous namespace
+
+
+/***********************************************************************/
+static int is_leap_year(long year) {
+/***********************************************************************/
+    return ((year % 4 == 0 && year % 100 != 0) || year % 400 ==0);
+}
+
+
+/***********************************************************************/
+static int days_in_year(long year) {
+/***********************************************************************/
+    return is_leap_year(year) ? 366 : 365;
+}
+
+/***********************************************************************/
+static unsigned int days_in_month(unsigned long year, unsigned long month) {
+/***********************************************************************/
+    const unsigned int month_table[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
+    unsigned int days;
+
+    if (month > 12) month = 12;
+    if (month == 0) month = 1;
+
+    days = month_table[month-1];
+    if (is_leap_year(year) && month == 2) days++;
+
+    return days;
+}
+
+
+/***********************************************************************/
+static int daynumber_in_year(unsigned long year, unsigned long month, unsigned long day) {
+/***********************************************************************/
+    unsigned int daynumber=0, i;
+
+    if (month > 12) month = 12;
+    if (month == 0) month = 1;
+    if (day > days_in_month(year, month)) day = days_in_month(year, month);
+
+    for (i = 1; i < month; i++)
+        daynumber += days_in_month(year, i);
+
+    daynumber += day;
+
+    return daynumber;
+
+}
+
+/***********************************************************************/
+static double mjd_to_mjd(double mjd) {
+/***********************************************************************
+    Modified julian date no-op function.
+
+    The Julian date is defined as (fractional) days since midnight
+    on 16th of November in 1858.
+************************************************************************/
+    return mjd;
+}
+
+
+/***********************************************************************/
+static double decimalyear_to_mjd(double decimalyear) {
+/***********************************************************************
+    Epoch of modified julian date is 1858-11-16 00:00
+************************************************************************/
+    long year;
+    double fractional_year;
+    double mjd;
+
+    if( decimalyear < -10000 || decimalyear > 10000 )
+        return 0;
+
+    year = lround(floor(decimalyear));
+    fractional_year = decimalyear - year;
+    mjd = (year - 1859)*365 + 14 + 31;
+    mjd += (double)fractional_year*(double)days_in_year(year);
+
+    /* take care of leap days */
+    year--;
+    for (; year > 1858; year--)
+        if (is_leap_year(year))
+            mjd++;
+
+    return mjd;
+}
+
+
+/***********************************************************************/
+static double mjd_to_decimalyear(double mjd) {
+/***********************************************************************
+    Epoch of modified julian date is 1858-11-16 00:00
+************************************************************************/
+    double decimalyear = mjd;
+    double mjd_iter = 14 + 31;
+    int year = 1859;
+
+    /* a smarter brain than mine could probably to do this more elegantly
+       - I'll just brute-force my way out of this... */
+    for (; mjd >= mjd_iter; year++) {
+        mjd_iter += days_in_year(year);
+    }
+    year--;
+    mjd_iter -= days_in_year(year);
+
+    decimalyear = year + (mjd-mjd_iter)/days_in_year(year);
+    return decimalyear;
+}
+
+
+/***********************************************************************/
+static double gps_week_to_mjd(double gps_week) {
+/***********************************************************************
+    GPS weeks are defined as the number of weeks since January the 6th
+    1980.
+
+    Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian
+    date is 44244.
+************************************************************************/
+    return 44244.0 + gps_week*7.0;
+}
+
+
+/***********************************************************************/
+static double mjd_to_gps_week(double mjd) {
+/***********************************************************************
+    GPS weeks are defined as the number of weeks since January the 6th
+    1980.
+
+    Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian
+    date is 44244.
+************************************************************************/
+    return (mjd - 44244.0) / 7.0;
+}
+
+
+/***********************************************************************/
+static double yyyymmdd_to_mjd(double yyyymmdd) {
+/************************************************************************
+    Date given in YYYY-MM-DD format.
+************************************************************************/
+
+    long year  = lround(floor( yyyymmdd / 10000 ));
+    long month = lround(floor((yyyymmdd - year*10000) / 100));
+    long day   = lround(floor( yyyymmdd - year*10000 - month*100));
+    double mjd = daynumber_in_year(year, month, day);
+
+    for (year -= 1; year > 1858; year--)
+        mjd += days_in_year(year);
+
+    return mjd + 13 + 31;
+}
+
+
+/***********************************************************************/
+static double mjd_to_yyyymmdd(double mjd) {
+/************************************************************************
+    Date given in YYYY-MM-DD format.
+************************************************************************/
+    double mjd_iter = 14 + 31;
+    int year = 1859, month=0, day=0;
+
+    for (; mjd >= mjd_iter; year++) {
+        mjd_iter += days_in_year(year);
+    }
+    year--;
+    mjd_iter -= days_in_year(year);
+
+    for (month=1; mjd_iter + days_in_month(year, month) <= mjd; month++)
+        mjd_iter += days_in_month(year, month);
+
+    day = (int)(mjd - mjd_iter + 1);
+
+    return year*10000.0 + month*100.0 + day;
+}
+
+static const struct TIME_UNITS time_units[] = {
+    {"mjd",         mjd_to_mjd,         mjd_to_mjd,         "Modified julian date"},
+    {"decimalyear", decimalyear_to_mjd, mjd_to_decimalyear, "Decimal year"},
+    {"gps_week",    gps_week_to_mjd,    mjd_to_gps_week,    "GPS Week"},
+    {"yyyymmdd",    yyyymmdd_to_mjd,    mjd_to_yyyymmdd,    "YYYYMMDD date"},
+    {nullptr,          nullptr,               nullptr,               nullptr}
+};
+
+
+/***********************************************************************/
+static XY forward_2d(LP lp, PJ *P) {
+/************************************************************************
+    Forward unit conversions in the plane
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD point = {{0,0,0,0}};
+    point.lp = lp;
+
+    point.xy.x *= Q->xy_factor;
+    point.xy.y *= Q->xy_factor;
+
+    return point.xy;
+}
+
+
+/***********************************************************************/
+static LP reverse_2d(XY xy, PJ *P) {
+/************************************************************************
+    Reverse unit conversions in the plane
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD point = {{0,0,0,0}};
+    point.xy = xy;
+
+    point.xy.x /= Q->xy_factor;
+    point.xy.y /= Q->xy_factor;
+
+    return point.lp;
+}
+
+
+/***********************************************************************/
+static XYZ forward_3d(LPZ lpz, PJ *P) {
+/************************************************************************
+    Forward unit conversions the vertical component
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD point = {{0,0,0,0}};
+    point.lpz = lpz;
+
+    /* take care of the horizontal components in the 2D function */
+    point.xy = forward_2d(point.lp, P);
+
+    point.xyz.z *= Q->z_factor;
+
+    return point.xyz;
+}
+
+/***********************************************************************/
+static LPZ reverse_3d(XYZ xyz, PJ *P) {
+/************************************************************************
+    Reverse unit conversions the vertical component
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD point = {{0,0,0,0}};
+    point.xyz = xyz;
+
+    /* take care of the horizontal components in the 2D function */
+    point.lp = reverse_2d(point.xy, P);
+
+    point.xyz.z /= Q->z_factor;
+
+    return point.lpz;
+}
+
+
+/***********************************************************************/
+static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+/************************************************************************
+    Forward conversion of time units
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD out = obs;
+
+    /* delegate unit conversion of physical dimensions to the 3D function */
+    out.xyz = forward_3d(obs.lpz, P);
+
+    if (Q->t_in_id >= 0)
+        out.xyzt.t = time_units[Q->t_in_id].t_in( obs.xyzt.t );
+    if (Q->t_out_id >= 0)
+        out.xyzt.t = time_units[Q->t_out_id].t_out( out.xyzt.t );
+
+    return out;
+}
+
+
+/***********************************************************************/
+static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+/************************************************************************
+    Reverse conversion of time units
+************************************************************************/
+    struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
+    PJ_COORD out = obs;
+
+    /* delegate unit conversion of physical dimensions to the 3D function */
+    out.lpz = reverse_3d(obs.xyz, P);
+
+    if (Q->t_out_id >= 0)
+        out.xyzt.t = time_units[Q->t_out_id].t_in( obs.xyzt.t );
+    if (Q->t_in_id >= 0)
+        out.xyzt.t = time_units[Q->t_in_id].t_out( out.xyzt.t );
+
+    return out;
+}
+
+/***********************************************************************/
+static double get_unit_conversion_factor(const char* name,
+                                         int* p_is_linear,
+                                         const char** p_normalized_name) {
+/***********************************************************************/
+    int i;
+    const char* s;
+    const PJ_UNITS *units;
+
+    units = proj_list_units();
+
+    /* Try first with linear units */
+    for (i = 0; (s = units[i].id) ; ++i) {
+        if ( strcmp(s, name) == 0 ) {
+            if( p_normalized_name ) {
+                *p_normalized_name = units[i].name;
+            }
+            if( p_is_linear ) {
+                *p_is_linear = 1;
+            }
+            return units[i].factor;
+        }
+    }
+
+    /* And then angular units */
+    units = proj_list_angular_units();
+    for (i = 0; (s = units[i].id) ; ++i) {
+        if ( strcmp(s, name) == 0 ) {
+            if( p_normalized_name ) {
+                *p_normalized_name = units[i].name;
+            }
+            if( p_is_linear ) {
+                *p_is_linear = 0;
+            }
+            return units[i].factor;
+        }
+    }
+    if( p_normalized_name ) {
+        *p_normalized_name = nullptr;
+    }
+    if( p_is_linear ) {
+        *p_is_linear = -1;
+    }
+    return 0.0;
+}
+
+/***********************************************************************/
+PJ *CONVERSION(unitconvert,0) {
+/***********************************************************************/
+    struct pj_opaque_unitconvert *Q = static_cast<struct pj_opaque_unitconvert*>(pj_calloc (1, sizeof (struct pj_opaque_unitconvert)));
+    const char *s, *name;
+    int i;
+    double f;
+    int xy_in_is_linear = -1; /* unknown */
+    int xy_out_is_linear = -1; /* unknown */
+    int z_in_is_linear = -1; /* unknown */
+    int z_out_is_linear = -1; /* unknown */
+
+    if (nullptr==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = (void *) Q;
+
+    P->fwd4d  = forward_4d;
+    P->inv4d  = reverse_4d;
+    P->fwd3d  = forward_3d;
+    P->inv3d  = reverse_3d;
+    P->fwd    = forward_2d;
+    P->inv    = reverse_2d;
+
+    P->left  = PJ_IO_UNITS_WHATEVER;
+    P->right = PJ_IO_UNITS_WHATEVER;
+
+    /* if no time input/output unit is specified we can skip them */
+    Q->t_in_id = -1;
+    Q->t_out_id = -1;
+
+    Q->xy_factor = 1.0;
+    Q->z_factor  = 1.0;
+
+    if ((name = pj_param (P->ctx, P->params, "sxy_in").s) != nullptr) {
+        const char* normalized_name = nullptr;
+        f = get_unit_conversion_factor(name, &xy_in_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "xy_in unit: %s", normalized_name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dxy_in").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->xy_factor *= f;
+    }
+
+    if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != nullptr) {
+        const char* normalized_name = nullptr;
+        f = get_unit_conversion_factor(name, &xy_out_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "xy_out unit: %s", normalized_name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dxy_out").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->xy_factor /= f;
+    }
+
+    if( xy_in_is_linear >= 0 && xy_out_is_linear >= 0 &&
+        xy_in_is_linear != xy_out_is_linear ) {
+        proj_log_debug(P, "inconsistent unit type between xy_in and xy_out");
+        return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT);
+    }
+
+    if ((name = pj_param (P->ctx, P->params, "sz_in").s) != nullptr) {
+        const char* normalized_name = nullptr;
+        f = get_unit_conversion_factor(name, &z_in_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "z_in unit: %s", normalized_name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dz_in").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->z_factor *= f;
+    }
+
+    if ((name = pj_param (P->ctx, P->params, "sz_out").s) != nullptr) {
+        const char* normalized_name = nullptr;
+        f = get_unit_conversion_factor(name, &z_out_is_linear, &normalized_name);
+        if (f != 0.0) {
+            proj_log_debug(P, "z_out unit: %s", normalized_name);
+        } else {
+            if ( (f = pj_param (P->ctx, P->params, "dz_out").f) == 0.0)
+                return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID);
+        }
+        if (f != 0.0)
+            Q->z_factor /= f;
+    }
+
+    if( z_in_is_linear >= 0 && z_out_is_linear >= 0 &&
+        z_in_is_linear != z_out_is_linear ) {
+        proj_log_debug(P, "inconsistent unit type between z_in and z_out");
+        return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT);
+    }
+
+    if ((name = pj_param (P->ctx, P->params, "st_in").s) != nullptr) {
+        for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i);
+
+        if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */
+
+        Q->t_in_id = i;
+        proj_log_debug(P, "t_in unit: %s", time_units[i].name);
+    }
+
+    s = nullptr;
+    if ((name = pj_param (P->ctx, P->params, "st_out").s) != nullptr) {
+        for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i);
+
+        if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */
+
+        Q->t_out_id = i;
+        proj_log_debug(P, "t_out unit: %s", time_units[i].name);
+    }
+
+    return P;
+}
diff --git a/src/conversions/pj_geocent.cpp b/src/conversions/pj_geocent.cpp
new file mode 100644
index 00000000..0e9d725e
--- /dev/null
+++ b/src/conversions/pj_geocent.cpp
@@ -0,0 +1,62 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Stub projection for geocentric.  The transformation isn't
+ *           really done here since this code is 2D.  The real transformation
+ *           is handled by pj_transform.c.
+ * Author:   Frank Warmerdam, warmerdam@pobox.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2002, Frank Warmerdam
+ *
+ * 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_LIB__
+#include "projects.h"
+
+PROJ_HEAD(geocent, "Geocentric")  "\n\t";
+
+static XY forward(LP lp, PJ *P) {
+    XY xy = {0.0,0.0};
+    (void) P;
+    xy.x = lp.lam;
+    xy.y = lp.phi;
+    return xy;
+}
+
+static LP inverse(XY xy, PJ *P) {
+    LP lp = {0.0,0.0};
+    (void) P;
+    lp.phi = xy.y;
+    lp.lam = xy.x;
+    return lp;
+}
+
+PJ *CONVERSION (geocent, 0) {
+    P->is_geocent = 1;
+    P->x0 = 0.0;
+    P->y0 = 0.0;
+    P->inv = inverse;
+    P->fwd = forward;
+    P->left = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_CARTESIAN;
+
+    return P;
+}
+