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

17 files changed, 6319 insertions, 0 deletions

diff --git a/src/apps/cct.cpp b/src/apps/cct.cpp
new file mode 100644
index 00000000..046257da
--- /dev/null
+++ b/src/apps/cct.cpp
@@ -0,0 +1,472 @@
+/***********************************************************************
+
+                 The cct 4D Transformation program
+
+************************************************************************
+
+cct is a 4D equivalent to the "proj" projection program.
+
+cct is an acronym meaning "Coordinate Conversion and Transformation".
+
+The acronym refers to definitions given in the OGC 08-015r2/ISO-19111
+standard "Geographical Information -- Spatial Referencing by Coordinates",
+which defines two different classes of coordinate operations:
+
+*Coordinate Conversions*, which are coordinate operations where input
+and output datum are identical (e.g. conversion from geographical to
+cartesian coordinates) and
+
+*Coordinate Transformations*, which are coordinate operations where
+input and output datums differ (e.g. change of reference frame).
+
+cct, however, also refers to Carl Christian Tscherning (1942--2014),
+professor of Geodesy at the University of Copenhagen, mentor and advisor
+for a generation of Danish geodesists, colleague and collaborator for
+two generations of global geodesists, Secretary General for the
+International Association of Geodesy, IAG (1995--2007), fellow of the
+American Geophysical Union (1991), recipient of the IAG Levallois Medal
+(2007), the European Geosciences Union Vening Meinesz Medal (2008), and
+of numerous other honours.
+
+cct, or Christian, as he was known to most of us, was recognized for his
+good mood, his sharp wit, his tireless work, and his great commitment to
+the development of geodesy - both through his scientific contributions,
+comprising more than 250 publications, and by his mentoring and teaching
+of the next generations of geodesists.
+
+As Christian was an avid Fortran programmer, and a keen Unix connoisseur,
+he would have enjoyed to know that his initials would be used to name a
+modest Unix style transformation filter, hinting at the tireless aspect
+of his personality, which was certainly one of the reasons he accomplished
+so much, and meant so much to so many people.
+
+Hence, in honour of cct (the geodesist) this is cct (the program).
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
+
+************************************************************************
+
+* Copyright (c) 2016, 2017 Thomas Knudsen
+* Copyright (c) 2017, 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.
+
+***********************************************************************/
+
+#include <ctype.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdarg.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+#include "proj_strtod.h"
+#include "projects.h"
+#include "optargpm.h"
+
+
+static void logger(void *data, int level, const char *msg);
+static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...);
+
+/* Prototypes from functions in this file */
+char *column (char *buf, int n);
+PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time);
+
+
+static const char usage[] = {
+    "--------------------------------------------------------------------------------\n"
+    "Usage: %s [-options]... [+operator_specs]... infile...\n"
+    "--------------------------------------------------------------------------------\n"
+    "Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    -c x,y,z,t        Specify input columns for (up to) 4 input parameters.\n"
+    "                      Defaults to 1,2,3,4\n"
+    "    -d n              Specify number of decimals in output.\n"
+    "    -I                Do the inverse transformation\n"
+    "    -o /path/to/file  Specify output file name\n"
+    "    -t value          Provide a fixed t value for all input data (e.g. -t 0)\n"
+    "    -z value          Provide a fixed z value for all input data (e.g. -z 0)\n"
+    "    -s n              Skip n first lines of a infile\n"
+    "    -v                Verbose: Provide non-essential informational output.\n"
+    "                      Repeat -v for more verbosity (e.g. -vv)\n"
+    "--------------------------------------------------------------------------------\n"
+    "Long Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    --output          Alias for -o\n"
+    "    --columns         Alias for -c\n"
+    "    --decimals        Alias for -d\n"
+    "    --height          Alias for -z\n"
+    "    --time            Alias for -t\n"
+    "    --verbose         Alias for -v\n"
+    "    --inverse         Alias for -I\n"
+    "    --skip-lines      Alias for -s\n"
+    "    --help            Alias for -h\n"
+    "    --version         Print version number\n"
+    "--------------------------------------------------------------------------------\n"
+    "Operator Specs:\n"
+    "--------------------------------------------------------------------------------\n"
+    "The operator specs describe the action to be performed by cct, e.g:\n"
+    "\n"
+    "        +proj=utm  +ellps=GRS80  +zone=32\n"
+    "\n"
+    "instructs cct to convert input data to Universal Transverse Mercator, zone 32\n"
+    "coordinates, based on the GRS80 ellipsoid.\n"
+    "\n"
+    "Hence, the command\n"
+    "\n"
+    "        echo 12 55 | cct -z0 -t0 +proj=utm +zone=32 +ellps=GRS80\n"
+    "\n"
+    "Should give results comparable to the classic proj command\n"
+    "\n"
+    "        echo 12 55 | proj +proj=utm +zone=32 +ellps=GRS80\n"
+    "--------------------------------------------------------------------------------\n"
+    "Examples:\n"
+    "--------------------------------------------------------------------------------\n"
+    "1. convert geographical input to UTM zone 32 on the GRS80 ellipsoid:\n"
+    "    cct +proj=utm +ellps=GRS80 +zone=32\n"
+    "2. roundtrip accuracy check for the case above:\n"
+    "    cct +proj=pipeline +proj=utm +ellps=GRS80 +zone=32 +step +step +inv\n"
+    "3. as (1) but specify input columns for longitude, latitude, height and time:\n"
+    "    cct -c 5,2,1,4  +proj=utm +ellps=GRS80 +zone=32\n"
+    "4. as (1) but specify fixed height and time, hence needing only 2 cols in input:\n"
+    "    cct -t 0 -z 0  +proj=utm  +ellps=GRS80  +zone=32\n"
+    "--------------------------------------------------------------------------------\n"
+};
+
+
+static void logger(void *data, int level, const char *msg) {
+    FILE *stream;
+    int log_tell = proj_log_level(PJ_DEFAULT_CTX, PJ_LOG_TELL);
+
+    stream  = (FILE *) data;
+
+    /* if we use PJ_LOG_NONE we always want to print stuff to stream */
+    if (level == PJ_LOG_NONE) {
+        fprintf(stream, "%s", msg);
+        return;
+    }
+
+    /* should always print to stderr if level == PJ_LOG_ERROR */
+    if (level == PJ_LOG_ERROR) {
+        fprintf(stderr, "%s", msg);
+        return;
+    }
+
+    /* otherwise only print if log level set by user is high enough */
+    if (level <= log_tell)
+        fprintf(stream, "%s", msg);
+}
+
+FILE *fout;
+
+static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...) {
+
+    va_list args;
+    char *msg_buf;
+
+    va_start( args, fmt );
+
+    msg_buf = (char *) malloc(100000);
+    if( msg_buf == nullptr ) {
+        va_end( args );
+        return;
+    }
+
+    vsprintf( msg_buf, fmt, args );
+
+    logger((void *) fout, log_level, msg_buf);
+
+    va_end( args );
+    free( msg_buf );
+}
+
+
+int main(int argc, char **argv) {
+    PJ *P;
+    PJ_COORD point;
+    PJ_PROJ_INFO info;
+    OPTARGS *o;
+    char blank_comment[] = "";
+    char whitespace[] = " ";
+    char *comment;
+    char *comment_delimiter;
+    char *buf;
+    int i, nfields = 4, skip_lines = 0, verbose;
+    double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL;
+    int decimals_angles = 10;
+    int decimals_distances = 4;
+    int columns_xyzt[] = {1, 2, 3, 4};
+    const char *longflags[]  = {"v=verbose", "h=help", "I=inverse", "version", nullptr};
+    const char *longkeys[]   = {
+        "o=output",
+        "c=columns",
+        "d=decimals",
+        "z=height",
+        "t=time",
+        "s=skip-lines",
+        nullptr};
+
+    fout = stdout;
+
+    o = opt_parse (argc, argv, "hvI", "cdozts", longflags, longkeys);
+    if (nullptr==o)
+        return 0;
+
+    if (opt_given (o, "h") || argc==1) {
+        printf (usage, o->progname);
+        return 0;
+    }
+
+    PJ_DIRECTION direction = opt_given (o, "I")? PJ_INV: PJ_FWD;
+
+    verbose   = MIN(opt_given (o, "v"), 3); /* log level can't be larger than 3 */
+    proj_log_level (PJ_DEFAULT_CTX, static_cast<PJ_LOG_LEVEL>(verbose));
+    proj_log_func  (PJ_DEFAULT_CTX, (void *) fout, logger);
+
+    if (opt_given (o, "version")) {
+        print (PJ_LOG_NONE, "%s: %s\n", o->progname, pj_get_release ());
+        return 0;
+    }
+
+    if (opt_given (o, "o"))
+        fout = fopen (opt_arg (o, "output"), "wt");
+    if (nullptr==fout) {
+        print (PJ_LOG_ERROR, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
+        free (o);
+        return 1;
+    }
+
+    print (PJ_LOG_TRACE, "%s: Running in very verbose mode\n", o->progname);
+
+    if (opt_given (o, "z")) {
+        fixed_z = proj_atof (opt_arg (o, "z"));
+        nfields--;
+    }
+
+    if (opt_given (o, "t")) {
+        fixed_time = proj_atof (opt_arg (o, "t"));
+        nfields--;
+    }
+
+     if (opt_given (o, "d")) {
+        int dec = atoi (opt_arg (o, "d"));
+        decimals_angles = dec;
+        decimals_distances = dec;
+    }
+
+    if (opt_given (o, "s")) {
+        skip_lines = atoi (opt_arg(o, "s"));
+    }
+
+    if (opt_given (o, "c")) {
+        int ncols;
+        /* reset column numbers to ease comment output later on */
+        for (i=0; i<4; i++)
+            columns_xyzt[i] = 0;
+
+        /* cppcheck-suppress invalidscanf */
+        ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3);
+        if (ncols != nfields) {
+            print (PJ_LOG_ERROR, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c"));
+            free (o);
+            if (stdout != fout)
+                fclose (fout);
+            return 1;
+        }
+    }
+
+    /* Setup transformation */
+    P = proj_create_argv (nullptr, o->pargc, o->pargv);
+    if ((nullptr==P) || (0==o->pargc)) {
+        print (PJ_LOG_ERROR, "%s: Bad transformation arguments - (%s)\n    '%s -h' for help\n",
+                 o->progname, pj_strerrno (proj_errno(P)), o->progname);
+        free (o);
+        if (stdout != fout)
+            fclose (fout);
+        return 1;
+    }
+
+    info = proj_pj_info (P);
+    print (PJ_LOG_TRACE, "Final: %s argc=%d pargc=%d\n", info.definition, argc, o->pargc);
+
+    if (direction== PJ_INV) {
+        /* fail if an inverse operation is not available */
+        if (!info.has_inverse) {
+            print (PJ_LOG_ERROR, "Inverse operation not available\n");
+            if (stdout != fout)
+                fclose (fout);
+            return 1;
+        }
+        /* We have no API call for inverting an operation, so we brute force it. */
+        P->inverted = !(P->inverted);
+    }
+    direction = PJ_FWD;
+
+    /* Allocate input buffer */
+    buf = static_cast<char*>(calloc (1, 10000));
+    if (nullptr==buf) {
+        print (PJ_LOG_ERROR, "%s: Out of memory\n", o->progname);
+        pj_free (P);
+        free (o);
+        if (stdout != fout)
+            fclose (fout);
+        return 1;
+    }
+
+
+    /* Loop over all records of all input files */
+    while (opt_input_loop (o, optargs_file_format_text)) {
+        int err;
+        void *ret = fgets (buf, 10000, o->input);
+        char *c = column (buf, 1);
+        opt_eof_handler (o);
+        if (nullptr==ret) {
+            print (PJ_LOG_ERROR, "Read error in record %d\n", (int) o->record_index);
+            continue;
+        }
+        point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time);
+        if (skip_lines > 0) {
+            skip_lines--;
+            continue;
+        }
+
+        /* if it's a comment or blank line, we reflect it */
+        if (c && ((*c=='\0') || (*c=='#'))) {
+            fprintf (fout, "%s", buf);
+            continue;
+        }
+
+        if (HUGE_VAL==point.xyzt.x) {
+            /* otherwise, it must be a syntax error */
+            print (PJ_LOG_NONE, "# Record %d UNREADABLE: %s", (int) o->record_index, buf);
+            print (PJ_LOG_ERROR, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
+            continue;
+        }
+
+        if (proj_angular_input (P, direction)) {
+            point.lpzt.lam = proj_torad (point.lpzt.lam);
+            point.lpzt.phi = proj_torad (point.lpzt.phi);
+        }
+        err = proj_errno_reset (P);
+        point = proj_trans (P, direction, point);
+
+        if (HUGE_VAL==point.xyzt.x) {
+            /* transformation error */
+            print (PJ_LOG_NONE, "# Record %d TRANSFORMATION ERROR: %s (%s)",
+                   (int) o->record_index, buf, pj_strerrno (proj_errno(P)));
+            proj_errno_restore (P, err);
+            continue;
+        }
+        proj_errno_restore (P, err);
+
+        /* handle comment string */
+        comment = column(buf, nfields+1);
+        if (opt_given(o, "c")) {
+            /* what number is the last coordinate column in the input data? */
+            int colmax = 0;
+            for (i=0; i<4; i++)
+                colmax = MAX(colmax, columns_xyzt[i]);
+            comment = column(buf, colmax+1);
+        }
+        comment_delimiter = (comment && *comment) ? whitespace : blank_comment;
+
+        /* Time to print the result */
+        if (proj_angular_output (P, direction)) {
+            point.lpzt.lam = proj_todeg (point.lpzt.lam);
+            point.lpzt.phi = proj_todeg (point.lpzt.phi);
+            print (PJ_LOG_NONE, "%14.*f  %14.*f  %12.*f  %12.4f%s%s\n",
+                   decimals_angles, point.xyzt.x,
+                   decimals_angles, point.xyzt.y,
+                   decimals_distances, point.xyzt.z,
+                   point.xyzt.t, comment_delimiter, comment
+            );
+        }
+        else
+            print (PJ_LOG_NONE, "%13.*f  %13.*f  %12.*f  %12.4f%s%s\n",
+                   decimals_distances, point.xyzt.x,
+                   decimals_distances, point.xyzt.y,
+                   decimals_distances, point.xyzt.z,
+                   point.xyzt.t, comment_delimiter, comment
+            );
+    }
+
+    if (stdout != fout)
+        fclose (fout);
+    free (o);
+    free (buf);
+    return 0;
+}
+
+
+
+
+
+/* return a pointer to the n'th column of buf */
+char *column (char *buf, int n) {
+    int i;
+    if (n <= 0)
+        return buf;
+    for (i = 0;  i < n;  i++) {
+        while (isspace(*buf))
+            buf++;
+        if (i == n - 1)
+            break;
+        while ((0 != *buf) && !isspace(*buf))
+            buf++;
+    }
+    return buf;
+}
+
+/* column to double */
+static double cold (char *args, int col) {
+    char *endp;
+    char *target;
+    double d;
+    target = column (args, col);
+    d = proj_strtod (target, &endp);
+    if (endp==target)
+        return HUGE_VAL;
+    return d;
+}
+
+PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time) {
+    PJ_COORD err = proj_coord (HUGE_VAL, HUGE_VAL, HUGE_VAL, HUGE_VAL);
+    PJ_COORD result = err;
+    int prev_errno = errno;
+    errno = 0;
+
+    result.xyzt.z = fixed_height;
+    result.xyzt.t = fixed_time;
+    result.xyzt.x = cold (buf, columns[0]);
+    result.xyzt.y = cold (buf, columns[1]);
+    if (result.xyzt.z==HUGE_VAL)
+        result.xyzt.z = cold (buf, columns[2]);
+    if (result.xyzt.t==HUGE_VAL)
+        result.xyzt.t = cold (buf, columns[3]);
+
+    if (0!=errno)
+        return err;
+
+    errno = prev_errno;
+    return result;
+}
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
new file mode 100644
index 00000000..f63bedcc
--- /dev/null
+++ b/src/apps/cs2cs.cpp
@@ -0,0 +1,640 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Mainline program sort of like ``proj'' for converting between
+ *           two coordinate systems.
+ * Author:   Frank Warmerdam, warmerda@home.com
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, 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 FROM_PROJ_CPP
+
+#include <ctype.h>
+#include <locale.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <cassert>
+#include <string>
+
+#include <proj/internal/internal.hpp>
+
+// PROJ include order is sensitive
+// clang-format off
+#include "proj.h"
+#include "projects.h"
+#include "emess.h"
+// clang-format on
+
+#define MAX_LINE 1000
+
+static PJ *transformation = nullptr;
+
+static bool srcIsGeog = false;
+static double srcToRadians = 0.0;
+
+static bool destIsGeog = false;
+static double destToRadians = 0.0;
+static bool destIsLatLong = false;
+
+static int reversein = 0, /* != 0 reverse input arguments */
+    reverseout = 0,       /* != 0 reverse output arguments */
+    echoin = 0,           /* echo input data to output line */
+    tag = '#';            /* beginning of line tag character */
+
+static const char *oform =
+    nullptr;                  /* output format for x-y or decimal degrees */
+static char oform_buffer[16]; /* buffer for oform when using -d */
+static const char *oterr = "*\t*"; /* output line for unprojectable input */
+static const char *usage =
+    "%s\nusage: %s [ -dDeEfIlrstvwW [args] ] [ +opts[=arg] ]\n"
+    "                   [+to [+opts[=arg] [ files ]\n";
+
+static double (*informat)(const char *,
+                          char **); /* input data deformatter function */
+
+/************************************************************************/
+/*                              process()                               */
+/*                                                                      */
+/*      File processing function.                                       */
+/************************************************************************/
+static void process(FILE *fid)
+
+{
+    char line[MAX_LINE + 3], *s, pline[40];
+    projUV data;
+
+    for (;;) {
+        double z;
+
+        ++emess_dat.File_line;
+        if (!(s = fgets(line, MAX_LINE, fid)))
+            break;
+        if (!strchr(s, '\n')) { /* overlong line */
+            int c;
+            (void)strcat(s, "\n");
+            /* gobble up to newline */
+            while ((c = fgetc(fid)) != EOF && c != '\n')
+                ;
+        }
+        if (*s == tag) {
+            fputs(line, stdout);
+            continue;
+        }
+
+        if (reversein) {
+            data.v = (*informat)(s, &s);
+            data.u = (*informat)(s, &s);
+        } else {
+            data.u = (*informat)(s, &s);
+            data.v = (*informat)(s, &s);
+        }
+
+        z = strtod(s, &s);
+
+        if (data.v == HUGE_VAL)
+            data.u = HUGE_VAL;
+
+        if (!*s && (s > line))
+            --s; /* assumed we gobbled \n */
+
+        if (echoin) {
+            char t;
+            t = *s;
+            *s = '\0';
+            (void)fputs(line, stdout);
+            *s = t;
+            putchar('\t');
+        }
+
+        if (data.u != HUGE_VAL) {
+
+            if (srcIsGeog) {
+                /* dmstor gives values to radians. Convert now to the SRS unit
+                 */
+                data.u /= srcToRadians;
+                data.v /= srcToRadians;
+            }
+
+            PJ_COORD coord;
+            coord.xyzt.x = data.u;
+            coord.xyzt.y = data.v;
+            coord.xyzt.z = z;
+            coord.xyzt.t = HUGE_VAL;
+            coord = proj_trans(transformation, PJ_FWD, coord);
+            data.u = coord.xyz.x;
+            data.v = coord.xyz.y;
+            z = coord.xyz.z;
+        }
+
+        if (data.u == HUGE_VAL) /* error output */
+            fputs(oterr, stdout);
+
+        else if (destIsGeog && !oform) { /*ascii DMS output */
+
+            // rtodms() expect radians: convert from the output SRS unit
+            data.u *= destToRadians;
+            data.v *= destToRadians;
+
+            if (destIsLatLong) {
+                if (reverseout) {
+                    fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+                    putchar('\t');
+                    fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+                } else {
+                    fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+                    putchar('\t');
+                    fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+                }
+            } else if (reverseout) {
+                fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
+                putchar('\t');
+                fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
+            } else {
+                fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
+                putchar('\t');
+                fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
+            }
+
+        } else { /* x-y or decimal degree ascii output */
+            if (destIsGeog) {
+                data.v *= destToRadians * RAD_TO_DEG;
+                data.u *= destToRadians * RAD_TO_DEG;
+            }
+            if (reverseout) {
+                printf(oform, data.v);
+                putchar('\t');
+                printf(oform, data.u);
+            } else {
+                printf(oform, data.u);
+                putchar('\t');
+                printf(oform, data.v);
+            }
+        }
+
+        putchar(' ');
+        if (oform != nullptr)
+            printf(oform, z);
+        else
+            printf("%.3f", z);
+        if (s)
+            printf("%s", s);
+        else
+            printf("\n");
+    }
+}
+
+/************************************************************************/
+/*                          instanciate_crs()                           */
+/************************************************************************/
+
+static PJ_OBJ *instanciate_crs(const std::string &definition,
+                               const char *const *optionsImportCRS,
+                               bool &isGeog, double &toRadians,
+                               bool &isLatFirst) {
+    PJ_OBJ *crs = proj_obj_create_from_user_input(nullptr, definition.c_str(),
+                                                  optionsImportCRS);
+    if (!crs) {
+        return nullptr;
+    }
+
+    isGeog = false;
+    toRadians = 0.0;
+    isLatFirst = false;
+
+    auto type = proj_obj_get_type(crs);
+    if (type == PJ_OBJ_TYPE_BOUND_CRS) {
+        auto base = proj_obj_get_source_crs(nullptr, crs);
+        proj_obj_destroy(crs);
+        crs = base;
+        type = proj_obj_get_type(crs);
+    }
+    if (type == PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS ||
+        type == PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+        auto cs = proj_obj_crs_get_coordinate_system(nullptr, crs);
+        assert(cs);
+
+        isGeog = true;
+        const char *axisName = "";
+        proj_obj_cs_get_axis_info(nullptr, cs, 0,
+                                  &axisName, // name,
+                                  nullptr,   // abbrev
+                                  nullptr,   // direction
+                                  &toRadians,
+                                  nullptr, // unit name
+                                  nullptr, // unit authority
+                                  nullptr  // unit code
+                                  );
+        isLatFirst =
+            NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
+            std::string::npos;
+
+        proj_obj_destroy(cs);
+    }
+
+    return crs;
+}
+
+/************************************************************************/
+/*               get_geog_crs_proj_string_from_proj_crs()               */
+/************************************************************************/
+
+static std::string get_geog_crs_proj_string_from_proj_crs(PJ_OBJ *src,
+                                                          double &toRadians,
+                                                          bool &isLatFirst) {
+    auto srcType = proj_obj_get_type(src);
+    if (srcType == PJ_OBJ_TYPE_BOUND_CRS) {
+        auto base = proj_obj_get_source_crs(nullptr, src);
+        assert(base);
+        proj_obj_destroy(src);
+        src = base;
+        srcType = proj_obj_get_type(src);
+    }
+    if (srcType != PJ_OBJ_TYPE_PROJECTED_CRS) {
+        return std::string();
+    }
+
+    auto base = proj_obj_get_source_crs(nullptr, src);
+    assert(base);
+    auto baseType = proj_obj_get_type(base);
+    if (baseType != PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS &&
+        baseType != PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+        proj_obj_destroy(base);
+        return std::string();
+    }
+
+    auto cs = proj_obj_crs_get_coordinate_system(nullptr, base);
+    assert(cs);
+
+    const char *axisName = "";
+    proj_obj_cs_get_axis_info(nullptr, cs, 0,
+                              &axisName, // name,
+                              nullptr,   // abbrev
+                              nullptr,   // direction
+                              &toRadians,
+                              nullptr, // unit name
+                              nullptr, // unit authority
+                              nullptr  // unit code
+                              );
+    isLatFirst = NS_PROJ::internal::ci_find(std::string(axisName),
+                                            "latitude") != std::string::npos;
+
+    proj_obj_destroy(cs);
+
+    auto retCStr = proj_obj_as_proj_string(nullptr, base, PJ_PROJ_5, nullptr);
+    std::string ret(retCStr ? retCStr : "");
+    proj_obj_destroy(base);
+    return ret;
+}
+
+/************************************************************************/
+/*                                main()                                */
+/************************************************************************/
+
+int main(int argc, char **argv) {
+    char *arg;
+    char **eargv = argv;
+    std::string fromStr;
+    std::string toStr;
+    FILE *fid;
+    int eargc = 0, mon = 0;
+    int have_to_flag = 0, inverse = 0;
+    int use_env_locale = 0;
+
+    /* This is just to check that pj_init() is locale-safe */
+    /* Used by nad/testvarious */
+    if (getenv("PROJ_USE_ENV_LOCALE") != nullptr)
+        use_env_locale = 1;
+
+    /* Enable compatibility mode for init=epsg:XXXX by default */
+    if (getenv("PROJ_USE_PROJ4_INIT_RULES") == nullptr) {
+        proj_context_use_proj4_init_rules(nullptr, true);
+    }
+
+    if ((emess_dat.Prog_name = strrchr(*argv, DIR_CHAR)) != nullptr)
+        ++emess_dat.Prog_name;
+    else
+        emess_dat.Prog_name = *argv;
+    inverse = !strncmp(emess_dat.Prog_name, "inv", 3);
+    if (argc <= 1) {
+        (void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name);
+        exit(0);
+    }
+
+    // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST>" syntax
+    int countNonOptionArg = 0;
+    for (int i = 1; i < argc; i++) {
+        if (argv[i][0] == '-') {
+            if (argv[i][1] == '\0') {
+                countNonOptionArg++;
+            }
+        } else {
+            if (strcmp(argv[i], "+to") == 0) {
+                countNonOptionArg = -1;
+                break;
+            }
+            countNonOptionArg++;
+        }
+    }
+    const bool isSrcDestSyntax = (countNonOptionArg == 2);
+
+    /* process run line arguments */
+    while (--argc > 0) { /* collect run line arguments */
+        if (**++argv == '-') {
+            for (arg = *argv;;) {
+                switch (*++arg) {
+                case '\0': /* position of "stdin" */
+                    if (arg[-1] == '-')
+                        eargv[eargc++] = const_cast<char *>("-");
+                    break;
+                case 'v': /* monitor dump of initialization */
+                    mon = 1;
+                    continue;
+                case 'I': /* alt. method to spec inverse */
+                    inverse = 1;
+                    continue;
+                case 'E': /* echo ascii input to ascii output */
+                    echoin = 1;
+                    continue;
+                case 't': /* set col. one char */
+                    if (arg[1])
+                        tag = *++arg;
+                    else
+                        emess(1, "missing -t col. 1 tag");
+                    continue;
+                case 'l': /* list projections, ellipses or units */
+                    if (!arg[1] || arg[1] == 'p' || arg[1] == 'P') {
+                        /* list projections */
+                        const struct PJ_LIST *lp;
+                        int do_long = arg[1] == 'P', c;
+                        const char *str;
+
+                        for (lp = proj_list_operations(); lp->id; ++lp) {
+                            (void)printf("%s : ", lp->id);
+                            if (do_long) /* possibly multiline description */
+                                (void)puts(*lp->descr);
+                            else { /* first line, only */
+                                str = *lp->descr;
+                                while ((c = *str++) && c != '\n')
+                                    putchar(c);
+                                putchar('\n');
+                            }
+                        }
+                    } else if (arg[1] == '=') { /* list projection 'descr' */
+                        const struct PJ_LIST *lp;
+
+                        arg += 2;
+                        for (lp = proj_list_operations(); lp->id; ++lp)
+                            if (!strcmp(lp->id, arg)) {
+                                (void)printf("%9s : %s\n", lp->id, *lp->descr);
+                                break;
+                            }
+                    } else if (arg[1] == 'e') { /* list ellipses */
+                        const struct PJ_ELLPS *le;
+
+                        for (le = proj_list_ellps(); le->id; ++le)
+                            (void)printf("%9s %-16s %-16s %s\n", le->id,
+                                         le->major, le->ell, le->name);
+                    } else if (arg[1] == 'u') { /* list units */
+                        const struct PJ_UNITS *lu;
+
+                        for (lu = proj_list_units(); lu->id; ++lu)
+                            (void)printf("%12s %-20s %s\n", lu->id,
+                                         lu->to_meter, lu->name);
+                    } else if (arg[1] == 'd') { /* list datums */
+                        const struct PJ_DATUMS *ld;
+
+                        printf("__datum_id__ __ellipse___ "
+                               "__definition/"
+                               "comments______________________________\n");
+                        for (ld = pj_get_datums_ref(); ld->id; ++ld) {
+                            printf("%12s %-12s %-30s\n", ld->id, ld->ellipse_id,
+                                   ld->defn);
+                            if (ld->comments != nullptr &&
+                                strlen(ld->comments) > 0)
+                                printf("%25s %s\n", " ", ld->comments);
+                        }
+                    } else if (arg[1] == 'm') { /* list prime meridians */
+                        const struct PJ_PRIME_MERIDIANS *lpm;
+
+                        for (lpm = proj_list_prime_meridians(); lpm->id; ++lpm)
+                            (void)printf("%12s %-30s\n", lpm->id, lpm->defn);
+                    } else
+                        emess(1, "invalid list option: l%c", arg[1]);
+                    exit(0);
+                    /* cppcheck-suppress duplicateBreak */
+                    continue; /* artificial */
+                case 'e':     /* error line alternative */
+                    if (--argc <= 0)
+                    noargument:
+                    emess(1, "missing argument for -%c", *arg);
+                    oterr = *++argv;
+                    continue;
+                case 'W': /* specify seconds precision */
+                case 'w': /* -W for constant field width */
+                {
+                    char c = arg[1];
+                    if (c != 0 && isdigit(c)) {
+                        set_rtodms(c - '0', *arg == 'W');
+                        ++arg;
+                    } else
+                        emess(1, "-W argument missing or non-digit");
+                    continue;
+                }
+                case 'f': /* alternate output format degrees or xy */
+                    if (--argc <= 0)
+                        goto noargument;
+                    oform = *++argv;
+                    continue;
+                case 'r': /* reverse input */
+                    reversein = 1;
+                    continue;
+                case 's': /* reverse output */
+                    reverseout = 1;
+                    continue;
+                case 'D': /* set debug level */
+                    if (--argc <= 0)
+                        goto noargument;
+                    pj_ctx_set_debug(pj_get_default_ctx(), atoi(*++argv));
+                    continue;
+                case 'd':
+                    if (--argc <= 0)
+                        goto noargument;
+                    sprintf(oform_buffer, "%%.%df", atoi(*++argv));
+                    oform = oform_buffer;
+                    break;
+                default:
+                    emess(1, "invalid option: -%c", *arg);
+                    break;
+                }
+                break;
+            }
+        } else if (isSrcDestSyntax) {
+            if (fromStr.empty())
+                fromStr = *argv;
+            else
+                toStr = *argv;
+        } else if (strcmp(*argv, "+to") == 0) {
+            have_to_flag = 1;
+
+        } else if (**argv == '+') { /* + argument */
+            if (have_to_flag) {
+                if (!toStr.empty())
+                    toStr += ' ';
+                toStr += *argv;
+            } else {
+                if (!fromStr.empty())
+                    fromStr += ' ';
+                fromStr += *argv;
+            }
+        } else if (!have_to_flag) {
+            fromStr = *argv;
+        } else if (toStr.empty()) {
+            toStr = *argv;
+        } else /* assumed to be input file name(s) */
+            eargv[eargc++] = *argv;
+    }
+    if (eargc == 0) /* if no specific files force sysin */
+        eargv[eargc++] = const_cast<char *>("-");
+
+    /*
+     * If the user has requested inverse, then just reverse the
+     * coordinate systems.
+     */
+    if (inverse) {
+        std::swap(fromStr, toStr);
+    }
+
+    if (use_env_locale) {
+        /* Set locale from environment */
+        setlocale(LC_ALL, "");
+    }
+
+    if (fromStr.empty() && toStr.empty()) {
+        emess(3, "missing source and target coordinate systems");
+    }
+
+    const char *const optionsProj4Mode[] = {"USE_PROJ4_INIT_RULES=YES",
+                                            nullptr};
+    const char *const *optionsImportCRS =
+        proj_context_get_use_proj4_init_rules(nullptr, TRUE) ? optionsProj4Mode
+                                                             : nullptr;
+
+    PJ_OBJ *src = nullptr;
+    if (!fromStr.empty()) {
+        bool ignored;
+        src = instanciate_crs(fromStr, optionsImportCRS, srcIsGeog,
+                              srcToRadians, ignored);
+        if (!src) {
+            emess(3, "cannot instanciate source coordinate system");
+        }
+    }
+
+    PJ_OBJ *dst = nullptr;
+    if (!toStr.empty()) {
+        dst = instanciate_crs(toStr, optionsImportCRS, destIsGeog,
+                              destToRadians, destIsLatLong);
+        if (!dst) {
+            emess(3, "cannot instanciate target coordinate system");
+        }
+    }
+
+    if (toStr.empty()) {
+        assert(src);
+        toStr = get_geog_crs_proj_string_from_proj_crs(src, destToRadians,
+                                                       destIsLatLong);
+        if (toStr.empty()) {
+            emess(3,
+                  "missing target CRS and source CRS is not a projected CRS");
+        }
+        destIsGeog = true;
+    } else if (fromStr.empty()) {
+        assert(dst);
+        bool ignored;
+        fromStr =
+            get_geog_crs_proj_string_from_proj_crs(dst, srcToRadians, ignored);
+        if (fromStr.empty()) {
+            emess(3,
+                  "missing source CRS and target CRS is not a projected CRS");
+        }
+        srcIsGeog = true;
+    }
+
+    proj_obj_destroy(src);
+    proj_obj_destroy(dst);
+
+    transformation = proj_create_crs_to_crs(nullptr, fromStr.c_str(),
+                                            toStr.c_str(), nullptr);
+    if (!transformation) {
+        emess(3, "cannot initialize transformation\ncause: %s",
+              pj_strerrno(pj_errno));
+    }
+
+    if (use_env_locale) {
+        /* Restore C locale to avoid issues in parsing/outputting numbers*/
+        setlocale(LC_ALL, "C");
+    }
+
+    if (mon) {
+        printf("%c ---- From Coordinate System ----\n", tag);
+        printf("%s\n", fromStr.c_str());
+        printf("%c ---- To Coordinate System ----\n", tag);
+        printf("%s\n", toStr.c_str());
+    }
+
+    /* set input formatting control */
+    if (!srcIsGeog)
+        informat = strtod;
+    else {
+        informat = dmstor;
+    }
+
+    if (!destIsGeog && !oform)
+        oform = "%.2f";
+
+    /* process input file list */
+    for (; eargc--; ++eargv) {
+        if (**eargv == '-') {
+            fid = stdin;
+            emess_dat.File_name = const_cast<char *>("<stdin>");
+
+        } else {
+            if ((fid = fopen(*eargv, "rt")) == nullptr) {
+                emess(-2, *eargv, "input file");
+                continue;
+            }
+            emess_dat.File_name = *eargv;
+        }
+        emess_dat.File_line = 0;
+        process(fid);
+        fclose(fid);
+        emess_dat.File_name = nullptr;
+    }
+
+    proj_destroy(transformation);
+
+    pj_deallocate_grids();
+
+    exit(0); /* normal completion */
+}
diff --git a/src/apps/emess.cpp b/src/apps/emess.cpp
new file mode 100644
index 00000000..144e9e23
--- /dev/null
+++ b/src/apps/emess.cpp
@@ -0,0 +1,68 @@
+/* Error message processing */
+
+#ifdef _MSC_VER
+#  ifndef _CRT_SECURE_NO_DEPRECATE
+#    define _CRT_SECURE_NO_DEPRECATE
+#  endif
+#  ifndef _CRT_NONSTDC_NO_DEPRECATE
+#    define _CRT_NONSTDC_NO_DEPRECATE
+#  endif
+#endif
+
+#ifndef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
+#define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
+#endif
+
+#include <errno.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "proj_api.h"
+#define EMESS_ROUTINE
+#include "emess.h"
+
+	void
+emess(int code, const char *fmt, ...) {
+	va_list args;
+
+	va_start(args, fmt);
+	/* prefix program name, if given */
+	if (fmt != nullptr)
+		(void)fprintf(stderr,"%s\n<%s>: ",pj_get_release(),
+                              emess_dat.Prog_name);
+	/* print file name and line, if given */
+	if (emess_dat.File_name != nullptr && *emess_dat.File_name) {
+		(void)fprintf(stderr,"while processing file: %s", emess_dat.File_name);
+		if (emess_dat.File_line > 0)
+			(void)fprintf(stderr,", line %d\n", emess_dat.File_line);
+		else
+			(void)fputc('\n', stderr);
+	} else
+		putc('\n', stderr);
+	/* if |code|==2, print errno code data */
+	if (code == 2 || code == -2)
+        {
+                int my_errno = errno;
+#ifdef HAVE_STRERROR
+                const char* my_strerror = strerror(my_errno);
+#endif
+#ifndef HAVE_STRERROR
+                const char* my_strerror = "<system mess. texts unavail.>";
+#endif
+		(void)fprintf(stderr, "Sys errno: %d: %s\n",
+			my_errno, my_strerror);
+        }
+
+	/* post remainder of call data */
+	(void)vfprintf(stderr,fmt,args);
+	va_end(args);
+	/* die if code positive */
+	if (code > 0) {
+		(void)fputs("\nprogram abnormally terminated\n", stderr);
+		exit(code);
+	}
+	else
+		putc('\n', stderr);
+}
diff --git a/src/apps/emess.h b/src/apps/emess.h
new file mode 100644
index 00000000..4c6f6783
--- /dev/null
+++ b/src/apps/emess.h
@@ -0,0 +1,29 @@
+/* Error message processing header file */
+#ifndef EMESS_H
+#define EMESS_H
+
+struct EMESS {
+	char	*File_name,	/* input file name */
+			*Prog_name;	/* name of program */
+	int		File_line;	/* approximate line read
+							where error occurred */
+};
+
+#ifdef EMESS_ROUTINE	/* use type */
+/* for emess procedure */
+struct EMESS emess_dat = { nullptr, nullptr, 0 };
+
+#ifdef sun /* Archaic SunOs 4.1.1, etc. */
+extern char *sys_errlist[];
+#define strerror(n) (sys_errlist[n])
+#endif
+
+#else	/* for for calling procedures */
+
+extern struct EMESS emess_dat;
+
+#endif /* use type */
+
+void emess(int, const char *, ...);
+
+#endif /* end EMESS_H */
diff --git a/src/apps/gen_cheb.cpp b/src/apps/gen_cheb.cpp
new file mode 100644
index 00000000..4ba514d4
--- /dev/null
+++ b/src/apps/gen_cheb.cpp
@@ -0,0 +1,106 @@
+/* generates 'T' option output */
+#define PJ_LIB__
+#include "projects.h"
+#include <stdio.h>
+#include <string.h>
+#include <errno.h>
+#include "emess.h"
+#ifndef COEF_LINE_MAX
+#define COEF_LINE_MAX 50
+#endif
+
+static double strtod_type_safe(const char *s, const char ** endptr)
+{
+    char* l_endptr = nullptr;
+    double ret= strtod(s, &l_endptr);
+    if( endptr )
+        *endptr = static_cast<const char*>(l_endptr);
+    return ret;
+}
+
+static double dmstor_type_safe(const char *s, const char ** endptr)
+{
+    char* l_endptr = nullptr;
+    double ret= dmstor(s, &l_endptr);
+    if( endptr )
+        *endptr = static_cast<const char*>(l_endptr);
+    return ret;
+}
+
+static long strtol_type_safe(const char *s, const char ** endptr, int base)
+{
+    char* l_endptr = nullptr;
+    long ret = strtol(s, &l_endptr, base);
+    if( endptr )
+        *endptr = static_cast<const char*>(l_endptr);
+    return ret;
+}
+
+
+/* FIXME: put the declaration in a header. Also used in proj.c */
+void gen_cheb(int inverse, projUV (*proj)(projUV), const char *s, PJ *P,
+              int iargc, char **iargv);
+extern void p_series(Tseries *, FILE *, char *);
+
+void gen_cheb(int inverse, projUV (*proj)(projUV), const char *s, PJ *P,
+              int iargc, char **iargv) {
+	long NU = 15, NV = 15;
+	int errin = 0, pwr;
+	long res = -1;
+	char *arg, fmt[32];
+	projUV low, upp, resid;
+	Tseries *F;
+	double (*input)(const char *, const char **);
+
+	input = inverse ? strtod_type_safe : dmstor_type_safe;
+	if (*s) low.u = input(s, &s); else { low.u = 0; ++errin; }
+	if (*s == ',') upp.u = input(s+1, &s); else { upp.u = 0; ++errin; }
+	if (*s == ',') low.v = input(s+1, &s); else { low.v = 0; ++errin; }
+	if (*s == ',') upp.v = input(s+1, &s); else { upp.v = 0; ++errin; }
+	if (errin)
+		emess(16,"null or absent -T parameters");
+	if (*s == ',') if (*++s != ',') res = strtol_type_safe(s, &s, 10);
+	if (*s == ',') if (*++s != ',') NU = strtol_type_safe(s, &s, 10);
+	if (*s == ',') if (*++s != ',') NV = strtol_type_safe(s, &s, 10);
+	pwr = s && *s && !strcmp(s, ",P");
+	(void)printf("#proj_%s\n#    run-line:\n",
+		pwr ? "Power" : "Chebyshev");
+	if (iargc > 0) { /* proj execution audit trail */
+		int n = 0, L;
+
+		for ( ; iargc ; --iargc) {
+			arg = *iargv++;
+			if (*arg != '+') {
+				if (!n) { putchar('#'); ++n; }
+				(void)printf(" %s%n",arg, &L);
+				if ((n += L) > COEF_LINE_MAX) { putchar('\n'); n = 0; }
+			}
+		}
+		if (n) putchar('\n');
+	}
+	(void)printf("# projection parameters\n");
+	pj_pr_list(P);
+	if (low.u == upp.u || low.v >= upp.v)
+		emess(16,"approx. argument range error");
+	if (low.u > upp.u)
+		low.u -= M_TWOPI;
+	if (NU < 2 || NV < 2 || NU > INT_MAX || NV > INT_MAX)
+		emess(16,"approx. work dimensions (%ld %ld) too small or large",NU,NV);
+	if (!(F = mk_cheby(low, upp, pow(10., (double)res)*.5, &resid, proj,
+		(int)NU, (int)NV, pwr)))
+		emess(16,"generation of approx failed\nreason: %s\n",
+			pj_strerrno(errno));
+	(void)printf("%c,%.12g,%.12g,%.12g,%.12g,%.12g\n",inverse?'I':'F',
+		P->lam0*RAD_TO_DEG,
+		low.u*(inverse?1.:RAD_TO_DEG),upp.u*(inverse?1.:RAD_TO_DEG),
+		low.v*(inverse?1.:RAD_TO_DEG),upp.v*(inverse?1.:RAD_TO_DEG));
+	if (pwr)
+		strcpy(fmt, "%.15g");
+	else if (res <= 0)
+		(void)sprintf(fmt,"%%.%ldf",-res+1);
+	else
+		(void)strcpy(fmt,"%.0f");
+	p_series(F, stdout, fmt);
+	(void)printf("# |u,v| sums %g %g\n#end_proj_%s\n",
+		resid.u, resid.v, pwr ? "Power" : "Chebyshev");
+}
diff --git a/src/apps/geod.cpp b/src/apps/geod.cpp
new file mode 100644
index 00000000..7b6367c6
--- /dev/null
+++ b/src/apps/geod.cpp
@@ -0,0 +1,241 @@
+/* <<<< Geodesic filter program >>>> */
+
+#include "proj.h"
+# include "projects.h"
+# include "geod_interface.h"
+# include "emess.h"
+# include <ctype.h>
+# include <stdio.h>
+# include <string.h>
+
+# define MAXLINE 200
+# define MAX_PARGS 50
+# define TAB putchar('\t')
+	static int
+fullout = 0,	/* output full set of geodesic values */
+tag = '#',	/* beginning of line tag character */
+pos_azi = 0,	/* output azimuths as positive values */
+inverse = 0;	/* != 0 then inverse geodesic */
+
+static const char *oform = nullptr; /* output format for decimal degrees */
+static const char *osform = "%.3f"; /* output format for S */
+
+static char pline[50];              /* work string */
+static const char *usage =
+"%s\nusage: %s [ -afFIlptwW [args] ] [ +opts[=arg] ] [ files ]\n";
+
+	static void
+printLL(double p, double l) {
+	if (oform) {
+		(void)printf(oform, p * RAD_TO_DEG); TAB;
+		(void)printf(oform, l * RAD_TO_DEG);
+	} else {
+		(void)fputs(rtodms(pline, p, 'N', 'S'),stdout); TAB;
+		(void)fputs(rtodms(pline, l, 'E', 'W'),stdout);
+	}
+}
+	static void
+do_arc(void) {
+	double az;
+
+	printLL(phi2, lam2); putchar('\n');
+	for (az = al12; n_alpha--; ) {
+		al12 = az = adjlon(az + del_alpha);
+		geod_pre();
+		geod_for();
+		printLL(phi2, lam2); putchar('\n');
+	}
+}
+	static void	/* generate intermediate geodesic coordinates */
+do_geod(void) {
+	double phil, laml, del_S;
+
+	phil = phi2;
+	laml = lam2;
+	printLL(phi1, lam1); putchar('\n');
+	for ( geod_S = del_S = geod_S / n_S; --n_S; geod_S += del_S) {
+		geod_for();
+		printLL(phi2, lam2); putchar('\n');
+	}
+	printLL(phil, laml); putchar('\n');
+}
+	static void	/* file processing function */
+process(FILE *fid) {
+	char line[MAXLINE+3], *s;
+
+	for (;;) {
+		++emess_dat.File_line;
+		if (!(s = fgets(line, MAXLINE, fid)))
+			break;
+		if (!strchr(s, '\n')) { /* overlong line */
+			int c;
+			strcat(s, "\n");
+			/* gobble up to newline */
+			while ((c = fgetc(fid)) != EOF && c != '\n') ;
+		}
+		if (*s == tag) {
+			fputs(line, stdout);
+			continue;
+		}
+		phi1 = dmstor(s, &s);
+		lam1 = dmstor(s, &s);
+		if (inverse) {
+			phi2 = dmstor(s, &s);
+			lam2 = dmstor(s, &s);
+			geod_inv();
+		} else {
+			al12 = dmstor(s, &s);
+			geod_S = strtod(s, &s) * to_meter;
+			geod_pre();
+			geod_for();
+		}
+		if (!*s && (s > line)) --s; /* assumed we gobbled \n */
+		if (pos_azi) {
+			if (al12 < 0.) al12 += M_TWOPI;
+			if (al21 < 0.) al21 += M_TWOPI;
+		}
+		if (fullout) {
+			printLL(phi1, lam1); TAB;
+			printLL(phi2, lam2); TAB;
+			if (oform) {
+				(void)printf(oform, al12 * RAD_TO_DEG); TAB;
+				(void)printf(oform, al21 * RAD_TO_DEG); TAB;
+				(void)printf(osform, geod_S * fr_meter);
+			}  else {
+				(void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
+				(void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
+				(void)printf(osform, geod_S * fr_meter);
+			}
+		} else if (inverse)
+			if (oform) {
+				(void)printf(oform, al12 * RAD_TO_DEG); TAB;
+				(void)printf(oform, al21 * RAD_TO_DEG); TAB;
+				(void)printf(osform, geod_S * fr_meter);
+			} else {
+				(void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
+				(void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
+				(void)printf(osform, geod_S * fr_meter);
+			}
+		else {
+			printLL(phi2, lam2); TAB;
+			if (oform)
+				(void)printf(oform, al21 * RAD_TO_DEG);
+			else
+				(void)fputs(rtodms(pline, al21, 0, 0), stdout);
+		}
+		(void)fputs(s, stdout);
+	}
+}
+
+static char *pargv[MAX_PARGS];
+static int   pargc = 0;
+
+int main(int argc, char **argv) {
+	char *arg, **eargv = argv;
+	FILE *fid;
+	static int eargc = 0, c;
+
+	if ((emess_dat.Prog_name = strrchr(*argv,'/')) != nullptr) ++emess_dat.Prog_name;
+	else emess_dat.Prog_name = *argv;
+	inverse = ! strncmp(emess_dat.Prog_name, "inv", 3);
+	if (argc <= 1 ) {
+		(void)fprintf(stderr, usage, pj_get_release(),
+                              emess_dat.Prog_name);
+		exit (0);
+	}
+		/* process run line arguments */
+	while (--argc > 0) { /* collect run line arguments */
+		if(**++argv == '-') for(arg = *argv;;) {
+			switch(*++arg) {
+			case '\0': /* position of "stdin" */
+				if (arg[-1] == '-') eargv[eargc++] = const_cast<char*>("-");
+				break;
+			case 'a': /* output full set of values */
+				fullout = 1;
+				continue;
+			case 'I': /* alt. inverse spec. */
+				inverse = 1;
+				continue;
+			case 't': /* set col. one char */
+				if (arg[1]) tag = *++arg;
+				else emess(1,"missing -t col. 1 tag");
+				continue;
+			case 'W': /* specify seconds precision */
+			case 'w': /* -W for constant field width */
+				if ((c = arg[1]) && isdigit(c)) {
+					set_rtodms(c - '0', *arg == 'W');
+					++arg;
+				} else
+				    emess(1,"-W argument missing or non-digit");
+				continue;
+			case 'f': /* alternate output format degrees or xy */
+				if (--argc <= 0)
+noargument:		   emess(1,"missing argument for -%c",*arg);
+				oform = *++argv;
+				continue;
+			case 'F': /* alternate output format degrees or xy */
+				if (--argc <= 0) goto noargument;
+				osform = *++argv;
+				continue;
+			case 'l':
+				if (!arg[1] || arg[1] == 'e') { /* list of ellipsoids */
+                                    const struct PJ_ELLPS *le;
+                                    
+                                    for (le=proj_list_ellps(); le->id ; ++le)
+                                        (void)printf("%9s %-16s %-16s %s\n",
+                                                     le->id, le->major, le->ell, le->name);
+				} else if (arg[1] == 'u') { /* list of units */
+                                    const struct PJ_UNITS *lu;
+                                    
+                                    for (lu = proj_list_units();lu->id ; ++lu)
+                                        (void)printf("%12s %-20s %s\n",
+                                                     lu->id, lu->to_meter, lu->name);
+				} else
+                                    emess(1,"invalid list option: l%c",arg[1]);
+                                exit( 0 );
+			case 'p': /* output azimuths as positive */
+				pos_azi = 1;
+				continue;
+			default:
+				emess(1, "invalid option: -%c",*arg);
+				break;
+			}
+			break;
+		} else if (**argv == '+') /* + argument */
+			if (pargc < MAX_PARGS)
+				pargv[pargc++] = *argv + 1;
+			else
+				emess(1,"overflowed + argument table");
+		else /* assumed to be input file name(s) */
+			eargv[eargc++] = *argv;
+	}
+	/* done with parameter and control input */
+	geod_set(pargc, pargv); /* setup projection */
+	if ((n_alpha || n_S) && eargc)
+		emess(1,"files specified for arc/geodesic mode");
+	if (n_alpha)
+		do_arc();
+	else if (n_S)
+		do_geod();
+	else { /* process input file list */
+		if (eargc == 0) /* if no specific files force sysin */
+			eargv[eargc++] = const_cast<char*>("-");
+		for ( ; eargc-- ; ++eargv) {
+			if (**eargv == '-') {
+				fid = stdin;
+				emess_dat.File_name = const_cast<char*>("<stdin>");
+			} else {
+				if ((fid = fopen(*eargv, "r")) == nullptr) {
+					emess(-2, *eargv, "input file");
+					continue;
+				}
+				emess_dat.File_name = *eargv;
+			}
+			emess_dat.File_line = 0;
+			process(fid);
+			(void)fclose(fid);
+			emess_dat.File_name = (char *)nullptr;
+		}
+	}
+	exit(0); /* normal completion */
+}
diff --git a/src/apps/geod_interface.cpp b/src/apps/geod_interface.cpp
new file mode 100644
index 00000000..a30377ac
--- /dev/null
+++ b/src/apps/geod_interface.cpp
@@ -0,0 +1,33 @@
+#include "projects.h"
+#include "geod_interface.h"
+
+void geod_ini(void) {
+  geod_init(&GlobalGeodesic, geod_a, geod_f);
+}
+
+void geod_pre(void) {
+  double
+    lat1 = phi1 / DEG_TO_RAD, lon1 = lam1 / DEG_TO_RAD,
+    azi1 = al12 / DEG_TO_RAD;
+  geod_lineinit(&GlobalGeodesicLine, &GlobalGeodesic, lat1, lon1, azi1, 0U);
+}
+
+void geod_for(void) {
+  double
+    s12 = geod_S, lat2, lon2, azi2;
+  geod_position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
+  azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
+  phi2 = lat2 * DEG_TO_RAD;
+  lam2 = lon2 * DEG_TO_RAD;
+  al21 = azi2 * DEG_TO_RAD;
+}
+
+void geod_inv(void) {
+  double
+    lat1 = phi1 / DEG_TO_RAD, lon1 = lam1 / DEG_TO_RAD,
+    lat2 = phi2 / DEG_TO_RAD, lon2 = lam2 / DEG_TO_RAD,
+    azi1, azi2, s12;
+  geod_inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
+  azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
+  al12 = azi1 * DEG_TO_RAD; al21 = azi2 * DEG_TO_RAD; geod_S = s12;
+}
diff --git a/src/apps/geod_interface.h b/src/apps/geod_interface.h
new file mode 100644
index 00000000..255d505a
--- /dev/null
+++ b/src/apps/geod_interface.h
@@ -0,0 +1,45 @@
+#if !defined(GEOD_INTERFACE_H)
+#define GEOD_INTERFACE_H
+
+#include "geodesic.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef _IN_GEOD_SET
+#  define GEOD_EXTERN extern
+#else
+#  define GEOD_EXTERN
+#endif
+
+GEOD_EXTERN struct geodesic {
+  double A, FLAT, LAM1, PHI1, ALPHA12, LAM2, PHI2, ALPHA21, DIST;
+} GEODESIC;
+
+# define geod_a	GEODESIC.A
+# define geod_f	GEODESIC.FLAT
+# define lam1	GEODESIC.LAM1
+# define phi1	GEODESIC.PHI1
+# define al12	GEODESIC.ALPHA12
+# define lam2	GEODESIC.LAM2
+# define phi2	GEODESIC.PHI2
+# define al21	GEODESIC.ALPHA21
+# define geod_S	GEODESIC.DIST
+    
+GEOD_EXTERN struct geod_geodesic GlobalGeodesic;
+GEOD_EXTERN struct geod_geodesicline GlobalGeodesicLine;
+GEOD_EXTERN int n_alpha, n_S;
+GEOD_EXTERN double to_meter, fr_meter, del_alpha;
+	
+void geod_set(int, char **);
+void geod_ini(void);
+void geod_pre(void);
+void geod_for(void);
+void geod_inv(void);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/src/apps/geod_set.cpp b/src/apps/geod_set.cpp
new file mode 100644
index 00000000..b9e9c42f
--- /dev/null
+++ b/src/apps/geod_set.cpp
@@ -0,0 +1,75 @@
+#define _IN_GEOD_SET
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "proj.h"
+#include "projects.h"
+#include "geod_interface.h"
+#include "emess.h"
+
+	void
+geod_set(int argc, char **argv) {
+	paralist *start = nullptr, *curr;
+	double es;
+	char *name;
+	int i;
+
+    /* put arguments into internal linked list */
+	if (argc <= 0)
+		emess(1, "no arguments in initialization list");
+	start = curr = pj_mkparam(argv[0]);
+	if (!curr)
+		emess(1, "memory allocation failed");
+	for (i = 1; curr != nullptr && i < argc; ++i) {
+		curr->next = pj_mkparam(argv[i]);
+		if (!curr->next)
+			emess(1, "memory allocation failed");
+		curr = curr->next;
+	}
+	/* set elliptical parameters */
+	if (pj_ell_set(pj_get_default_ctx(),start, &geod_a, &es)) emess(1,"ellipse setup failure");
+	/* set units */
+	if ((name = pj_param(nullptr,start, "sunits").s) != nullptr) {
+		const char *s;
+                const struct PJ_UNITS *unit_list = proj_list_units();
+		for (i = 0; (s = unit_list[i].id) && strcmp(name, s) ; ++i) ;
+		if (!s)
+			emess(1,"%s unknown unit conversion id", name);
+		to_meter = unit_list[i].factor;
+		fr_meter = 1 / to_meter;
+	} else
+		to_meter = fr_meter = 1;
+	geod_f = es/(1 + sqrt(1 - es));
+	geod_ini();
+	/* check if line or arc mode */
+	if (pj_param(nullptr,start, "tlat_1").i) {
+		double del_S;
+#undef f
+		phi1 = pj_param(nullptr,start, "rlat_1").f;
+		lam1 = pj_param(nullptr,start, "rlon_1").f;
+		if (pj_param(nullptr,start, "tlat_2").i) {
+			phi2 = pj_param(nullptr,start, "rlat_2").f;
+			lam2 = pj_param(nullptr,start, "rlon_2").f;
+			geod_inv();
+			geod_pre();
+		} else if ((geod_S = pj_param(nullptr,start, "dS").f) != 0.) {
+			al12 = pj_param(nullptr,start, "rA").f;
+			geod_pre();
+			geod_for();
+		} else emess(1,"incomplete geodesic/arc info");
+		if ((n_alpha = pj_param(nullptr,start, "in_A").i) > 0) {
+			if ((del_alpha = pj_param(nullptr,start, "rdel_A").f) == 0.0)
+				emess(1,"del azimuth == 0");
+		} else if ((del_S = fabs(pj_param(nullptr,start, "ddel_S").f)) != 0.) {
+			n_S = (int)(geod_S / del_S + .5);
+		} else if ((n_S = pj_param(nullptr,start, "in_S").i) <= 0)
+			emess(1,"no interval divisor selected");
+	}
+	/* free up linked list */
+	for ( ; start; start = curr) {
+		curr = start->next;
+		pj_dalloc(start);
+	}
+}
diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp
new file mode 100644
index 00000000..3e4770a2
--- /dev/null
+++ b/src/apps/gie.cpp
@@ -0,0 +1,1458 @@
+/***********************************************************************
+
+       gie - The Geospatial Integrity Investigation Environment
+
+************************************************************************
+
+The Geospatial Integrity Investigation Environment "gie" is a modest
+regression testing environment for the PROJ.4 transformation library.
+
+Its primary design goal was to be able to replace those thousands of
+lines of regression testing code that are (at time of writing) part
+of PROJ.4, while not requiring any other kind of tooling than the same
+C compiler already employed for compiling the library.
+
+The basic functionality of the gie command language is implemented
+through just 3 command verbs:
+
+operation,     which defines the PROJ.4 operation to test,
+accept,        which defines the input coordinate to read, and
+expect,        which defines the result to expect.
+
+E.g:
+
+operation  +proj=utm  +zone=32  +ellps=GRS80
+accept     12  55
+expect     691_875.632_14   6_098_907.825_05
+
+Note that gie accepts the underscore ("_") as a thousands separator.
+It is not required (in fact, it is entirely ignored by the input
+routine), but it significantly improves the readability of the very
+long strings of numbers typically required in projected coordinates.
+
+By default, gie considers the EXPECTation met, if the result comes to
+within 0.5 mm of the expected. This default can be changed using the
+'tolerance' command verb (and yes, I know, linguistically speaking, both
+"operation" and "tolerance" are nouns, not verbs). See the first
+few hundred lines of the "builtins.gie" test file for more details of
+the command verbs available (verbs of both the VERBal and NOUNistic
+persuation).
+
+--
+
+But more importantly than being an acronym for "Geospatial Integrity
+Investigation Environment", gie were also the initials, user id, and
+USGS email address of Gerald Ian Evenden (1935--2016), the geospatial
+visionary, who, already in the 1980s, started what was to become the
+PROJ.4 of today.
+
+Gerald's clear vision was that map projections are *just special
+functions*. Some of them rather complex, most of them of two variables,
+but all of them *just special functions*, and not particularly more
+special than the sin(), cos(), tan(), and hypot() already available in
+the C standard library.
+
+And hence, according to Gerald, *they should not be particularly much
+harder to use*, for a programmer, than the sin()s, tan()s and hypot()s
+so readily available.
+
+Gerald's ingenuity also showed in the implementation of the vision,
+where he devised a comprehensive, yet simple, system of key-value
+pairs for parameterising a map projection, and the highly flexible
+PJ struct, storing run-time compiled versions of those key-value pairs,
+hence making a map projection function call, pj_fwd(PJ, point), as easy
+as a traditional function call like hypot(x,y).
+
+While today, we may have more formally well defined metadata systems
+(most prominent the OGC WKT2 representation), nothing comes close being
+as easily readable ("human compatible") as Gerald's key-value system.
+This system in particular, and the PROJ.4 system in general, was
+Gerald's great gift to anyone using and/or communicating about geodata.
+
+It is only reasonable to name a program, keeping an eye on the integrity
+of the PROJ.4 system, in honour of Gerald.
+
+So in honour, and hopefully also in the spirit, of Gerald Ian Evenden
+(1935--2016), this is the Geospatial Integrity Investigation Environment.
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
+
+************************************************************************
+
+* Copyright (c) 2017 Thomas Knudsen
+* Copyright (c) 2017, 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.
+
+***********************************************************************/
+
+#include <ctype.h>
+#include <errno.h>
+#include <math.h>
+#include <stdarg.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+#include "proj_math.h"
+#include "proj_strtod.h"
+#include "projects.h"
+
+#include "optargpm.h"
+
+/* Package for flexible format I/O - ffio */
+typedef struct ffio {
+    FILE *f;
+    const char **tags;
+    const char *tag;
+    char *args;
+    char *next_args;
+    size_t n_tags;
+    size_t args_size;
+    size_t next_args_size;
+    size_t argc;
+    size_t lineno, next_lineno;
+    size_t level;
+}  ffio;
+
+static int get_inp (ffio *G);
+static int skip_to_next_tag (ffio *G);
+static int step_into_gie_block (ffio *G);
+static int locate_tag (ffio *G, const char *tag);
+static int nextline (ffio *G);
+static int at_end_delimiter (ffio *G);
+static const char *at_tag (ffio *G);
+static int at_decorative_element (ffio *G);
+static ffio *ffio_destroy (ffio *G);
+static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_size);
+
+static const char *gie_tags[] = {
+    "<gie>", "operation", "use_proj4_init_rules",
+    "accept", "expect", "roundtrip", "banner", "verbose",
+    "direction", "tolerance", "ignore", "require_grid", "echo", "skip", "</gie>"
+};
+
+static const size_t n_gie_tags = sizeof gie_tags / sizeof gie_tags[0];
+
+
+int   main(int argc, char **argv);
+
+static int   dispatch (const char *cmnd, const char *args);
+static int   errmsg (int errlev, const char *msg, ...);
+static int   errno_from_err_const (const char *err_const);
+static int   list_err_codes (void);
+static int   process_file (const char *fname);
+
+static const char *column (const char *buf, int n);
+static const char *err_const_from_errno (int err);
+
+
+
+#define SKIP -1
+
+#define MAX_OPERATION 10000
+
+typedef struct {
+    char operation[MAX_OPERATION+1];
+    PJ *P;
+    PJ_COORD a, b, c, e;
+    PJ_DIRECTION dir;
+    int verbosity;
+    int skip;
+    int op_id;
+    int op_ok,    op_ko,    op_skip;
+    int total_ok, total_ko, total_skip;
+    int grand_ok, grand_ko, grand_skip;
+    size_t operation_lineno;
+    size_t dimensions_given, dimensions_given_at_last_accept;
+    double tolerance;
+    int use_proj4_init_rules;
+    int ignore;
+    int skip_test;
+    const char *curr_file;
+    FILE *fout;
+} gie_ctx;
+
+ffio *F = nullptr;
+
+static gie_ctx T;
+int tests=0, succs=0, succ_fails=0, fail_fails=0, succ_rtps=0, fail_rtps=0;
+
+static const char delim[] = {"-------------------------------------------------------------------------------\n"};
+
+
+static const char usage[] = {
+    "--------------------------------------------------------------------------------\n"
+    "Usage: %s [-options]... infile...\n"
+    "--------------------------------------------------------------------------------\n"
+    "Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    -h                Help: print this usage information\n"
+    "    -o /path/to/file  Specify output file name\n"
+    "    -v                Verbose: Provide non-essential informational output.\n"
+    "                      Repeat -v for more verbosity (e.g. -vv)\n"
+    "    -q                Quiet: Opposite of verbose. In quiet mode not even errors\n"
+    "                      are reported. Only interaction is through the return code\n"
+    "                      (0 on success, non-zero indicates number of FAILED tests)\n"
+    "    -l                List the PROJ internal system error codes\n"
+    "--------------------------------------------------------------------------------\n"
+    "Long Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    --output          Alias for -o\n"
+    "    --verbose         Alias for -v\n"
+    "    --help            Alias for -h\n"
+    "    --list            Alias for -l\n"
+    "    --version         Print version number\n"
+    "--------------------------------------------------------------------------------\n"
+    "Examples:\n"
+    "--------------------------------------------------------------------------------\n"
+    "1. Run all tests in file \"corner-cases.gie\", providing much extra information\n"
+    "       gie -vvvv corner-cases.gie\n"
+    "2. Run all tests in files \"foo\" and \"bar\", providing info on failures only\n"
+    "       gie foo bar\n"
+    "--------------------------------------------------------------------------------\n"
+};
+
+int main (int argc, char **argv) {
+    int  i;
+    const char *longflags[]  = {"v=verbose", "q=quiet", "h=help", "l=list", "version", nullptr};
+    const char *longkeys[]   = {"o=output", nullptr};
+    OPTARGS *o;
+
+    memset (&T, 0, sizeof (T));
+    T.dir = PJ_FWD;
+    T.verbosity = 1;
+    T.tolerance = 5e-4;
+    T.ignore = 5555; /* Error code that will not be issued by proj_create() */
+    T.use_proj4_init_rules = FALSE;
+
+    o = opt_parse (argc, argv, "hlvq", "o", longflags, longkeys);
+    if (nullptr==o)
+        return 0;
+
+    if (opt_given (o, "h") || argc==1) {
+        printf (usage, o->progname);
+        free (o);
+        return 0;
+    }
+
+
+    if (opt_given (o, "version")) {
+        fprintf (stdout, "%s: %s\n", o->progname, pj_get_release ());
+        free (o);
+        return 0;
+    }
+
+    T.verbosity = opt_given (o, "q");
+    if (T.verbosity)
+        T.verbosity = -1;
+    if (T.verbosity != -1)
+        T.verbosity = opt_given (o, "v") + 1;
+
+    T.fout = stdout;
+    if (opt_given (o, "o"))
+        T.fout = fopen (opt_arg (o, "output"), "rt");
+
+    if (nullptr==T.fout) {
+        fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
+        free (o);
+        return 1;
+    }
+
+    if (opt_given (o, "l")) {
+        free (o);
+        return list_err_codes ();
+    }
+
+    if (0==o->fargc) {
+        if (T.verbosity==-1)
+            return -1;
+        fprintf (T.fout, "Nothing to do\n");
+        free (o);
+        return 0;
+    }
+
+    F = ffio_create (gie_tags, n_gie_tags, 1000);
+    if (nullptr==F) {
+        fprintf (stderr, "%s: No memory\n", o->progname);
+        free (o);
+        return 1;
+    }
+
+    for (i = 0;  i < o->fargc;  i++)
+        process_file (o->fargv[i]);
+
+    if (T.verbosity > 0) {
+        if (o->fargc > 1) {
+            fprintf (T.fout, "%sGrand total: %d. Success: %d, Skipped: %d, Failure: %d\n",
+                     delim, T.grand_ok+T.grand_ko+T.grand_skip, T.grand_ok, T.grand_skip,
+                     T.grand_ko);
+        }
+        fprintf (T.fout, "%s", delim);
+        if (T.verbosity > 1) {
+            fprintf (T.fout, "Failing roundtrips: %4d,    Succeeding roundtrips: %4d\n", fail_rtps, succ_rtps);
+            fprintf (T.fout, "Failing failures:   %4d,    Succeeding failures:   %4d\n", fail_fails, succ_fails);
+            fprintf (T.fout, "Internal counters:                            %4.4d(%4.4d)\n", tests, succs);
+            fprintf (T.fout, "%s", delim);
+        }
+    }
+    else
+        if (T.grand_ko)
+            fprintf (T.fout, "Failures: %d", T.grand_ko);
+
+    if (stdout != T.fout)
+        fclose (T.fout);
+
+    free (o);
+    ffio_destroy (F);
+    return T.grand_ko;
+}
+
+static int another_failure (void) {
+    T.op_ko++;
+    T.total_ko++;
+    proj_errno_reset (T.P);
+    return 0;
+}
+
+static int another_skip (void) {
+    T.op_skip++;
+    T.total_skip++;
+    return 0;
+}
+
+static int another_success (void) {
+    T.op_ok++;
+    T.total_ok++;
+    proj_errno_reset (T.P);
+    return 0;
+}
+
+static int another_succeeding_failure (void) {
+    succ_fails++;
+    return another_success ();
+}
+
+static int another_failing_failure (void) {
+    fail_fails++;
+    return another_failure ();
+}
+
+static int another_succeeding_roundtrip (void) {
+    succ_rtps++;
+    return another_success ();
+}
+
+static int another_failing_roundtrip (void) {
+    fail_rtps++;
+    return another_failure ();
+}
+
+static int process_file (const char *fname) {
+    FILE *f;
+
+    F->lineno = F->next_lineno = F->level = 0;
+    T.op_ok = T.total_ok = 0;
+    T.op_ko = T.total_ko = 0;
+    T.op_skip = T.total_skip = 0;
+
+    if (T.skip) {
+        proj_destroy (T.P);
+        T.P = nullptr;
+        return 0;
+    }
+
+    f = fopen (fname, "rt");
+    if (nullptr==f) {
+        if (T.verbosity > 0) {
+            fprintf (T.fout, "%sCannot open spec'd input file '%s' - bye!\n", delim, fname);
+            return 2;
+        }
+        errmsg (2, "Cannot open spec'd input file '%s' - bye!\n", fname);
+    }
+    F->f = f;
+
+    if (T.verbosity > 0)
+        fprintf (T.fout, "%sReading file '%s'\n", delim, fname);
+    T.curr_file = fname;
+
+    while (get_inp(F)) {
+        if (SKIP==dispatch (F->tag, F->args)) {
+            proj_destroy (T.P);
+            T.P = nullptr;
+            return 0;
+        }
+    }
+
+    fclose (f);
+    F->lineno = F->next_lineno = 0;
+
+    T.grand_ok   += T.total_ok;
+    T.grand_ko   += T.total_ko;
+    T.grand_skip += T.grand_skip;
+    if (T.verbosity > 0) {
+        fprintf (T.fout, "%stotal: %2d tests succeeded, %2d tests skipped, %2d tests %s\n",
+                 delim, T.total_ok, T.total_skip, T.total_ko,
+                 T.total_ko? "FAILED!": "failed.");
+    }
+    if (F->level==0)
+        return errmsg (-3, "File '%s':Missing '<gie>' cmnd - bye!\n", fname);
+    if (F->level && F->level%2)
+        return errmsg (-4, "File '%s':Missing '</gie>' cmnd - bye!\n", fname);
+    return 0;
+}
+
+
+/*****************************************************************************/
+const char *column (const char *buf, int n) {
+/*****************************************************************************
+Return a pointer to the n'th column of buf. Column numbers start at 0.
+******************************************************************************/
+    int i;
+    if (n <= 0)
+        return buf;
+    for (i = 0;  i < n;  i++) {
+        while (isspace(*buf))
+            buf++;
+        if (i == n - 1)
+            break;
+        while ((0 != *buf) && !isspace(*buf))
+            buf++;
+    }
+    return buf;
+}
+
+
+/*****************************************************************************/
+static double strtod_scaled (const char *args, double default_scale) {
+/*****************************************************************************
+Interpret <args> as a numeric followed by a linear decadal prefix.
+Return the properly scaled numeric
+******************************************************************************/
+    double s;
+    const double GRS80_DEG = 111319.4908; /* deg-to-m at equator of GRS80 */
+    const char *endp = args;
+    s = proj_strtod (args, (char **) &endp);
+    if (args==endp)
+        return HUGE_VAL;
+
+    endp = column (args, 2);
+
+    if (0==strcmp(endp, "km"))
+        s *= 1000;
+    else if (0==strcmp(endp, "m"))
+        s *= 1;
+    else if (0==strcmp(endp, "dm"))
+        s /= 10;
+    else if (0==strcmp(endp, "cm"))
+        s /= 100;
+    else if (0==strcmp(endp, "mm"))
+        s /= 1000;
+    else if (0==strcmp(endp, "um"))
+        s /= 1e6;
+    else if (0==strcmp(endp, "nm"))
+        s /= 1e9;
+    else if (0==strcmp(endp, "rad"))
+        s = GRS80_DEG * proj_todeg (s);
+    else if (0==strcmp(endp, "deg"))
+        s = GRS80_DEG * s;
+    else
+        s *= default_scale;
+    return s;
+}
+
+
+static int banner (const char *args) {
+    char dots[] = {"..."}, nodots[] = {""}, *thedots = nodots;
+    if (strlen(args) > 70)
+        thedots = dots;
+    fprintf (T.fout, "%s%-70.70s%s\n", delim, args, thedots);
+    return 0;
+}
+
+
+static int tolerance (const char *args) {
+    T.tolerance = strtod_scaled (args, 1);
+    if (HUGE_VAL==T.tolerance) {
+        T.tolerance = 0.0005;
+        return 1;
+    }
+    return 0;
+}
+
+
+static int use_proj4_init_rules (const char *args) {
+    T.use_proj4_init_rules = strcmp(args, "true") == 0;
+    return 0;
+}
+
+static int ignore (const char *args) {
+    T.ignore = errno_from_err_const (column (args, 1));
+    return 0;
+}
+
+static int require_grid (const char *args) {
+    PJ_GRID_INFO grid_info;
+    const char* grid_filename = column (args, 1);
+    grid_info = proj_grid_info(grid_filename);
+    if( strlen(grid_info.filename) == 0 ) {
+        if (T.verbosity > 1) {
+            fprintf (T.fout, "Test skipped because of missing grid %s\n",
+                    grid_filename);
+        }
+        T.skip_test = 1;
+    }
+    return 0;
+}
+
+static int direction (const char *args) {
+    const char *endp = args;
+    while (isspace (*endp))
+        endp++;
+    switch (*endp) {
+        case 'F':
+        case 'f':
+        T.dir = PJ_FWD;
+            break;
+        case 'I':
+        case 'i':
+        case 'R':
+        case 'r':
+            T.dir = PJ_INV;
+            break;
+        default:
+            return 1;
+    }
+
+    return 0;
+}
+
+
+static void finish_previous_operation (const char *args) {
+    if (T.verbosity > 1 && T.op_id > 1 && T.op_ok+T.op_ko)
+        fprintf (T.fout, "%s     %d tests succeeded,  %d tests skipped, %d tests %s\n",
+                 delim, T.op_ok, T.op_skip, T.op_ko, T.op_ko? "FAILED!": "failed.");
+    (void) args;
+}
+
+
+
+/*****************************************************************************/
+static int operation (char *args) {
+/*****************************************************************************
+Define the operation to apply to the input data (in ISO 19100 lingo,
+an operation is the general term describing something that can be
+either a conversion or a transformation)
+******************************************************************************/
+    T.op_id++;
+
+    T.operation_lineno = F->lineno;
+
+    strncpy (&(T.operation[0]), F->args, MAX_OPERATION);
+    T.operation[MAX_OPERATION] = '\0';
+
+    if (T.verbosity > 1) {
+        finish_previous_operation (F->args);
+        banner (args);
+    }
+
+
+    T.op_ok = 0;
+    T.op_ko = 0;
+    T.op_skip = 0;
+    T.skip_test = 0;
+
+    direction ("forward");
+    tolerance ("0.5 mm");
+    ignore ("pjd_err_dont_skip");
+
+    proj_errno_reset (T.P);
+
+    if (T.P)
+        proj_destroy (T.P);
+    proj_errno_reset (nullptr);
+    proj_context_use_proj4_init_rules(nullptr, T.use_proj4_init_rules);
+
+    T.P = proj_create (nullptr, F->args);
+
+    /* Checking that proj_create succeeds is first done at "expect" time, */
+    /* since we want to support "expect"ing specific error codes */
+
+    return 0;
+}
+
+static PJ_COORD torad_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    const char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_INV)
+        axis = l->param + strlen ("axis=");
+    n = strlen (axis);
+    for (i = 0;  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_torad (a.v[i]);
+    return a;
+}
+
+
+static PJ_COORD todeg_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    const char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_FWD)
+        axis = l->param + strlen ("axis=");
+    n = strlen (axis);
+    for (i = 0;  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_todeg (a.v[i]);
+    return a;
+}
+
+
+
+/*****************************************************************************/
+static PJ_COORD parse_coord (const char *args) {
+/*****************************************************************************
+Attempt to interpret args as a PJ_COORD.
+******************************************************************************/
+    int i;
+    const char *endp;
+    const char *dmsendp;
+    const char *prev = args;
+    PJ_COORD a = proj_coord (0,0,0,0);
+
+    T.dimensions_given = 0;
+    for (i = 0;   i < 4;   i++) {
+        /* proj_strtod doesn't read values like 123d45'678W so we need a bit */
+        /* of help from proj_dmstor. proj_strtod effectively ignores what    */
+        /* comes after "d", so we use that fact that when dms is larger than */
+        /* d the value was stated in "dms" form.                             */
+        /* This could be avoided if proj_dmstor used the same proj_strtod()  */
+        /* as gie, but that is not the case (yet). When we remove projects.h */
+        /* from the public API we can change that.                           */
+        double d = proj_strtod(prev,  (char **) &endp);
+        double dms = PJ_TODEG(proj_dmstor (prev, (char **) &dmsendp));
+       /* TODO: When projects.h is removed, call proj_dmstor() in all cases */
+        if (d != dms && fabs(d) < fabs(dms) && fabs(dms) < fabs(d) + 1) {
+            d = dms;
+            endp = dmsendp;
+        }
+        /* A number like -81d00'00.000 will be parsed correctly by both */
+        /* proj_strtod and proj_dmstor but only the latter will return  */
+        /* the correct end-pointer.                                     */
+        if (d == dms && endp != dmsendp)
+            endp = dmsendp;
+
+        /* Break out if there were no more numerals */
+        if (prev==endp)
+            return i > 1? a: proj_coord_error ();
+
+        a.v[i] = d;
+        prev = endp;
+        T.dimensions_given++;
+    }
+
+    return a;
+}
+
+
+/*****************************************************************************/
+static int accept (const char *args) {
+/*****************************************************************************
+Read ("ACCEPT") a 2, 3, or 4 dimensional input coordinate.
+******************************************************************************/
+    T.a = parse_coord (args);
+    if (T.verbosity > 3)
+        fprintf (T.fout, "#  %s\n", args);
+    T.dimensions_given_at_last_accept = T.dimensions_given;
+    return 0;
+}
+
+
+/*****************************************************************************/
+static int roundtrip (const char *args) {
+/*****************************************************************************
+Check how far we go from the ACCEPTed point when doing successive
+back/forward transformation pairs.
+
+Without args, roundtrip defaults to 100 iterations:
+
+  roundtrip
+
+With one arg, roundtrip will default to a tolerance of T.tolerance:
+
+  roundtrip ntrips
+
+With two args:
+
+  roundtrip ntrips tolerance
+
+Always returns 0.
+******************************************************************************/
+    int ntrips;
+    double d, r, ans;
+    char *endp;
+    PJ_COORD coo;
+
+    if (nullptr==T.P) {
+        if (T.ignore == proj_errno(T.P))
+            return another_skip();
+
+        return another_failure ();
+    }
+
+    ans = proj_strtod (args, &endp);
+    if (endp==args) {
+        /* Default to 100 iterations if not args. */
+        ntrips = 100;
+    } else {
+        if (ans < 1.0 || ans > 1000000.0) {
+            errmsg (2, "Invalid number of roundtrips: %lf\n", ans);
+            return another_failing_roundtrip ();
+        }
+        ntrips = (int)ans;
+    }
+
+    d = strtod_scaled (endp, 1);
+    d = d==HUGE_VAL?  T.tolerance:  d;
+
+    /* input ("accepted") values - probably in degrees */
+    coo = proj_angular_input  (T.P, T.dir)? torad_coord (T.P, T.dir, T.a):  T.a;
+
+    r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
+    if (r <= d)
+        return another_succeeding_roundtrip ();
+
+    if (T.verbosity > -1) {
+        if (0==T.op_ko && T.verbosity < 2)
+            banner (T.operation);
+        fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+        fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) F->lineno);
+        fprintf (T.fout, "     roundtrip deviation: %.6f mm, expected: %.6f mm\n", 1000*r, 1000*d);
+    }
+    return another_failing_roundtrip ();
+}
+
+
+static int expect_message (double d, const char *args) {
+    another_failure ();
+
+    if (T.verbosity < 0)
+        return 1;
+    if (d > 1e6)
+        d = 999999.999999;
+    if (0==T.op_ko && T.verbosity < 2)
+        banner (T.operation);
+    fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+
+    fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) F->lineno);
+    fprintf (T.fout, "     expected: %s\n", args);
+    fprintf (T.fout, "     got:      %.12f   %.12f", T.b.xy.x,  T.b.xy.y);
+    if (T.b.xyzt.t!=0 || T.b.xyzt.z!=0)
+        fprintf (T.fout, "   %.9f", T.b.xyz.z);
+    if (T.b.xyzt.t!=0)
+        fprintf (T.fout, "   %.9f", T.b.xyzt.t);
+    fprintf (T.fout, "\n");
+    fprintf (T.fout, "     deviation:  %.6f mm,  expected:  %.6f mm\n", 1000*d, 1000*T.tolerance);
+    return 1;
+}
+
+
+static int expect_message_cannot_parse (const char *args) {
+    another_failure ();
+    if (T.verbosity > -1) {
+        if (0==T.op_ko && T.verbosity < 2)
+            banner (T.operation);
+        fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+        fprintf (T.fout, "     FAILURE in %s(%d):\n     Too few args: %s\n", opt_strip_path (T.curr_file), (int) F->lineno, args);
+    }
+    return 1;
+}
+
+static int expect_failure_with_errno_message (int expected, int got) {
+    another_failing_failure ();
+
+    if (T.verbosity < 0)
+        return 1;
+    if (0==T.op_ko && T.verbosity < 2)
+        banner (T.operation);
+    fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+    fprintf (T.fout, "     FAILURE in %s(%d):\n",    opt_strip_path (T.curr_file), (int) F->lineno);
+    fprintf (T.fout, "     got errno %s (%d): %s\n", err_const_from_errno(got), got,  pj_strerrno (got));
+    fprintf (T.fout, "     expected %s (%d):  %s",   err_const_from_errno(expected), expected, pj_strerrno (expected));
+    fprintf (T.fout, "\n");
+    return 1;
+}
+
+
+/* For test purposes, we want to call a transformation of the same */
+/* dimensionality as the number of dimensions given in accept */
+static PJ_COORD expect_trans_n_dim (PJ_COORD ci) {
+    if (4==T.dimensions_given_at_last_accept)
+        return proj_trans (T.P, T.dir, ci);
+
+    if (3==T.dimensions_given_at_last_accept)
+        return pj_approx_3D_trans (T.P, T.dir, ci);
+
+    return pj_approx_2D_trans (T.P, T.dir, ci);
+}
+
+
+/*****************************************************************************/
+static int expect (const char *args) {
+/*****************************************************************************
+Tell GIE what to expect, when transforming the ACCEPTed input
+******************************************************************************/
+    PJ_COORD ci, co, ce;
+    double d;
+    int expect_failure = 0;
+    int expect_failure_with_errno = 0;
+
+    if (0==strncmp (args, "failure", 7)) {
+        expect_failure = 1;
+
+        /* Option: Fail with an expected errno (syntax: expect failure errno -33) */
+        if (0==strncmp (column (args, 2), "errno", 5))
+            expect_failure_with_errno = errno_from_err_const (column (args, 3));
+    }
+
+    if (T.ignore==proj_errno(T.P))
+        return another_skip ();
+
+    if (nullptr==T.P) {
+        /* If we expect failure, and fail, then it's a success... */
+        if (expect_failure) {
+            /* Failed to fail correctly? */
+            if (expect_failure_with_errno && proj_errno (T.P)!=expect_failure_with_errno)
+                return expect_failure_with_errno_message (expect_failure_with_errno, proj_errno(T.P));
+
+            return another_succeeding_failure ();
+        }
+
+        /* Otherwise, it's a true failure */
+        banner (T.operation);
+        errmsg (3, "%sInvalid operation definition in line no. %d:\n       %s (errno=%s/%d)\n",
+            delim, (int) T.operation_lineno, pj_strerrno(proj_errno(T.P)),
+            err_const_from_errno (proj_errno(T.P)), proj_errno(T.P)
+        );
+        return another_failing_failure ();
+    }
+
+    /* We may still successfully fail even if the proj_create succeeded */
+    if (expect_failure) {
+        proj_errno_reset (T.P);
+
+        /* Try to carry out the operation - and expect failure */
+        ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
+        co = expect_trans_n_dim (ci);
+
+        if (expect_failure_with_errno) {
+            if (proj_errno (T.P)==expect_failure_with_errno)
+                return another_succeeding_failure ();
+            fprintf (T.fout, "errno=%d, expected=%d\n", proj_errno (T.P), expect_failure_with_errno);
+            return another_failing_failure ();
+        }
+
+
+        /* Succeeded in failing? - that's a success */
+        if (co.xyz.x==HUGE_VAL)
+            return another_succeeding_failure ();
+
+        /* Failed to fail? - that's a failure */
+        banner (T.operation);
+        errmsg (3, "%sFailed to fail. Operation definition in line no. %d\n",
+            delim, (int) T.operation_lineno
+        );
+        return another_failing_failure ();
+    }
+
+
+    if (T.verbosity > 3) {
+        fprintf (T.fout, "%s\n", T.P->inverted? "INVERTED": "NOT INVERTED");
+        fprintf (T.fout, "%s\n", T.dir== 1? "forward": "reverse");
+        fprintf (T.fout, "%s\n", proj_angular_input (T.P, T.dir)?  "angular in":  "linear in");
+        fprintf (T.fout, "%s\n", proj_angular_output (T.P, T.dir)? "angular out": "linear out");
+        fprintf (T.fout, "left: %d   right:  %d\n", T.P->left, T.P->right);
+    }
+
+    tests++;
+    T.e  =  parse_coord (args);
+    if (HUGE_VAL==T.e.v[0])
+        return expect_message_cannot_parse (args);
+
+
+    /* expected angular values, probably in degrees */
+    ce = proj_angular_output (T.P, T.dir)? torad_coord (T.P, T.dir, T.e): T.e;
+    if (T.verbosity > 3)
+        fprintf (T.fout, "EXPECTS  %.12f  %.12f  %.12f  %.12f\n",
+                 ce.v[0],ce.v[1],ce.v[2],ce.v[3]);
+
+    /* input ("accepted") values, also probably in degrees */
+    ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
+    if (T.verbosity > 3)
+        fprintf (T.fout, "ACCEPTS  %.12f  %.12f  %.12f  %.12f\n",
+                 ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
+
+    /* do the transformation, but mask off dimensions not given in expect-ation */
+    co = expect_trans_n_dim (ci);
+    if (T.dimensions_given < 4)
+        co.v[3] = 0;
+    if (T.dimensions_given < 3)
+        co.v[2] = 0;
+
+    /* angular output from proj_trans comes in radians */
+    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (T.P, T.dir, co): co;
+    if (T.verbosity > 3)
+        fprintf (T.fout, "GOT      %.12f  %.12f  %.12f  %.12f\n",
+                 co.v[0],co.v[1],co.v[2],co.v[3]);
+
+#if 0
+    /* We need to handle unusual axis orders - that'll be an item for version 5.1 */
+    if (T.P->axisswap) {
+        ce = proj_trans (T.P->axisswap, T.dir, ce);
+        co = proj_trans (T.P->axisswap, T.dir, co);
+    }
+#endif
+    if (proj_angular_output (T.P, T.dir))
+        d = proj_lpz_dist (T.P, ce, co);
+    else
+        d = proj_xyz_dist (co, ce);
+
+    if (d > T.tolerance)
+        return expect_message (d, args);
+    succs++;
+
+    another_success ();
+    return 0;
+}
+
+
+
+/*****************************************************************************/
+static int verbose (const char *args) {
+/*****************************************************************************
+Tell the system how noisy it should be
+******************************************************************************/
+    int i = (int) proj_atof (args);
+
+    /* if -q/--quiet flag has been given, we do nothing */
+    if (T.verbosity < 0)
+        return 0;
+
+    if (strlen (args))
+        T.verbosity = i;
+    else
+        T.verbosity++;
+    return 0;
+}
+
+
+/*****************************************************************************/
+static int echo (const char *args) {
+/*****************************************************************************
+Add user defined noise to the output stream
+******************************************************************************/
+fprintf (T.fout, "%s\n", args);
+    return 0;
+}
+
+
+
+/*****************************************************************************/
+static int skip (const char *args) {
+/*****************************************************************************
+Indicate that the remaining material should be skipped. Mostly for debugging.
+******************************************************************************/
+    T.skip = 1;
+    (void) args;
+    F->level = 2; /* Silence complaints about missing </gie> element */
+    return 0;
+}
+
+
+static int dispatch (const char *cmnd, const char *args) {
+    if (T.skip)
+        return SKIP;
+    if  (0==strcmp (cmnd, "operation")) return  operation ((char *) args);
+    if (T.skip_test)
+    {
+        if  (0==strcmp (cmnd, "expect"))    return  another_skip();
+        return 0;
+    }
+    if  (0==strcmp (cmnd, "accept"))    return  accept    (args);
+    if  (0==strcmp (cmnd, "expect"))    return  expect    (args);
+    if  (0==strcmp (cmnd, "roundtrip")) return  roundtrip (args);
+    if  (0==strcmp (cmnd, "banner"))    return  banner    (args);
+    if  (0==strcmp (cmnd, "verbose"))   return  verbose   (args);
+    if  (0==strcmp (cmnd, "direction")) return  direction (args);
+    if  (0==strcmp (cmnd, "tolerance")) return  tolerance (args);
+    if  (0==strcmp (cmnd, "ignore"))    return  ignore    (args);
+    if  (0==strcmp (cmnd, "require_grid"))   return  require_grid   (args);
+    if  (0==strcmp (cmnd, "echo"))      return  echo      (args);
+    if  (0==strcmp (cmnd, "skip"))      return  skip      (args);
+    if  (0==strcmp (cmnd, "use_proj4_init_rules"))
+                                        return  use_proj4_init_rules    (args);
+
+    return 0;
+}
+
+
+
+
+namespace { // anonymous namespace
+struct errno_vs_err_const {const char *the_err_const; int the_errno;};
+static const struct errno_vs_err_const lookup[] = {
+    {"pjd_err_no_args"                  ,  -1},
+    {"pjd_err_no_option_in_init_file"   ,  -2},
+    {"pjd_err_no_colon_in_init_string"  ,  -3},
+    {"pjd_err_proj_not_named"           ,  -4},
+    {"pjd_err_unknown_projection_id"    ,  -5},
+    {"pjd_err_eccentricity_is_one"      ,  -6},
+    {"pjd_err_unknown_unit_id"          ,  -7},
+    {"pjd_err_invalid_boolean_param"    ,  -8},
+    {"pjd_err_unknown_ellp_param"       ,  -9},
+    {"pjd_err_rev_flattening_is_zero"   ,  -10},
+    {"pjd_err_ref_rad_larger_than_90"   ,  -11},
+    {"pjd_err_es_less_than_zero"        ,  -12},
+    {"pjd_err_major_axis_not_given"     ,  -13},
+    {"pjd_err_lat_or_lon_exceed_limit"  ,  -14},
+    {"pjd_err_invalid_x_or_y"           ,  -15},
+    {"pjd_err_wrong_format_dms_value"   ,  -16},
+    {"pjd_err_non_conv_inv_meri_dist"   ,  -17},
+    {"pjd_err_non_con_inv_phi2"         ,  -18},
+    {"pjd_err_acos_asin_arg_too_large"  ,  -19},
+    {"pjd_err_tolerance_condition"      ,  -20},
+    {"pjd_err_conic_lat_equal"          ,  -21},
+    {"pjd_err_lat_larger_than_90"       ,  -22},
+    {"pjd_err_lat1_is_zero"             ,  -23},
+    {"pjd_err_lat_ts_larger_than_90"    ,  -24},
+    {"pjd_err_control_point_no_dist"    ,  -25},
+    {"pjd_err_no_rotation_proj"         ,  -26},
+    {"pjd_err_w_or_m_zero_or_less"      ,  -27},
+    {"pjd_err_lsat_not_in_range"        ,  -28},
+    {"pjd_err_path_not_in_range"        ,  -29},
+    {"pjd_err_h_less_than_zero"         ,  -30},
+    {"pjd_err_k_less_than_zero"         ,  -31},
+    {"pjd_err_lat_1_or_2_zero_or_90"    ,  -32},
+    {"pjd_err_lat_0_or_alpha_eq_90"     ,  -33},
+    {"pjd_err_ellipsoid_use_required"   ,  -34},
+    {"pjd_err_invalid_utm_zone"         ,  -35},
+    {"pjd_err_tcheby_val_out_of_range"  ,  -36},
+    {"pjd_err_failed_to_find_proj"      ,  -37},
+    {"pjd_err_failed_to_load_grid"      ,  -38},
+    {"pjd_err_invalid_m_or_n"           ,  -39},
+    {"pjd_err_n_out_of_range"           ,  -40},
+    {"pjd_err_lat_1_2_unspecified"      ,  -41},
+    {"pjd_err_abs_lat1_eq_abs_lat2"     ,  -42},
+    {"pjd_err_lat_0_half_pi_from_mean"  ,  -43},
+    {"pjd_err_unparseable_cs_def"       ,  -44},
+    {"pjd_err_geocentric"               ,  -45},
+    {"pjd_err_unknown_prime_meridian"   ,  -46},
+    {"pjd_err_axis"                     ,  -47},
+    {"pjd_err_grid_area"                ,  -48},
+    {"pjd_err_invalid_sweep_axis"       ,  -49},
+    {"pjd_err_malformed_pipeline"       ,  -50},
+    {"pjd_err_unit_factor_less_than_0"  ,  -51},
+    {"pjd_err_invalid_scale"            ,  -52},
+    {"pjd_err_non_convergent"           ,  -53},
+    {"pjd_err_missing_args"             ,  -54},
+    {"pjd_err_lat_0_is_zero"            ,  -55},
+    {"pjd_err_ellipsoidal_unsupported"  ,  -56},
+    {"pjd_err_too_many_inits"           ,  -57},
+    {"pjd_err_invalid_arg"              ,  -58},
+    {"pjd_err_dont_skip"                ,  5555},
+    {"pjd_err_unknown"                  ,  9999},
+    {"pjd_err_enomem"                   ,  ENOMEM},
+};
+} // anonymous namespace
+
+static const struct errno_vs_err_const unknown = {"PJD_ERR_UNKNOWN", 9999};
+
+
+static int list_err_codes (void) {
+    int i;
+    const int n = sizeof lookup / sizeof lookup[0];
+
+    for (i = 0;  i < n;  i++) {
+        if (9999==lookup[i].the_errno)
+            break;
+        fprintf (T.fout, "%25s  (%2.2d):  %s\n", lookup[i].the_err_const + 8,
+                 lookup[i].the_errno, pj_strerrno(lookup[i].the_errno));
+    }
+    return 0;
+}
+
+
+static const char *err_const_from_errno (int err) {
+    size_t i;
+    const size_t n = sizeof lookup / sizeof lookup[0];
+
+    for (i = 0;  i < n;  i++) {
+        if (err==lookup[i].the_errno)
+            return lookup[i].the_err_const + 8;
+    }
+    return unknown.the_err_const;
+}
+
+
+static int errno_from_err_const (const char *err_const) {
+    const size_t n = sizeof lookup / sizeof lookup[0];
+    size_t i, len;
+    int ret;
+    char tolower_err_const[100];
+
+    /* Make a lower case copy for matching */
+    for (i = 0;  i < 99; i++) {
+        if (0==err_const[i] || isspace (err_const[i]))
+             break;
+        tolower_err_const[i] = (char) tolower (err_const[i]);
+    }
+    tolower_err_const[i] = 0;
+
+    /* If it looks numeric, return that numeric */
+    ret = (int) pj_atof (err_const);
+    if (0!=ret)
+        return ret;
+
+    /* Else try to find a matching identifier */
+    len = strlen (tolower_err_const);
+
+    /* First try to find a match excluding the PJD_ERR_ prefix */
+    for (i = 0;  i < n;  i++) {
+        if (0==strncmp (lookup[i].the_err_const + 8, err_const, len))
+            return lookup[i].the_errno;
+    }
+
+    /* If that did not work, try with the full name */
+    for (i = 0;  i < n;  i++) {
+        if (0==strncmp (lookup[i].the_err_const, err_const, len))
+            return lookup[i].the_errno;
+    }
+
+    /* On failure, return something unlikely */
+    return 9999;
+}
+
+
+static int errmsg (int errlev, const char *msg, ...) {
+    va_list args;
+    va_start(args, msg);
+    vfprintf(stdout, msg, args);
+    va_end(args);
+    if (errlev)
+        errno = errlev;
+    return errlev;
+}
+
+
+
+
+
+
+
+
+/****************************************************************************************
+
+FFIO - Flexible format I/O
+
+FFIO provides functionality for reading proj style instruction strings written
+in a less strict format than usual:
+
+*  Whitespace is generally allowed everywhere
+*  Comments can be written inline, '#' style
+*  ... or as free format blocks
+
+The overall mission of FFIO is to facilitate communications of geodetic
+parameters and test material in a format that is highly human readable,
+and provides ample room for comment, documentation, and test material.
+
+See the PROJ ".gie" test suites for examples of supported formatting.
+
+****************************************************************************************/
+
+
+/***************************************************************************************/
+static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_size) {
+/****************************************************************************************
+Constructor for the ffio object.
+****************************************************************************************/
+    ffio *G = static_cast<ffio*>(calloc (1, sizeof (ffio)));
+    if (nullptr==G)
+        return nullptr;
+
+    if (0==max_record_size)
+        max_record_size = 1000;
+
+    G->args = static_cast<char*>(calloc (1, 5*max_record_size));
+    if (nullptr==G->args) {
+        free (G);
+        return nullptr;
+    }
+
+    G->next_args = static_cast<char*>(calloc (1, max_record_size));
+    if (nullptr==G->args) {
+        free (G->args);
+        free (G);
+        return nullptr;
+    }
+
+    G->args_size = 5*max_record_size;
+    G->next_args_size = max_record_size;
+
+    G->tags = tags;
+    G->n_tags = n_tags;
+    return G;
+}
+
+
+
+/***************************************************************************************/
+static ffio *ffio_destroy (ffio *G) {
+/****************************************************************************************
+Free all allocated associated memory, then free G itself. For extra RAII compliancy,
+the file object should also be closed if still open, but this will require additional
+control logic, and ffio is a gie tool specific package, so we fall back to asserting that
+fclose has been called prior to ffio_destroy.
+****************************************************************************************/
+    free (G->args);
+    free (G->next_args);
+    free (G);
+    return nullptr;
+}
+
+
+
+/***************************************************************************************/
+static int at_decorative_element (ffio *G) {
+/****************************************************************************************
+A decorative element consists of a line of at least 5 consecutive identical chars,
+starting at buffer position 0:
+"-----", "=====", "*****", etc.
+
+A decorative element serves as a end delimiter for the current element, and
+continues until a gie command verb is found at the start of a line
+****************************************************************************************/
+    int i;
+    char *c;
+    if (nullptr==G)
+        return 0;
+    c = G->next_args;
+    if (nullptr==c)
+        return 0;
+    if (0==c[0])
+        return 0;
+    for (i = 1; i < 5; i++)
+        if (c[i]!=c[0])
+            return 0;
+    return 1;
+}
+
+
+
+/***************************************************************************************/
+static const char *at_tag (ffio *G) {
+/****************************************************************************************
+A start of a new command serves as an end delimiter for the current command
+****************************************************************************************/
+    size_t j;
+    for (j = 0;  j < G->n_tags;  j++)
+        if (strncmp (G->next_args, G->tags[j], strlen(G->tags[j]))==0)
+            return G->tags[j];
+    return nullptr;
+}
+
+
+
+/***************************************************************************************/
+static int at_end_delimiter (ffio *G) {
+/****************************************************************************************
+An instruction consists of everything from its introductory tag to its end
+delimiter.  An end delimiter can be either the introductory tag of the next
+instruction, or a "decorative element", i.e. one of the "ascii art" style
+block delimiters typically used to mark up block comments in a free format
+file.
+****************************************************************************************/
+    if (G==nullptr)
+        return 0;
+    if (at_decorative_element (G))
+        return 1;
+    if (at_tag (G))
+        return 1;
+    return 0;
+}
+
+
+
+/***************************************************************************************/
+static int nextline (ffio *G) {
+/****************************************************************************************
+Read next line of input file. Returns 1 on success, 0 on failure.
+****************************************************************************************/
+    G->next_args[0] = 0;
+    if (T.skip)
+        return 0;
+    if (nullptr==fgets (G->next_args, (int) G->next_args_size - 1, G->f))
+        return 0;
+    if (feof (G->f))
+        return 0;
+    pj_chomp (G->next_args);
+    G->next_lineno++;
+    return 1;
+}
+
+
+
+/***************************************************************************************/
+static int locate_tag (ffio *G, const char *tag) {
+/****************************************************************************************
+Find start-of-line tag (currently only used to search for for <gie>, but any tag
+valid).
+
+Returns 1 on success, 0 on failure.
+****************************************************************************************/
+    size_t n = strlen (tag);
+    while (0!=strncmp (tag, G->next_args, n))
+        if (0==nextline (G))
+            return 0;
+    return 1;
+}
+
+
+
+/***************************************************************************************/
+static int step_into_gie_block (ffio *G) {
+/****************************************************************************************
+Make sure we're inside a <gie>-block. Return 1 on success, 0 otherwise.
+****************************************************************************************/
+    /* Already inside */
+    if (G->level % 2)
+        return 1;
+
+    if (0==locate_tag (G, "<gie>"))
+        return 0;
+
+    while (0!=strncmp ("<gie>", G->next_args, 5)) {
+        G->next_args[0] = 0;
+        if (feof (G->f))
+            return 0;
+        if (nullptr==fgets (G->next_args, (int) G->next_args_size - 1, G->f))
+            return 0;
+        pj_chomp (G->next_args);
+        G->next_lineno++;
+    }
+    G->level++;
+
+    /* We're ready at the start - now step into the block */
+    return nextline (G);
+}
+
+
+
+/***************************************************************************************/
+static int skip_to_next_tag (ffio *G) {
+/****************************************************************************************
+Skip forward to the next command tag. Return 1 on success, 0 otherwise.
+****************************************************************************************/
+    const char *c;
+    if (0==step_into_gie_block (G))
+        return 0;
+
+    c = at_tag (G);
+
+    /* If not already there - get there */
+    while (!c) {
+        if (0==nextline (G))
+            return 0;
+        c = at_tag (G);
+    }
+
+    /* If we reached the end of a <gie> block, locate the next and retry */
+    if (0==strcmp (c, "</gie>")) {
+        G->level++;
+        if (feof (G->f))
+            return 0;
+        if (0==step_into_gie_block (G))
+            return 0;
+        G->args[0] = 0;
+        return skip_to_next_tag (G);
+    }
+    G->lineno = G->next_lineno;
+
+    return 1;
+}
+
+/* Add the most recently read line of input to the block already stored. */
+static int append_args (ffio *G) {
+    size_t skip_chars = 0;
+    size_t next_len = strlen (G->next_args);
+    size_t args_len = strlen (G->args);
+    const char *tag = at_tag (G);
+
+    if (tag)
+        skip_chars = strlen (tag);
+
+    /* +2: 1 for the space separator and 1 for the NUL termination. */
+    if (G->args_size < args_len + next_len - skip_chars + 2) {
+        char *p = static_cast<char*>(realloc (G->args, 2 * G->args_size));
+        if (nullptr==p)
+            return 0;
+        G->args = p;
+        G->args_size = 2 * G->args_size;
+    }
+
+    G->args[args_len] = ' ';
+    strcpy (G->args + args_len + 1,  G->next_args + skip_chars);
+
+    G->next_args[0] = 0;
+    return 1;
+}
+
+
+
+
+
+/***************************************************************************************/
+static int get_inp (ffio *G) {
+/****************************************************************************************
+The primary command reader for gie. Reads a block of gie input, cleans up repeated
+whitespace etc. The block is stored in G->args. Returns 1 on success, 0 otherwise.
+****************************************************************************************/
+    G->args[0] = 0;
+
+    if (0==skip_to_next_tag (G))
+        return 0;
+    G->tag = at_tag (G);
+
+    if (nullptr==G->tag)
+        return 0;
+
+    do {
+        append_args (G);
+        if (0==nextline (G))
+            return 0;
+    } while (!at_end_delimiter (G));
+
+    pj_shrink (G->args);
+    return 1;
+}
diff --git a/src/apps/nad2bin.cpp b/src/apps/nad2bin.cpp
new file mode 100644
index 00000000..ff8f2ebd
--- /dev/null
+++ b/src/apps/nad2bin.cpp
@@ -0,0 +1,382 @@
+/* Convert bivariate ASCII NAD27 to NAD83 tables to NTv2 binary structure */
+#include <stdio.h>
+#include <stdlib.h>
+
+#define PJ_LIB__
+#include "proj_internal.h"
+#include "projects.h"
+#define U_SEC_TO_RAD 4.848136811095359935899141023e-12
+
+/************************************************************************/
+/*                             swap_words()                             */
+/*                                                                      */
+/*      Convert the byte order of the given word(s) in place.           */
+/************************************************************************/
+
+static const int  byte_order_test = 1;
+#define IS_LSB	(((const unsigned char *) (&byte_order_test))[0] == 1)
+
+static void swap_words( void *data_in, int word_size, int word_count )
+
+{
+    int	word;
+    unsigned char *data = (unsigned char *) data_in;
+
+    for( word = 0; word < word_count; word++ )
+    {
+        int	i;
+        
+        for( i = 0; i < word_size/2; i++ )
+        {
+            unsigned char	t;
+            
+            t = data[i];
+            data[i] = data[word_size-i-1];
+            data[word_size-i-1] = t;
+        }
+        
+        data += word_size;
+    }
+}
+
+/************************************************************************/
+/*                               Usage()                                */
+/************************************************************************/
+
+static void Usage()
+{
+    fprintf(stderr,
+            "usage: nad2bin [-f ctable/ctable2/ntv2] binary_output < ascii_source\n" );
+    exit(1);
+}
+
+/************************************************************************/
+/*                                main()                                */
+/************************************************************************/
+int main(int argc, char **argv) {
+    struct CTABLE ct;
+    FLP *p, t;
+    size_t tsize;
+    int i, j, ichk;
+    long lam, laml, phi, phil;
+    FILE *fp;
+
+    const char *output_file = nullptr;
+
+    const char *format   = "ctable2";
+    const char *GS_TYPE  = "SECONDS";
+    const char *VERSION  = "";
+    const char *SYSTEM_F = "NAD27";
+    const char *SYSTEM_T = "NAD83";
+    const char *SUB_NAME = "";
+    const char *CREATED  = "";
+    const char *UPDATED  = "";
+
+/* ==================================================================== */
+/*      Process arguments.                                              */
+/* ==================================================================== */
+    for( i = 1; i < argc; i++ )
+    {
+        if( i < argc-1 && strcmp(argv[i],"-f") == 0 ) 
+        {
+            format = argv[++i];
+        }
+        else if( output_file == nullptr )
+        {
+            output_file = argv[i];
+        }
+        else
+            Usage();
+    }
+
+    if( output_file == nullptr )
+        Usage();
+
+    fprintf( stdout, "Output Binary File Format: %s\n", format );
+
+/* ==================================================================== */
+/*      Read the ASCII Table                                            */
+/* ==================================================================== */
+
+    memset(ct.id,0,MAX_TAB_ID);
+    if ( nullptr == fgets(ct.id, MAX_TAB_ID, stdin) ) {
+        perror("fgets");
+        exit(1);
+    }
+    /* cppcheck-suppress invalidscanf */
+    if ( EOF == scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
+          &ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi) ) {
+        perror("scanf");
+        exit(1);
+    }
+    if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi *
+                                 sizeof(FLP)))) {
+        perror("mem. alloc");
+        exit(1);
+    }
+    ct.ll.lam *= DEG_TO_RAD;
+    ct.ll.phi *= DEG_TO_RAD;
+    ct.del.lam *= DEG_TO_RAD;
+    ct.del.phi *= DEG_TO_RAD;
+    /* load table */
+    p = ct.cvs;
+    for (i = 0; i < ct.lim.phi; ++i) {
+        /* cppcheck-suppress invalidscanf */
+        if ( EOF == scanf("%d:%ld %ld", &ichk, &laml, &phil) ) {
+            perror("scanf on row");
+            exit(1);
+        }
+        if (ichk != i) {
+            fprintf(stderr,"format check on row\n");
+            exit(1);
+        }
+        t.lam = (float) (laml * U_SEC_TO_RAD);
+        t.phi = (float) (phil * U_SEC_TO_RAD);
+        *p++ = t;
+        for (j = 1; j < ct.lim.lam; ++j) {
+            /* cppcheck-suppress invalidscanf */
+            if ( EOF == scanf("%ld %ld", &lam, &phi) ) {
+                perror("scanf on column");
+                exit(1);
+            }
+            t.lam = (float) ((laml += lam) * U_SEC_TO_RAD);
+            t.phi = (float) ((phil += phi) * U_SEC_TO_RAD);
+            *p++ = t;
+        }
+    }
+    if (feof(stdin)) {
+        fprintf(stderr, "premature EOF\n");
+        exit(1);
+    }
+
+/* ==================================================================== */
+/*      Write out the old ctable format - this is machine and byte      */
+/*      order specific.                                                 */
+/* ==================================================================== */
+    if( strcmp(format,"ctable") == 0 ) 
+    {
+	if (!(fp = fopen(output_file, "wb"))) {
+            perror(output_file);
+            exit(2);
+	}
+	if (fwrite(&ct, sizeof(ct), 1, fp) != 1 ||
+            fwrite(ct.cvs, tsize, 1, fp) != 1) {
+            fprintf(stderr, "output failure\n");
+            exit(2);
+	}
+        fclose( fp );
+	exit(0); /* normal completion */
+    }
+
+/* ==================================================================== */
+/*      Write out the old ctable format - this is machine and byte      */
+/*      order specific.                                                 */
+/* ==================================================================== */
+    if( strcmp(format,"ctable2") == 0 ) 
+    {
+        char header[160];
+
+	if (!(fp = fopen(output_file, "wb"))) {
+            perror(output_file);
+            exit(2);
+	}
+
+	/* cppcheck-suppress sizeofCalculation */
+        STATIC_ASSERT( MAX_TAB_ID == 80 );
+        /* cppcheck-suppress sizeofCalculation */
+        STATIC_ASSERT( sizeof(pj_int32) == 4 ); /* for ct.lim.lam/phi */
+
+        memset( header, 0, sizeof(header) );
+
+        memcpy( header +   0, "CTABLE V2.0     ", 16 );
+        memcpy( header +  16, ct.id, 80 );
+        memcpy( header +  96, &ct.ll.lam, 8 );
+        memcpy( header + 104, &ct.ll.phi, 8 );
+        memcpy( header + 112, &ct.del.lam, 8 );
+        memcpy( header + 120, &ct.del.phi, 8 );
+        memcpy( header + 128, &ct.lim.lam, 4 );
+        memcpy( header + 132, &ct.lim.phi, 4 );
+
+        /* force into LSB format */
+        if( !IS_LSB ) 
+        {
+            swap_words( header +  96, 8, 4 );
+            swap_words( header + 128, 4, 2 );
+            swap_words( ct.cvs, 4, ct.lim.lam * 2 * ct.lim.phi );
+        }
+
+        if( fwrite( header, sizeof(header), 1, fp ) != 1 ) {
+            perror( "fwrite" );
+            exit( 2 );
+        }
+
+	if (fwrite(ct.cvs, tsize, 1, fp) != 1) {
+            perror( "fwrite" );
+            exit(2);
+	}
+
+        fclose( fp );
+	exit(0); /* normal completion */
+    }
+
+/* ==================================================================== */
+/*      Write out the NTv2 format grid shift file.                      */
+/* ==================================================================== */
+    if( strcmp(format,"ntv2") == 0 ) 
+    {
+        if (!(fp = fopen(output_file, "wb"))) 
+        {
+            perror(output_file);
+            exit(2);
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Write the file header.                                          */
+/* -------------------------------------------------------------------- */
+        {    
+            char achHeader[11*16];
+
+            memset( achHeader, 0, sizeof(achHeader) );
+        
+            memcpy( achHeader +  0*16, "NUM_OREC", 8 );
+            achHeader[ 0*16 + 8] = 0xb;
+
+            memcpy( achHeader +  1*16, "NUM_SREC", 8 );
+            achHeader[ 1*16 + 8] = 0xb;
+
+            memcpy( achHeader +  2*16, "NUM_FILE", 8 );
+            achHeader[ 2*16 + 8] = 0x1;
+
+            memcpy( achHeader +  3*16, "GS_TYPE         ", 16 );
+            memcpy( achHeader +  3*16+8, GS_TYPE, MIN(16,strlen(GS_TYPE)) );
+
+            memcpy( achHeader +  4*16, "VERSION         ", 16 );
+            memcpy( achHeader +  4*16+8, VERSION, MIN(16,strlen(VERSION)) );
+
+            memcpy( achHeader +  5*16, "SYSTEM_F        ", 16 );
+            memcpy( achHeader +  5*16+8, SYSTEM_F, MIN(16,strlen(SYSTEM_F)) );
+
+            memcpy( achHeader +  6*16, "SYSTEM_T        ", 16 );
+            memcpy( achHeader +  6*16+8, SYSTEM_T, MIN(16,strlen(SYSTEM_T)) );
+
+            memcpy( achHeader +  7*16, "MAJOR_F ", 8);
+            memcpy( achHeader +  8*16, "MINOR_F ", 8 );
+            memcpy( achHeader +  9*16, "MAJOR_T ", 8 );
+            memcpy( achHeader + 10*16, "MINOR_T ", 8 );
+
+            fwrite( achHeader, 1, sizeof(achHeader), fp );
+        }
+        
+/* -------------------------------------------------------------------- */
+/*      Write the grid header.                                          */
+/* -------------------------------------------------------------------- */
+        {
+            unsigned char achHeader[11*16];
+            double dfValue;
+            pj_int32 nGSCount = ct.lim.lam * ct.lim.phi;
+            LP ur;
+
+            ur.lam = ct.ll.lam + (ct.lim.lam-1) * ct.del.lam;
+            ur.phi = ct.ll.phi + (ct.lim.phi-1) * ct.del.phi;
+
+            /* cppcheck-suppress sizeofCalculation */
+            STATIC_ASSERT( sizeof(nGSCount) == 4 );
+
+            memset( achHeader, 0, sizeof(achHeader) );
+
+            memcpy( achHeader +  0*16, "SUB_NAME        ", 16 );
+            memcpy( achHeader +  0*16+8, SUB_NAME, MIN(16,strlen(SUB_NAME)) );
+    
+            memcpy( achHeader +  1*16, "PARENT          ", 16 );
+            memcpy( achHeader +  1*16+8, "NONE", MIN(16,strlen("NONE")) );
+    
+            memcpy( achHeader +  2*16, "CREATED         ", 16 );
+            memcpy( achHeader +  2*16+8, CREATED, MIN(16,strlen(CREATED)) );
+    
+            memcpy( achHeader +  3*16, "UPDATED         ", 16 );
+            memcpy( achHeader +  3*16+8, UPDATED, MIN(16,strlen(UPDATED)) );
+
+            memcpy( achHeader +  4*16, "S_LAT   ", 8 );
+            dfValue = ct.ll.phi * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  4*16 + 8, &dfValue, 8 );
+
+            memcpy( achHeader +  5*16, "N_LAT   ", 8 );
+            dfValue = ur.phi * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  5*16 + 8, &dfValue, 8 );
+
+            memcpy( achHeader +  6*16, "E_LONG  ", 8 );
+            dfValue = -1 * ur.lam * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  6*16 + 8, &dfValue, 8 );
+
+            memcpy( achHeader +  7*16, "W_LONG  ", 8 );
+            dfValue = -1 * ct.ll.lam * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  7*16 + 8, &dfValue, 8 );
+
+            memcpy( achHeader +  8*16, "LAT_INC ", 8 );
+            dfValue = ct.del.phi * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  8*16 + 8, &dfValue, 8 );
+    
+            memcpy( achHeader +  9*16, "LONG_INC", 8 );
+            dfValue = ct.del.lam * 3600.0 / DEG_TO_RAD;
+            memcpy( achHeader +  9*16 + 8, &dfValue, 8 );
+    
+            memcpy( achHeader + 10*16, "GS_COUNT", 8 );
+            memcpy( achHeader + 10*16+8, &nGSCount, 4 );
+    
+            if( !IS_LSB ) 
+            {
+                swap_words( achHeader +  4*16 + 8, 8, 1 );
+                swap_words( achHeader +  5*16 + 8, 8, 1 );
+                swap_words( achHeader +  6*16 + 8, 8, 1 );
+                swap_words( achHeader +  7*16 + 8, 8, 1 );
+                swap_words( achHeader +  8*16 + 8, 8, 1 );
+                swap_words( achHeader +  9*16 + 8, 8, 1 );
+                swap_words( achHeader + 10*16 + 8, 4, 1 );
+            }
+
+            fwrite( achHeader, 1, sizeof(achHeader), fp );
+        }
+
+/* -------------------------------------------------------------------- */
+/*      Write the actual grid cells.                                    */
+/* -------------------------------------------------------------------- */
+        {
+            float *row_buf;
+            int row;
+
+            row_buf = (float *) pj_malloc(ct.lim.lam * sizeof(float) * 4);
+            memset( row_buf, 0, sizeof(float)*4 );
+
+            for( row = 0; row < ct.lim.phi; row++ )
+            {
+                for( i = 0; i < ct.lim.lam; i++ )
+                {
+                    FLP *cvs = ct.cvs + (row) * ct.lim.lam
+                        + (ct.lim.lam - i - 1);
+
+                    /* convert radians to seconds */
+                    row_buf[i*4+0] = (float) (cvs->phi * (3600.0 / (M_PI/180.0)));
+                    row_buf[i*4+1] = (float) (cvs->lam * (3600.0 / (M_PI/180.0)));
+
+                    /* We leave the accuracy values as zero */
+                }
+
+                if( !IS_LSB )
+                    swap_words( row_buf, 4, ct.lim.lam * 4 );
+
+                if( fwrite( row_buf, sizeof(float), ct.lim.lam*4, fp ) 
+                    != (size_t)( 4 * ct.lim.lam ) )
+                {
+                    perror( "write()" );
+                    exit( 2 );
+                }
+            }
+        }
+
+        fclose( fp );
+        exit(0); /* normal completion */
+    }
+
+    fprintf( stderr, "Unsupported format, nothing written.\n" );
+    exit( 3 );
+}
diff --git a/src/apps/optargpm.h b/src/apps/optargpm.h
new file mode 100644
index 00000000..035c6f92
--- /dev/null
+++ b/src/apps/optargpm.h
@@ -0,0 +1,635 @@
+/***********************************************************************
+
+          OPTARGPM - a header-only library for decoding
+                PROJ.4 style command line options
+
+                   Thomas Knudsen, 2017-09-10
+
+************************************************************************
+
+For PROJ.4 command line programs, we have a somewhat complex option
+decoding situation, since we have to navigate in a cocktail of classic
+single letter style options, prefixed by "-", GNU style long options
+prefixed by "--", transformation specification elements prefixed by "+",
+and input file names prefixed by "" (i.e. nothing).
+
+Hence, classic getopt.h style decoding does not cut the mustard, so
+this is an attempt to catch up and chop the ketchup.
+
+Since optargpm (for "optarg plus minus") does not belong, in any
+obvious way, in any systems development library, it is provided as
+a "header only" library.
+
+While this is conventional in C++, it is frowned at in plain C.
+But frown away - "header only" has its places, and this is one of
+them.
+
+By convention, we expect a command line to consist of the following
+elements:
+
+        <operator/program name>
+        [short ("-")/long ("--") options}
+        [operator ("+") specs]
+        [operands/input files]
+
+or less verbose:
+
+        <operator>   [options]   [operator specs]   [operands]
+
+or less abstract:
+
+   proj  -I --output=foo  +proj=utm +zone=32 +ellps=GRS80   bar baz...
+
+Where
+
+Operator is              proj
+Options are             -I --output=foo
+Operator specs are      +proj=utm +zone=32 +ellps=GRS80
+Operands are             bar baz
+
+
+While neither claiming to save the world, nor to hint at the "shape of
+jazz to come", at least optargpm has shown useful in constructing cs2cs
+style transformation filters.
+
+Supporting a wide range of option syntax, the getoptpm API is somewhat
+quirky, but also compact, consisting of one data type, 3(+2) functions,
+and one enumeration:
+
+OPTARGS
+        Housekeeping data type. An instance of OPTARGS is conventionally
+        called o or opt
+opt_parse (opt, argc, argv ...):
+        The work horse: Define supported options; Split (argc, argv)
+        into groups (options, op specs, operands); Parse option
+        arguments.
+opt_given (o, option):
+        The number of times <option> was given on the command line.
+        (i.e. 0 if not given or option unsupported)
+opt_arg (o, option):
+        A char pointer to the argument for <option>
+
+The 2 additional functions (of which, one is really a macro) implements
+a "read all operands sequentially" functionality, eliminating the need to
+handle open/close of a sequence of input files:
+
+enum OPTARGS_FILE_MODE:
+        indicates whether to read operands in text (0) or binary (1) mode
+opt_input_loop (o, mode):
+        When used as condition in a while loop, traverses all operands,
+        giving the impression of reading just a single input file.
+opt_eof_handler (o):
+        Auxiliary macro, to be called inside the input loop after each
+        read operation
+
+Usage is probably easiest understood by a brief textbook style example:
+
+Consider a simple program taking the conventional "-v, -h, -o" options
+indicating "verbose output", "help please", and "output file specification",
+respectively.
+
+The "-v" and "-h" options are *flags*, taking no arguments, while the
+"-o" option is a *key*, taking a *value* argument, representing the
+output file name.
+
+The short options have long aliases: "--verbose", "--help" and "--output".
+Additionally, the long key "--hello", without any short counterpart, is
+supported.
+
+-------------------------------------------------------------------------------
+
+
+int main(int argc, char **argv) {
+    PJ *P;
+    OPTARGS *o;
+    FILE *out = stdout;
+    char *longflags[]  = {"v=verbose", "h=help", 0};
+    char *longkeys[]   = {"o=output", "hello", 0};
+
+    o = opt_parse (argc, argv, "hv", "o", longflags, longkeys);
+    if (nullptr==o)
+        return nullptr;
+
+
+    if (opt_given (o, "h")) {
+        printf ("Usage: %s [-v|--verbose] [-h|--help] [-o|--output <filename>] [--hello=<name>] infile...", o->progname);
+        exit (0);
+    }
+
+    if (opt_given (o, "v"))
+        puts ("Feeling chatty today?");
+
+    if (opt_given (o, "hello")) {
+        printf ("Hello, %s!\n", opt_arg(o, "hello"));
+        exit (0);
+    }
+
+    if (opt_given (o, "o"))
+        out = fopen (opt_arg (o, "output"), "rt"); // Note: "output" translates to "o" internally
+
+    // Setup transformation
+    P = proj_create_argv (0, o->pargc, o->pargv);
+
+    // Loop over all lines of all input files
+    while (opt_input_loop (o, optargs_file_format_text)) {
+        char buf[1000];
+        int ret = fgets (buf, 1000, o->input);
+        opt_eof_handler (o);
+        if (nullptr==ret) {
+            fprintf (stderr, "Read error in record %d\n", (int) o->record_index);
+            continue;
+        }
+        do_what_needs_to_be_done (buf);
+    }
+
+    return nullptr;
+}
+
+
+-------------------------------------------------------------------------------
+
+Note how short aliases for longflags and longkeys are defined by prefixing
+an "o=", "h=" or "v=", respectively. This also means that it is possible to
+have more than one alias for each short option, e.g.
+
+               longkeys = {"o=output", "o=banana", 0}
+
+would define both "--output" and "--banana" to be aliases for "-o".
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-10
+
+************************************************************************
+
+* Copyright (c) 2016, 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.
+
+***********************************************************************/
+#include <ctype.h>
+#include <errno.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+/**************************************************************************************************/
+struct OPTARGS;
+typedef struct OPTARGS OPTARGS;
+enum OPTARGS_FILE_FORMAT {optargs_file_format_text = 0, optargs_file_format_binary = 1};
+
+char *opt_filename (OPTARGS *opt);
+static int opt_eof (OPTARGS *opt);
+int opt_record (OPTARGS *opt);
+int opt_input_loop (OPTARGS *opt, int binary);
+static int opt_is_flag (OPTARGS *opt, int ordinal);
+static int opt_raise_flag (OPTARGS *opt, int ordinal);
+static int opt_ordinal (OPTARGS *opt, const char *option);
+int opt_given (OPTARGS *opt, const char *option);
+char *opt_arg (OPTARGS *opt, const char *option);
+const char *opt_strip_path (const char *full_name);
+OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys, const char **longflags, const char **longkeys);
+
+#define opt_eof_handler(opt) if (opt_eof (opt)) {continue;} else {;}
+/**************************************************************************************************/
+
+struct OPTARGS {
+    int    argc,   margc,   pargc,   fargc;
+    char **argv, **margv, **pargv, **fargv;
+    FILE *input;
+    int   input_index;
+    int   record_index;
+    const char  *progname;   /* argv[0], stripped from /path/to, if present */
+    char   flaglevel[21];    /* if flag -f is specified n times, its optarg pointer is set to flaglevel + n */
+    char  *optarg[256];      /* optarg[(int) 'f'] holds a pointer to the argument of option "-f" */
+    char  *flags;            /* a list of flag style options supported, e.g. "hv" (help and verbose) */
+    char  *keys;             /* a list of key/value style options supported, e.g. "o" (output) */
+    const char **longflags;  /* long flags, {"help", "verbose"}, or {"h=help", "v=verbose"}, to indicate homologous short options */
+    const char **longkeys;   /* e.g. {"output"} or {o=output"} to support --output=/path/to/output-file. In the latter case, */
+                             /* all operations on "--output" gets redirected to "-o", so user code need handle arguments to "-o" only */
+};
+
+
+/* name of file currently read from */
+char *opt_filename (OPTARGS *opt) {
+    if (nullptr==opt)
+        return nullptr;
+    if (0==opt->fargc)
+        return opt->flaglevel;
+    return opt->fargv[opt->input_index];
+}
+
+static int opt_eof (OPTARGS *opt) {
+    if (nullptr==opt)
+        return 1;
+    return feof (opt->input);
+}
+
+/* record number of most recently read record */
+int opt_record (OPTARGS *opt) {
+    if (nullptr==opt)
+        return 0;
+    return opt->record_index + 1;
+}
+
+
+/* handle closing/opening of a "stream-of-streams" */
+int opt_input_loop (OPTARGS *opt, int binary) {
+    if (nullptr==opt)
+        return 0;
+
+    /* most common case: increment record index and read on */
+    if ( (opt->input!=nullptr) && !feof (opt->input) ) {
+        opt->record_index++;
+        return 1;
+    }
+
+    opt->record_index = 0;
+
+    /* no input files specified - read from stdin */
+    if ((0==opt->fargc) && (nullptr==opt->input)) {
+        opt->input = stdin;
+        return 1;
+    }
+
+    /* if we're here, we have either reached eof on current input file.   */
+    /*  or not yet opened a file. If eof on stdin, we're done             */
+    if (opt->input==stdin)
+        return 0;
+
+    /* end if no more input */
+    if (nullptr!=opt->input)
+        fclose (opt->input);
+    if (opt->input_index >= opt->fargc)
+        return 0;
+
+    /* otherwise, open next input file */
+    opt->input = fopen (opt->fargv[opt->input_index++], binary? "rb": "rt");
+    if (nullptr != opt->input)
+        return 1;
+
+    /* ignore non-existing files - go on! */
+    if (nullptr==opt->input)
+        return opt_input_loop (opt, binary);
+    return 0;
+}
+
+
+/* return true if option with given ordinal is a flag, false if undefined or key=value */
+static int opt_is_flag (OPTARGS *opt, int ordinal) {
+    if (opt->optarg[ordinal] < opt->flaglevel)
+        return 0;
+    if (opt->optarg[ordinal] > opt->flaglevel + 20)
+        return 0;
+    return 1;
+}
+
+static int opt_raise_flag (OPTARGS *opt, int ordinal) {
+    if (opt->optarg[ordinal] < opt->flaglevel)
+        return 1;
+    if (opt->optarg[ordinal] > opt->flaglevel + 20)
+        return 1;
+
+    /* Max out at 20 */
+    if (opt->optarg[ordinal]==opt->flaglevel + 20)
+        return 0;
+    opt->optarg[ordinal]++;
+    return 0;
+}
+
+/* Find the ordinal value of any (short or long) option */
+static int opt_ordinal (OPTARGS *opt, const char *option) {
+    int i;
+    if (nullptr==opt)
+        return 0;
+    if (nullptr==option)
+        return 0;
+    if (0==option[0])
+        return 0;
+    /* An ordinary -o style short option */
+    if (strlen (option)==1) {
+        /* Undefined option? */
+        if (nullptr==opt->optarg[(int) option[0]])
+            return 0;
+        return (int) option[0];
+    }
+
+    /* --longname style long options are slightly harder */
+    for (i = 0; i < 64; i++) {
+        const char **f = opt->longflags;
+        if (nullptr==f)
+            break;
+        if (nullptr==f[i])
+            break;
+        if (0==strcmp(f[i], "END"))
+            break;
+        if (0==strcmp(f[i], option))
+            return 128 + i;
+
+        /* long alias? - return ordinal for corresponding short */
+        if ((strlen(f[i]) > 2) && (f[i][1]=='=') && (0==strcmp(f[i]+2, option))) {
+                /* Undefined option? */
+                if (nullptr==opt->optarg[(int) f[i][0]])
+                    return 0;
+                return (int) f[i][0];
+        }
+    }
+
+    for (i = 0; i < 64; i++) {
+        const char **v = opt->longkeys;
+        if (nullptr==v)
+            return 0;
+        if (nullptr==v[i])
+            return 0;
+        if (0==strcmp (v[i], "END"))
+            return 0;
+        if (0==strcmp(v[i], option))
+            return 192 + i;
+
+        /* long alias? - return ordinal for corresponding short */
+        if ((strlen(v[i]) > 2) && (v[i][1]=='=') && (0==strcmp(v[i]+2, option))) {
+            /* Undefined option? */
+            if (nullptr==opt->optarg[(int) v[i][0]])
+                return 0;
+            return (int) v[i][0];
+        }
+
+    }
+    /* kill some potential compiler warnings about unused functions */
+    (void) opt_eof (nullptr);
+    return 0;
+}
+
+
+/* Returns 0 if option was not given on command line, non-0 otherwise */
+int opt_given (OPTARGS *opt, const char *option) {
+    int ordinal = opt_ordinal (opt, option);
+    if (0==ordinal)
+        return 0;
+    /* For flags we return the number of times the flag was specified (mostly for repeated -v(erbose) flags) */
+    if (opt_is_flag (opt, ordinal))
+        return (int) (opt->optarg[ordinal] - opt->flaglevel);
+    return opt->argv[0] != opt->optarg[ordinal];
+}
+
+
+/* Returns the argument to a given option */
+char *opt_arg (OPTARGS *opt, const char *option) {
+    int ordinal = opt_ordinal (opt, option);
+    if (0==ordinal)
+        return nullptr;
+    return opt->optarg[ordinal];
+}
+
+const char *opt_strip_path (const char *full_name) {
+    const char *last_path_delim, *stripped_name = full_name;
+
+    last_path_delim = strrchr (stripped_name, '\\');
+    if (last_path_delim > stripped_name)
+        stripped_name = last_path_delim + 1;
+
+    last_path_delim = strrchr (stripped_name, '/');
+    if (last_path_delim > stripped_name)
+        stripped_name = last_path_delim + 1;
+    return stripped_name;
+}
+
+/* split command line options into options/flags ("-" style), projdefs ("+" style) and input file args */
+OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys, const char **longflags, const char **longkeys) {
+    int i, j;
+    int free_format;
+    OPTARGS *o;
+
+    o = (OPTARGS *) calloc (1, sizeof(OPTARGS));
+    if (nullptr==o)
+        return nullptr;
+
+    o->argc = argc;
+    o->argv = argv;
+    o->progname = opt_strip_path (argv[0]);
+
+    /* Reset all flags */
+    for (i = 0;  i < (int) strlen (flags);  i++)
+        o->optarg[(int) flags[i]] = o->flaglevel;
+
+     /* Flag args for all argument taking options as "unset" */
+    for (i = 0; i < (int) strlen (keys); i++)
+        o->optarg[(int) keys[i]] = argv[0];
+
+    /* Hence, undefined/illegal options have an argument of 0 */
+
+    /* long opts are handled similarly, but are mapped to the high bit character range (above 128) */
+    o->longflags  =  longflags;
+    o->longkeys   =  longkeys;
+
+
+    /* check aliases, An end user should never experience this, but */
+    /* the developer should make sure that aliases are valid */
+    for (i = 0;  longflags && longflags[i]; i++) {
+        /* Go on if it does not look like an alias */
+        if (strlen (longflags[i]) < 3)
+            continue;
+        if ('='!=longflags[i][1])
+            continue;
+        if (nullptr==strchr (flags, longflags[i][0])) {
+            fprintf (stderr, "%s: Invalid alias - '%s'. Valid short flags are '%s'\n", o->progname, longflags[i], flags);
+            free (o);
+            return nullptr;
+        }
+    }
+    for (i = 0;  longkeys && longkeys[i]; i++) {
+        /* Go on if it does not look like an alias */
+        if (strlen (longkeys[i]) < 3)
+            continue;
+        if ('='!=longkeys[i][1])
+            continue;
+        if (nullptr==strchr (keys, longkeys[i][0])) {
+            fprintf (stderr, "%s: Invalid alias - '%s'. Valid short flags are '%s'\n", o->progname, longkeys[i], keys);
+            free (o);
+            return nullptr;
+        }
+    }
+
+    /* aside from counting the number of times a flag has been specified, we also abuse the */
+    /* flaglevel array to provide a pseudo-filename for the case of reading from stdin      */
+    strcpy (o->flaglevel, "<stdin>");
+
+    for (i = 128; (longflags != nullptr) && (longflags[i - 128] != nullptr); i++) {
+        if (i==192) {
+            free (o);
+            fprintf (stderr, "Too many flag style long options\n");
+            return nullptr;
+        }
+        o->optarg[i] = o->flaglevel;
+    }
+
+    for (i = 192;  (longkeys != nullptr) && (longkeys[i - 192] != nullptr);  i++) {
+        if (i==256) {
+            free (o);
+            fprintf (stderr, "Too many value style long options\n");
+            return nullptr;
+        }
+        o->optarg[i] = argv[0];
+    }
+
+    /* Now, set up the agrc/argv pairs, and interpret args */
+    o->argc = argc;
+    o->argv = argv;
+
+    /* Process all '-' and '--'-style options */
+    for (i = 1;  i < argc;  i++) {
+        int arg_group_size = (int) strlen (argv[i]);
+
+        if ('-' != argv[i][0])
+            break;
+
+        if (nullptr==o->margv)
+            o->margv = argv + i;
+        o->margc++;
+
+        for (j = 1;  j < arg_group_size;  j++) {
+            int c =  argv[i][j];
+            char cstring[2], *crepr = cstring;
+            cstring[0] = (char) c;
+            cstring[1] = 0;
+
+
+            /* Long style flags and options (--long_opt_name, --long_opt_namr arg, --long_opt_name=arg) */
+            if (c== (int)'-') {
+                char *equals;
+                crepr = argv[i] + 2;
+
+                /* We need to manipulate a bit to support gnu style --foo=bar syntax.   */
+                /* NOTE: This will segfault for read-only (const char * style) storage, */
+                /* but since the canonical use case, int main (int argc, char **argv),  */
+                /* is non-const, we ignore this for now */
+                equals = strchr (crepr, '=');
+                if (equals)
+                    *equals = 0;
+                c = opt_ordinal (o, crepr);
+                if (0==c) {
+                    fprintf (stderr, "Invalid option \"%s\"\n", crepr);
+                    return nullptr;
+                }
+
+                /* inline (gnu) --foo=bar style arg */
+                if (equals) {
+                    *equals = '=';
+                    if (opt_is_flag (o, c)) {
+                        fprintf (stderr, "Option \"%s\" takes no arguments\n", crepr);
+                        return nullptr;
+                    }
+                    o->optarg[c] = equals + 1;
+                    break;
+                }
+
+                /* "outline" --foo bar style arg */
+                if (!opt_is_flag (o, c)) {
+                    if ((argc==i + 1) || ('+'==argv[i+1][0]) || ('-'==argv[i+1][0])) {
+                        fprintf (stderr, "Missing argument for option \"%s\"\n", crepr);
+                        return nullptr;
+                    }
+                    o->optarg[c] = argv[i + 1];
+                    i++; /* eat the arg */
+                    break;
+                }
+
+                if (!opt_is_flag (o, c)) {
+                    fprintf (stderr, "Expected flag style long option here, but got \"%s\"\n", crepr);
+                    return nullptr;
+                }
+
+                /* Flag style option, i.e. taking no arguments */
+                opt_raise_flag (o, c);
+                break;
+            }
+
+            /* classic short options */
+            if (nullptr==o->optarg[c]) {
+                fprintf (stderr, "Invalid option \"%s\"\n", crepr);
+                return nullptr;
+            }
+
+            /* Flag style option, i.e. taking no arguments */
+            if (opt_is_flag (o, c)) {
+                opt_raise_flag (o, c);
+                continue;
+            }
+
+            /* options taking argumants */
+
+            /* argument separate (i.e. "-i 10") */
+            if (j + 1==arg_group_size) {
+                if ((argc==i + 1) || ('+'==argv[i+1][0]) || ('-'==argv[i+1][0]))
+                {
+                    fprintf (stderr, "Bad or missing arg for option \"%s\"\n", crepr);
+                    return nullptr;
+                }
+                o->optarg[(int) c] = argv[i + 1];
+                i++;
+                break;
+            }
+
+            /* Option arg inline (i.e. "-i10") */
+            o->optarg[c] = argv[i] + j + 1;
+            break;
+        }
+    }
+
+    /* Process all '+'-style options, starting from where '-'-style processing ended */
+    o->pargv = argv + i;
+
+    /* Is free format in use, instead of plus-style? */
+    free_format = 0;
+    for (j = 1;  j < argc;  j++) {
+        if (0==strcmp ("--", argv[j])) {
+            free_format = j;
+            break;
+        }
+    }
+
+    if (free_format) {
+        o->pargc = free_format - (o->margc + 1);
+        o->fargc = argc - (free_format + 1);
+        if (0 != o->fargc)
+            o->fargv = argv + free_format + 1;
+        return o;
+    }
+
+    for (/* empty */; i < argc; i++) {
+        if ('-' == argv[i][0]) {
+            free (o);
+            fprintf (stderr, "+ and - style options must not be mixed\n");
+            return nullptr;
+        }
+
+        if ('+' != argv[i][0])
+            break;
+        o->pargc++;
+    }
+
+    /* Handle input file names */
+    o->fargc = argc - i;
+    if (0!=o->fargc)
+        o->fargv = argv + i;
+
+    return o;
+
+}
diff --git a/src/apps/p_series.cpp b/src/apps/p_series.cpp
new file mode 100644
index 00000000..cddea888
--- /dev/null
+++ b/src/apps/p_series.cpp
@@ -0,0 +1,42 @@
+/* print row coefficients of Tseries structure */
+#include "projects.h"
+#include <stdio.h>
+#include <string.h>
+#define NF 20 /* length of final format string */
+#define CUT 60 /* check length of line */
+
+/* FIXME: put the declaration in a header. Also used in gen_cheb.c */
+void p_series(Tseries *T, FILE *file, char *fmt);
+
+void p_series(Tseries *T, FILE *file, char *fmt) {
+	int i, j, n, L;
+	char format[NF+1];
+
+	*format = ' ';
+	strncpy(format + 1, fmt, NF - 3);
+	strcat(format, "%n");
+	fprintf(file, "u: %d\n", T->mu+1);
+	for (i = 0; i <= T->mu; ++i)
+		if (T->cu[i].m) {
+			fprintf(file, "%d %d%n", i, T->cu[i].m, &L);
+			n = 0;
+			for (j = 0; j < T->cu[i].m; ++j) {
+				if ((L += n) > CUT)
+					fprintf(file, "\n %n", &L);
+				fprintf(file, format, T->cu[i].c[j], &n);
+			}
+			fputc('\n', file);
+		}
+	fprintf(file, "v: %d\n", T->mv+1);
+	for (i = 0; i <= T->mv; ++i)
+		if (T->cv[i].m) {
+			fprintf(file, "%d %d%n", i, T->cv[i].m, &L);
+			n = 0;
+			for (j = 0; j < T->cv[i].m; ++j) {
+				if ((L += n) > 60)
+					fprintf(file, "\n %n", &L);
+				fprintf(file, format, T->cv[i].c[j], &n);
+			}
+			fputc('\n', file);
+		}
+}
diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp
new file mode 100644
index 00000000..b93fb04d
--- /dev/null
+++ b/src/apps/proj.cpp
@@ -0,0 +1,582 @@
+/* <<<< Cartographic projection filter program >>>> */
+#include "proj.h"
+#include "projects.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <math.h>
+#include "emess.h"
+
+#if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__WIN32__)
+#  include <fcntl.h>
+#  include <io.h>
+#  define SET_BINARY_MODE(file) _setmode(_fileno(file), O_BINARY)
+#else
+#  define SET_BINARY_MODE(file)
+#endif
+
+#define MAX_LINE 1000
+#define MAX_PARGS 100
+#define PJ_INVERS(P) (P->inv ? 1 : 0)
+
+extern void gen_cheb(int, projUV(*)(projUV), const char *, PJ *, int, char **);
+
+static PJ *Proj;
+static union {
+    projUV (*generic)(projUV, PJ *);
+    projXY (*fwd)(projLP, PJ *);
+    projLP (*inv)(projXY, PJ *);
+} proj;
+
+static int
+    reversein = 0,  /* != 0 reverse input arguments */
+    reverseout = 0, /* != 0 reverse output arguments */
+    bin_in = 0,     /* != 0 then binary input */
+    bin_out = 0,    /* != 0 then binary output */
+    echoin = 0,     /* echo input data to output line */
+    tag = '#',      /* beginning of line tag character */
+    inverse = 0,    /* != 0 then inverse projection */
+    prescale = 0,   /* != 0 apply cartesian scale factor */
+    dofactors = 0,  /* determine scale factors */
+    very_verby = 0, /* very verbose mode */
+    postscale = 0;
+
+static const char
+    *cheby_str,         /* string controlling Chebychev evaluation */
+    *oform = nullptr; /* output format for x-y or decimal degrees */
+static char oform_buffer[16];   /* Buffer for oform when using -d */
+
+static const char
+    *oterr = "*\t*",    /* output line for unprojectable input */
+    *usage = "%s\nusage: %s [ -bdeEfiIlmorsStTvVwW [args] ] [ +opts[=arg] ] [ files ]\n";
+
+static PJ_FACTORS facs;
+
+static double (*informat)(const char *, char **), /* input data deformatter function */
+              fscale = 0.;                        /* cartesian scale factor */
+
+static projUV int_proj(projUV data) {
+    if (prescale) {
+        data.u *= fscale;
+        data.v *= fscale;
+    }
+
+    data = (*proj.generic)(data, Proj);
+
+    if (postscale && data.u != HUGE_VAL) {
+        data.u *= fscale;
+        data.v *= fscale;
+    }
+
+    return data;
+}
+
+/* file processing function */
+static void process(FILE *fid) {
+    char line[MAX_LINE+3], *s = nullptr, pline[40];
+    PJ_COORD data;
+
+    for (;;) {
+        int facs_bad = 0;
+        ++emess_dat.File_line;
+
+        if (bin_in) {   /* binary input */
+            if (fread(&data, sizeof(projUV), 1, fid) != 1)
+                break;
+        } else {    /* ascii input */
+            if (!(s = fgets(line, MAX_LINE, fid)))
+                break;
+
+            if (!strchr(s, '\n')) { /* overlong line */
+                int c;
+                (void)strcat(s, "\n");
+                /* gobble up to newline */
+                while ((c = fgetc(fid)) != EOF && c != '\n') ;
+            }
+
+            if (*s == tag) {
+                if (!bin_out)
+                    (void)fputs(line, stdout);
+                continue;
+            }
+
+            if (reversein) {
+                data.uv.v = (*informat)(s, &s);
+                data.uv.u = (*informat)(s, &s);
+            } else {
+                data.uv.u = (*informat)(s, &s);
+                data.uv.v = (*informat)(s, &s);
+            }
+
+            if (data.uv.v == HUGE_VAL)
+                data.uv.u = HUGE_VAL;
+
+            if (!*s && (s > line)) --s; /* assumed we gobbled \n */
+            if (!bin_out && echoin) {
+                char t;
+                t = *s;
+                *s = '\0';
+                (void)fputs(line, stdout);
+                *s = t;
+                putchar('\t');
+            }
+        }
+
+        if (data.uv.u != HUGE_VAL) {
+            PJ_COORD coord;
+            coord.lp = data.lp;
+            if (prescale) { data.uv.u *= fscale; data.uv.v *= fscale; }
+            if (dofactors && !inverse) {
+                facs = proj_factors(Proj, coord);
+                facs_bad = proj_errno(Proj);
+            }
+
+            data.xy = (*proj.fwd)(data.lp, Proj);
+
+            if (dofactors && inverse) {
+                facs = proj_factors(Proj, coord);
+                facs_bad = proj_errno(Proj);
+            }
+
+            if (postscale && data.uv.u != HUGE_VAL)
+                { data.uv.u *= fscale; data.uv.v *= fscale; }
+        }
+
+        if (bin_out) { /* binary output */
+            (void)fwrite(&data, sizeof(projUV), 1, stdout);
+            continue;
+        } else if (data.uv.u == HUGE_VAL) /* error output */
+            (void)fputs(oterr, stdout);
+        else if (inverse && !oform) {   /*ascii DMS output */
+            if (reverseout) {
+                (void)fputs(rtodms(pline, data.uv.v, 'N', 'S'), stdout);
+                putchar('\t');
+                (void)fputs(rtodms(pline, data.uv.u, 'E', 'W'), stdout);
+            } else {
+                (void)fputs(rtodms(pline, data.uv.u, 'E', 'W'), stdout);
+                putchar('\t');
+                (void)fputs(rtodms(pline, data.uv.v, 'N', 'S'), stdout);
+            }
+        } else {    /* x-y or decimal degree ascii output, scale if warranted by output units */
+            if (inverse) {
+                if (proj_angular_input(Proj, PJ_FWD)) {
+                    data.uv.v *= RAD_TO_DEG;
+                    data.uv.u *= RAD_TO_DEG;
+                }
+            } else {
+                if (proj_angular_output(Proj, PJ_FWD)) {
+                    data.uv.v *= RAD_TO_DEG;
+                    data.uv.u *= RAD_TO_DEG;
+                }
+            }
+
+            if (reverseout) {
+                (void)printf(oform, data.uv.v); putchar('\t');
+                (void)printf(oform, data.uv.u);
+            } else {
+                (void)printf(oform, data.uv.u); putchar('\t');
+                (void)printf(oform, data.uv.v);
+            }
+        }
+
+        /* print scale factor data */
+        if (dofactors) {
+            if (!facs_bad)
+                (void)printf("\t<%g %g %g %g %g %g>",
+                    facs.meridional_scale, facs.parallel_scale, facs.areal_scale,
+                    facs.angular_distortion * RAD_TO_DEG, facs.tissot_semimajor, facs.tissot_semiminor);
+            else
+                (void)fputs("\t<* * * * * *>", stdout);
+        }
+        (void)fputs(bin_in ? "\n" : s, stdout);
+    }
+}
+
+/* file processing function --- verbosely */
+static void vprocess(FILE *fid) {
+    char line[MAX_LINE+3], *s, pline[40];
+    LP dat_ll;
+    projXY dat_xy;
+    int linvers;
+    PJ_COORD coord;
+
+
+    if (!oform)
+        oform = "%.3f";
+
+    if (bin_in || bin_out)
+        emess(1,"binary I/O not available in -V option");
+
+    for (;;) {
+        proj_errno_reset(Proj);
+        ++emess_dat.File_line;
+
+        if (!(s = fgets(line, MAX_LINE, fid)))
+            break;
+
+        if (!strchr(s, '\n')) { /* overlong line */
+            int c;
+            (void)strcat(s, "\n");
+            /* gobble up to newline */
+            while ((c = fgetc(fid)) != EOF && c != '\n') ;
+        }
+
+        if (*s == tag) { /* pass on data */
+            (void)fputs(s, stdout);
+            continue;
+        }
+
+        /* check to override default input mode */
+        if (*s == 'I' || *s == 'i') {
+            linvers = 1;
+            ++s;
+        } else
+            linvers = inverse;
+
+        if (linvers) {
+            if (!PJ_INVERS(Proj)) {
+                emess(-1,"inverse for this projection not avail.\n");
+                continue;
+            }
+            dat_xy.x = strtod(s, &s);
+            dat_xy.y = strtod(s, &s);
+            if (dat_xy.x == HUGE_VAL || dat_xy.y == HUGE_VAL) {
+                emess(-1,"lon-lat input conversion failure\n");
+                continue;
+            }
+            if (prescale) { dat_xy.x *= fscale; dat_xy.y *= fscale; }
+            if (reversein) {
+                projXY temp = dat_xy;
+                dat_xy.x = temp.y;
+                dat_xy.y = temp.x;
+            }
+            dat_ll = pj_inv(dat_xy, Proj);
+        } else {
+            dat_ll.lam = proj_dmstor(s, &s);
+            dat_ll.phi = proj_dmstor(s, &s);
+            if (dat_ll.lam == HUGE_VAL || dat_ll.phi == HUGE_VAL) {
+                emess(-1,"lon-lat input conversion failure\n");
+                continue;
+            }
+            if (reversein) {
+                LP temp = dat_ll;
+                dat_ll.lam = temp.phi;
+                dat_ll.phi = temp.lam;
+            }
+            dat_xy = pj_fwd(dat_ll, Proj);
+            if (postscale) { dat_xy.x *= fscale; dat_xy.y *= fscale; }
+        }
+
+        /* For some reason pj_errno does not work as expected in some   */
+        /* versions of Visual Studio, so using pj_get_errno_ref instead */
+        if (*pj_get_errno_ref()) {
+            emess(-1, pj_strerrno(*pj_get_errno_ref()));
+            continue;
+        }
+
+        if (!*s && (s > line)) --s; /* assumed we gobbled \n */
+        coord.lp = dat_ll;
+        facs = proj_factors(Proj, coord);
+        if (proj_errno(Proj)) {
+            emess(-1,"failed to compute factors\n\n");
+            continue;
+        }
+
+        if (*s != '\n')
+            (void)fputs(s, stdout);
+
+        (void)fputs("Longitude: ", stdout);
+        (void)fputs(proj_rtodms(pline, dat_ll.lam, 'E', 'W'), stdout);
+        (void)printf(" [ %.11g ]\n", dat_ll.lam * RAD_TO_DEG);
+        (void)fputs("Latitude:  ", stdout);
+        (void)fputs(proj_rtodms(pline, dat_ll.phi, 'N', 'S'), stdout);
+        (void)printf(" [ %.11g ]\n", dat_ll.phi * RAD_TO_DEG);
+        (void)fputs("Easting (x):   ", stdout);
+        (void)printf(oform, dat_xy.x); putchar('\n');
+        (void)fputs("Northing (y):  ", stdout);
+        (void)printf(oform, dat_xy.y); putchar('\n');
+        (void)printf("Meridian scale (h) : %.8f  ( %.4g %% error )\n", facs.meridional_scale, (facs.meridional_scale-1.)*100.);
+        (void)printf("Parallel scale (k) : %.8f  ( %.4g %% error )\n", facs.parallel_scale, (facs.parallel_scale-1.)*100.);
+        (void)printf("Areal scale (s):     %.8f  ( %.4g %% error )\n", facs.areal_scale, (facs.areal_scale-1.)*100.);
+        (void)printf("Angular distortion (w): %.3f\n", facs.angular_distortion * RAD_TO_DEG);
+        (void)printf("Meridian/Parallel angle: %.5f\n", facs.meridian_parallel_angle * RAD_TO_DEG);
+        (void)printf("Convergence : ");
+        (void)fputs(proj_rtodms(pline, facs.meridian_convergence, 0, 0), stdout);
+        (void)printf(" [ %.8f ]\n", facs.meridian_convergence * RAD_TO_DEG);
+        (void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n", facs.tissot_semimajor, facs.tissot_semiminor);
+    }
+}
+
+int main(int argc, char **argv) {
+    char *arg, *pargv[MAX_PARGS], **iargv = argv;
+    char **eargv = argv;
+    FILE *fid;
+    int pargc = 0, iargc = argc, eargc = 0, mon = 0;
+
+    if ( (emess_dat.Prog_name = strrchr(*argv,DIR_CHAR)) != nullptr)
+        ++emess_dat.Prog_name;
+    else emess_dat.Prog_name = *argv;
+    inverse = ! strncmp(emess_dat.Prog_name, "inv", 3);
+    if (argc <= 1 ) {
+        (void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name);
+        exit (0);
+    }
+
+    /* process run line arguments */
+    while (--argc > 0) { /* collect run line arguments */
+        if(**++argv == '-') for(arg = *argv;;) {
+            switch(*++arg) {
+              case '\0': /* position of "stdin" */
+                if (arg[-1] == '-') eargv[eargc++] = const_cast<char*>("-");
+                break;
+              case 'b': /* binary I/O */
+                bin_in = bin_out = 1;
+                continue;
+              case 'v': /* monitor dump of initialization */
+                mon = 1;
+                continue;
+              case 'i': /* input binary */
+                bin_in = 1;
+                continue;
+              case 'o': /* output binary */
+                bin_out = 1;
+                continue;
+              case 'I': /* alt. method to spec inverse */
+                inverse = 1;
+                continue;
+              case 'E': /* echo ascii input to ascii output */
+                echoin = 1;
+                continue;
+              case 'V': /* very verbose processing mode */
+                very_verby = 1;
+                mon = 1;
+                continue;
+              case 'S': /* compute scale factors */
+                dofactors = 1;
+                continue;
+              case 't': /* set col. one char */
+                if (arg[1]) tag = *++arg;
+                else emess(1,"missing -t col. 1 tag");
+                continue;
+              case 'l': /* list projections, ellipses or units */
+                if (!arg[1] || arg[1] == 'p' || arg[1] == 'P') {
+                    /* list projections */
+                    const struct PJ_LIST *lp;
+                    int do_long = arg[1] == 'P', c;
+                    const char *str;
+
+                    for (lp = proj_list_operations() ; lp->id ; ++lp) {
+                        if( strcmp(lp->id,"latlong") == 0
+                            || strcmp(lp->id,"longlat") == 0
+                            || strcmp(lp->id,"geocent") == 0 )
+                            continue;
+
+                        (void)printf("%s : ", lp->id);
+                        if (do_long)  /* possibly multiline description */
+                            (void)puts(*lp->descr);
+                        else { /* first line, only */
+                            str = *lp->descr;
+                            while ((c = *str++) && c != '\n')
+                                putchar(c);
+                            putchar('\n');
+                        }
+                    }
+                } else if (arg[1] == '=') { /* list projection 'descr' */
+                    const struct PJ_LIST *lp;
+
+                    arg += 2;
+                    for (lp = proj_list_operations(); lp->id ; ++lp)
+                        if (!strcmp(lp->id, arg)) {
+                            (void)printf("%9s : %s\n", lp->id, *lp->descr);
+                            break;
+                        }
+                } else if (arg[1] == 'e') { /* list ellipses */
+                    const struct PJ_ELLPS *le;
+
+                    for (le = proj_list_ellps(); le->id ; ++le)
+                        (void)printf("%9s %-16s %-16s %s\n",
+                                     le->id, le->major, le->ell, le->name);
+                } else if (arg[1] == 'u') { /* list units */
+                    const struct PJ_UNITS *lu;
+
+                    for (lu = proj_list_units(); lu->id ; ++lu)
+                        (void)printf("%12s %-20s %s\n",
+                                     lu->id, lu->to_meter, lu->name);
+                } else if (arg[1] == 'd') { /* list datums */
+                    const struct PJ_DATUMS *ld;
+
+                    printf("__datum_id__ __ellipse___ __definition/comments______________________________\n" );
+                    for (ld = pj_get_datums_ref(); ld->id ; ++ld)
+                    {
+                        printf("%12s %-12s %-30s\n",
+                               ld->id, ld->ellipse_id, ld->defn);
+                        if( ld->comments != nullptr && strlen(ld->comments) > 0 )
+                            printf( "%25s %s\n", " ", ld->comments );
+                    }
+                } else
+                    emess(1,"invalid list option: l%c",arg[1]);
+                exit(0);
+                /* cppcheck-suppress duplicateBreak */
+                continue; /* artificial */
+              case 'e': /* error line alternative */
+                if (--argc <= 0)
+                    noargument:
+                emess(1,"missing argument for -%c",*arg);
+                oterr = *++argv;
+                continue;
+              case 'T': /* generate Chebyshev coefficients */
+                if (--argc <= 0) goto noargument;
+                cheby_str = *++argv;
+                continue;
+              case 'm': /* cartesian multiplier */
+                if (--argc <= 0) goto noargument;
+                postscale = 1;
+                if (!strncmp("1/",*++argv,2) ||
+                    !strncmp("1:",*argv,2)) {
+                    if((fscale = atof((*argv)+2)) == 0.)
+                        goto badscale;
+                    fscale = 1. / fscale;
+                } else
+                    if ((fscale = atof(*argv)) == 0.) {
+                      badscale:
+                        emess(1,"invalid scale argument");
+                    }
+                continue;
+              case 'W': /* specify seconds precision */
+              case 'w': /* -W for constant field width */
+              {
+                int c = arg[1];
+                if (c != 0 && isdigit(c)) {
+                    set_rtodms(c - '0', *arg == 'W');
+                    ++arg;
+                } else
+                    emess(1,"-W argument missing or non-digit");
+                continue;
+              }
+              case 'f': /* alternate output format degrees or xy */
+                if (--argc <= 0) goto noargument;
+                oform = *++argv;
+                continue;
+              case 'd':
+                if (--argc <= 0) goto noargument;
+                sprintf(oform_buffer, "%%.%df", atoi(*++argv));
+                oform = oform_buffer;
+                break;
+              case 'r': /* reverse input */
+                reversein = 1;
+                continue;
+              case 's': /* reverse output */
+                reverseout = 1;
+                continue;
+              default:
+                emess(1, "invalid option: -%c",*arg);
+                break;
+            }
+            break;
+        } else if (**argv == '+') { /* + argument */
+            if (pargc < MAX_PARGS)
+                pargv[pargc++] = *argv + 1;
+            else
+                emess(1,"overflowed + argument table");
+        } else /* assumed to be input file name(s) */
+            eargv[eargc++] = *argv;
+    }
+    if (eargc == 0 && !cheby_str) /* if no specific files force sysin */
+        eargv[eargc++] = const_cast<char*>("-");
+    else if (eargc > 0 && cheby_str) /* warning */
+        emess(4, "data files when generating Chebychev prohibited");
+    /* done with parameter and control input */
+    if (inverse && postscale) {
+        prescale = 1;
+        postscale = 0;
+        fscale = 1./fscale;
+    }
+    if (!(Proj = pj_init(pargc, pargv)))
+        emess(3,"projection initialization failure\ncause: %s",
+              pj_strerrno(pj_errno));
+
+    if (!proj_angular_input(Proj, PJ_FWD)) {
+        emess(3, "can't initialize operations that take non-angular input coordinates");
+        exit(0);
+    }
+
+    if (proj_angular_output(Proj, PJ_FWD)) {
+        emess(3, "can't initialize operations that produce angular output coordinates");
+        exit(0);
+    }
+
+    if (inverse) {
+        if (!Proj->inv)
+            emess(3,"inverse projection not available");
+        proj.inv = pj_inv;
+    } else
+        proj.fwd = pj_fwd;
+    if (cheby_str) {
+        gen_cheb(inverse, int_proj, cheby_str, Proj, iargc, iargv);
+        exit(0);
+    }
+
+    /* set input formatting control */
+    if (mon) {
+        pj_pr_list(Proj);
+        if (very_verby) {
+            (void)printf("#Final Earth figure: ");
+            if (Proj->es != 0.0) {
+                (void)printf("ellipsoid\n#  Major axis (a): ");
+                (void)printf(oform ? oform : "%.3f", Proj->a);
+                (void)printf("\n#  1/flattening: %.6f\n",
+                             1./(1. - sqrt(1. - Proj->es)));
+                (void)printf("#  squared eccentricity: %.12f\n", Proj->es);
+            } else {
+                (void)printf("sphere\n#  Radius: ");
+                (void)printf(oform ? oform : "%.3f", Proj->a);
+                (void)putchar('\n');
+            }
+        }
+    }
+
+    if (inverse)
+        informat = strtod;
+    else {
+        informat = proj_dmstor;
+        if (!oform)
+            oform = "%.2f";
+    }
+
+    if (bin_out) {
+        SET_BINARY_MODE(stdout);
+    }
+
+    /* process input file list */
+    for ( ; eargc-- ; ++eargv) {
+        if (**eargv == '-') {
+            fid = stdin;
+            emess_dat.File_name = const_cast<char*>("<stdin>");
+
+            if (bin_in)
+            {
+                SET_BINARY_MODE(stdin);
+            }
+
+        } else {
+            if ((fid = fopen(*eargv, "rb")) == nullptr) {
+                emess(-2, *eargv, "input file");
+                continue;
+            }
+            emess_dat.File_name = *eargv;
+        }
+        emess_dat.File_line = 0;
+        if (very_verby)
+            vprocess(fid);
+        else
+            process(fid);
+        (void)fclose(fid);
+        emess_dat.File_name = nullptr;
+    }
+
+    if( Proj )
+        pj_free(Proj);
+
+    exit(0); /* normal completion */
+}
diff --git a/src/apps/proj_strtod.cpp b/src/apps/proj_strtod.cpp
new file mode 100644
index 00000000..b8edc6a3
--- /dev/null
+++ b/src/apps/proj_strtod.cpp
@@ -0,0 +1,440 @@
+/***********************************************************************
+
+ proj_strtod: Convert string to double, accepting underscore separators
+
+                    Thomas Knudsen, 2017-01-17/09-19
+
+************************************************************************
+
+Conventionally, PROJ.4 does not honor locale settings, consistently
+behaving as if LC_ALL=C.
+
+For this to work, we have, for many years, been using other solutions
+than the C standard library strtod/atof functions for converting strings
+to doubles.
+
+In the early versions of proj, iirc, a gnu version of strtod was used,
+mostly to work around cases where the same system library was used for
+C and Fortran linking, hence making strtod accept "D" and "d" as
+exponentiation indicators, following Fortran Double Precision constant
+syntax. This broke the proj angular syntax, accepting a "d" to mean
+"degree": 12d34'56", meaning 12 degrees 34 minutes and 56 seconds.
+
+With an explicit MIT licence, PROJ.4 could not include GPL code any
+longer, and apparently at some time, the GPL code was replaced by the
+current C port of a GDAL function (in pj_strtod.c), which reads the
+LC_NUMERIC setting and, behind the back of the user, momentarily changes
+the conventional '.' delimiter to whatever the locale requires, then
+calls the system supplied strtod.
+
+While this requires a minimum amount of coding, it only solves one
+problem, and not in a very generic way.
+
+Another problem, I would like to see solved, is the handling of underscores
+as generic delimiters. This is getting popular in a number of programming
+languages (Ada, C++, C#, D, Java, Julia, Perl 5, Python, Rust, etc.
+cf. e.g. https://www.python.org/dev/peps/pep-0515/), and in our case of
+handling numbers being in the order of magnitude of the Earth's dimensions,
+and a resolution of submillimetre, i.e. having 10 or more significant digits,
+splitting the "wall of digits" into smaller chunks is of immense value.
+
+Hence this reimplementation of strtod, which hardcodes '.' as indicator of
+numeric fractions, and accepts '_' anywhere in a numerical string sequence:
+So a typical northing value can be written
+
+                            6_098_907.8250 m
+rather than
+                             6098907.8250 m
+
+which, in my humble opinion, is well worth the effort.
+
+While writing this code, I took ample inspiration from Michael Ringgaard's
+strtod version over at http://www.jbox.dk/sanos/source/lib/strtod.c.html,
+and Yasuhiro Matsumoto's public domain version over at
+https://gist.github.com/mattn/1890186. The code below is, however, not
+copied from any of the two mentioned - it is a reimplementation, and
+probably suffers from its own set of bugs. So for now, it is intended
+not as a replacement of pj_strtod, but only as an experimental piece of
+code for use in an experimental new transformation program, cct.
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18
+
+************************************************************************
+
+* Copyright (c) 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.
+
+***********************************************************************/
+
+#include "proj_strtod.h"
+
+#include <stdlib.h> /* for abs */
+#include <string.h> /* for strchr */
+#include <errno.h>
+#include <ctype.h>
+#include <float.h>  /* for HUGE_VAL */
+#include <math.h>   /* for pow() */
+
+
+double proj_strtod(const char *str, char **endptr) {
+    double number = 0, integral_part = 0;
+    int exponent = 0;
+    int fraction_is_nonzero = 0;
+    int sign = 0;
+    char *p = (char *) str;
+    int n = 0;
+    int num_digits_total        = 0;
+    int num_digits_after_comma  = 0;
+    int num_prefixed_zeros      = 0;
+
+    if (nullptr==str) {
+        errno = EFAULT;
+        if (endptr)
+            *endptr = p;
+        return HUGE_VAL;
+    }
+
+    /* First skip leading whitespace */
+    while (isspace(*p))
+        p++;
+
+    /* Empty string? */
+    if (0==*p) {
+        if (endptr)
+            *endptr = (char *) str;
+        return 0;
+    }
+
+    /* non-numeric? */
+    if (nullptr==strchr("0123456789+-._", *p)) {
+        if (endptr)
+            *endptr = (char *) str;
+        return 0;
+    }
+
+    /* Then handle optional prefixed sign and skip prefix zeros */
+    switch (*p) {
+        case '-':
+            sign = -1;
+            p++;
+            break;
+        case '+':
+            sign =  1;
+            p++;
+            break;
+        default:
+            if (isdigit(*p) || '_'==*p || '.'==*p)
+                break;
+            if (endptr)
+                *endptr = (char *) str;
+            return 0;
+    }
+
+    /* stray sign, as in "+/-"? */
+    if (0!=sign && (nullptr==strchr ("0123456789._", *p) || 0==*p)) {
+        if (endptr)
+            *endptr = (char *) str;
+        return 0;
+    }
+
+    /* skip prefixed zeros before '.' */
+    while ('0'==*p || '_'==*p)
+        p++;
+
+    /* zero? */
+    if ((0==*p) || nullptr==strchr ("0123456789eE.", *p) || isspace(*p)) {
+        if (endptr)
+            *endptr = p;
+        return sign==-1? -0: 0;
+    }
+
+    /* Now expect a (potentially zero-length) string of digits */
+    while (isdigit(*p) || ('_'==*p)) {
+        if ('_'==*p) {
+            p++;
+            continue;
+        }
+        number = number * 10. + (*p - '0');
+        p++;
+        num_digits_total++;
+    }
+    integral_part = number;
+
+    /* Done? */
+    if (0==*p) {
+        if (endptr)
+            *endptr = p;
+        if (sign==-1)
+            return -number;
+        return number;
+    }
+
+    /* Do we have a fractional part? */
+    if ('.'==*p) {
+        p++;
+
+        /* keep on skipping prefixed zeros (i.e. allow writing 1e-20 */
+        /* as 0.00000000000000000001 without losing precision) */
+        if (0==integral_part)
+            while ('0'==*p || '_'==*p) {
+                if ('0'==*p)
+                    num_prefixed_zeros++;
+                p++;
+            }
+
+        /* if the next character is nonnumeric, we have reached the end */
+        if (0==*p || nullptr==strchr ("_0123456789eE+-", *p)) {
+            if (endptr)
+                *endptr = p;
+            if (sign==-1)
+                return -number;
+            return number;
+        }
+
+        while (isdigit(*p) || '_'==*p) {
+            /* Don't let pathologically long fractions destroy precision */
+            if ('_'==*p || num_digits_total > 17) {
+                p++;
+                continue;
+            }
+
+            number = number * 10. + (*p - '0');
+            if (*p!='0')
+                fraction_is_nonzero = 1;
+            p++;
+            num_digits_total++;
+            num_digits_after_comma++;
+        }
+
+        /* Avoid having long zero-tails (4321.000...000) destroy precision */
+        if (fraction_is_nonzero)
+            exponent = -(num_digits_after_comma + num_prefixed_zeros);
+        else
+            number = integral_part;
+    }  /* end of fractional part */
+
+
+    /* non-digit */
+    if (0==num_digits_total) {
+        errno = EINVAL;
+        if (endptr)
+            *endptr = p;
+        return HUGE_VAL;
+    }
+
+    if (sign==-1)
+        number = -number;
+
+    /* Do we have an exponent part? */
+    while (*p == 'e' || *p == 'E') {
+        p++;
+
+        /* Just a stray "e", as in 100elephants? */
+        if (0==*p || nullptr==strchr ("0123456789+-_", *p)) {
+            p--;
+            break;
+        }
+
+        while ('_'==*p)
+            p++;
+        /* Does it have a sign? */
+        sign = 0;
+        if ('-'==*p)
+            sign = -1;
+        if ('+'==*p)
+            sign = +1;
+        if (0==sign) {
+            if (!isdigit(*p) && *p!='_') {
+                if (endptr)
+                    *endptr = p;
+                return HUGE_VAL;
+            }
+        }
+        else
+            p++;
+
+
+        /* Go on and read the exponent */
+        n = 0;
+        while (isdigit(*p) || '_'==*p) {
+            if ('_'==*p) {
+                p++;
+                continue;
+            }
+            n = n * 10 + (*p - '0');
+            p++;
+        }
+
+        if (-1==sign)
+            n = -n;
+        exponent += n;
+        break;
+    }
+
+    if (endptr)
+        *endptr = p;
+
+    if ((exponent < DBL_MIN_EXP) || (exponent > DBL_MAX_EXP)) {
+        errno = ERANGE;
+        return HUGE_VAL;
+    }
+
+    /* on some platforms pow() is very slow - so don't call it if exponent is close to 0 */
+    if (0==exponent)
+        return number;
+    if (abs (exponent) < 20) {
+        double ex = 1;
+        int absexp = exponent < 0? -exponent: exponent;
+        while (absexp--)
+            ex *= 10;
+        number = exponent < 0? number / ex: number * ex;
+    }
+    else
+        number *= pow (10, exponent);
+
+    return number;
+}
+
+double proj_atof(const char *str) {
+    return proj_strtod(str, nullptr);
+}
+
+#ifdef TEST
+
+/* compile/run: gcc -DTEST -o proj_strtod_test proj_strtod.c  &&  proj_strtod_test */
+
+#include <stdio.h>
+
+char *un_underscore (char *s) {
+    static char u[1024];
+    int i, m,  n;
+    for (i = m = 0, n = strlen (s);  i < n;  i++) {
+        if (s[i]=='_') {
+            m++;
+            continue;
+        }
+        u[i - m] = s[i];
+    }
+    u[n-m] = 0;
+    return u;
+}
+
+int thetest (char *s, int line) {
+    char *endp, *endq, *u;
+    double p, q;
+    int errnop, errnoq, prev_errno;
+
+    prev_errno = errno;
+
+    u = un_underscore (s);
+
+    errno = 0;
+    p = proj_strtod (s, &endp);
+    errnop = errno;
+    errno = 0;
+    q = strtod (u, &endq);
+    errnoq = errno;
+
+    errno = prev_errno;
+
+    if (q==p && 0==strcmp (endp, endq) && errnop==errnoq)
+        return 0;
+
+    errno = line;
+    printf ("Line: %3.3d  -  [%s] [%s]\n", line, s, u);
+    printf ("proj_strtod: %2d %.17g  [%s]\n", errnop, p, endp);
+    printf ("libc_strtod: %2d %.17g  [%s]\n", errnoq, q, endq);
+    return 1;
+}
+
+#define test(s) thetest(s, __LINE__)
+
+int main (int argc, char **argv) {
+    double res;
+    char *endptr;
+
+    errno = 0;
+
+    test ("");
+    test ("     ");
+    test ("     abcde");
+    test ("     edcba");
+    test ("abcde");
+    test ("edcba");
+    test ("+");
+    test ("-");
+    test ("+ ");
+    test ("- ");
+    test (" + ");
+    test (" - ");
+    test ("e 1");
+    test ("e1");
+    test ("0 66");
+    test ("1.");
+    test ("0.");
+    test ("1.0");
+    test ("0.0");
+    test ("1 ");
+    test ("0 ");
+    test ("-0 ");
+    test ("0_ ");
+    test ("0_");
+    test ("1e");
+    test ("_1.0");
+    test ("_0.0");
+    test ("1_.0");
+    test ("0_.0");
+    test ("1__.0");
+    test ("0__.0");
+    test ("1.__0");
+    test ("0.__0");
+    test ("1.0___");
+    test ("0.0___");
+    test ("1e2");
+    test ("__123_456_789_._10_11_12");
+    test ("1______");
+    test ("1___e__2__");
+    test ("-1");
+    test ("-1.0");
+    test ("-0");
+    test ("-1e__-_2__rest");
+    test ("0.00002");
+    test ("0.00001");
+    test ("-0.00002");
+    test ("-0.00001");
+    test ("-0.00001e-2");
+    test ("-0.00001e2");
+    test ("1e9999");
+
+    /* We expect this one to differ */
+    test ("0.000000000000000000000000000000000000000000000000000000000000000000000000002");
+
+    if (errno)
+        printf ("First discrepancy in line %d\n", errno);
+
+    if (argc < 2)
+        return 0;
+    res =  proj_strtod (argv[1], &endptr);
+    printf ("res = %20.15g. Rest = [%s],  errno = %d\n", res, endptr, (int) errno);
+    return 0;
+}
+#endif
diff --git a/src/apps/proj_strtod.h b/src/apps/proj_strtod.h
new file mode 100644
index 00000000..38c2d1f4
--- /dev/null
+++ b/src/apps/proj_strtod.h
@@ -0,0 +1,4 @@
+/* Internal header for proj_strtod.c */
+
+double proj_strtod(const char *str, char **endptr);
+double proj_atof(const char *str);
diff --git a/src/apps/projinfo.cpp b/src/apps/projinfo.cpp
new file mode 100644
index 00000000..d604365a
--- /dev/null
+++ b/src/apps/projinfo.cpp
@@ -0,0 +1,1067 @@
+/******************************************************************************
+ *
+ * Project:  PROJ
+ * Purpose:  projinfo utility
+ * Author:   Even Rouault <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2018, Even Rouault <even dot rouault at spatialys dot com>
+ *
+ * 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.
+ ****************************************************************************/
+
+//! @cond Doxygen_Suppress
+
+#define FROM_PROJ_CPP
+
+#include <cstdlib>
+#include <fstream> // std::ifstream
+#include <iostream>
+#include <utility>
+
+#include "projects.h"
+
+#include <proj/common.hpp>
+#include <proj/coordinateoperation.hpp>
+#include <proj/crs.hpp>
+#include <proj/io.hpp>
+#include <proj/metadata.hpp>
+#include <proj/util.hpp>
+
+#include "proj/internal/internal.hpp" // for split
+
+using namespace NS_PROJ::common;
+using namespace NS_PROJ::crs;
+using namespace NS_PROJ::io;
+using namespace NS_PROJ::metadata;
+using namespace NS_PROJ::operation;
+using namespace NS_PROJ::util;
+using namespace NS_PROJ::internal;
+
+// ---------------------------------------------------------------------------
+
+namespace { // anonymous namespace
+struct OutputOptions {
+    bool quiet = false;
+    bool PROJ5 = false;
+    bool PROJ4 = false;
+    bool WKT2_2018 = false;
+    bool WKT2_2018_SIMPLIFIED = false;
+    bool WKT2_2015 = false;
+    bool WKT2_2015_SIMPLIFIED = false;
+    bool WKT1_GDAL = false;
+    bool WKT1_ESRI = false;
+    bool c_ify = false;
+    bool singleLine = false;
+};
+} // anonymous namespace
+
+// ---------------------------------------------------------------------------
+
+static void usage() {
+    std::cerr
+        << "usage: projinfo [-o formats] [-k crs|operation] [--summary] [-q]"
+        << std::endl
+        << "                ([--area name_or_code] | "
+           "[--bbox min_long,min_lat,max_long,max_lat]) "
+        << std::endl
+        << "                [--spatial-test contains|intersects]" << std::endl
+        << "                [--crs-extent-use none|both|intersection|smallest]"
+        << std::endl
+        << "                [--grid-check none|discard_missing|sort] "
+           "[--show-superseded]"
+        << std::endl
+        << "                [--pivot-crs none|{auth:code[,auth:code]*}]"
+        << std::endl
+        << "                [--boundcrs-to-wgs84]" << std::endl
+        << "                [--main-db-path path] [--aux-db-path path]*"
+        << std::endl
+        << "                [--identify]" << std::endl
+        << "                [--c-ify] [--single-line]" << std::endl
+        << "                {object_definition} | (-s {srs_def} -t {srs_def})"
+        << std::endl;
+    std::cerr << std::endl;
+    std::cerr << "-o: formats is a comma separated combination of: "
+                 "all,default,PROJ4,PROJ,WKT_ALL,WKT2_2015,WKT2_2018,WKT1_GDAL,"
+                 "WKT1_ESRI"
+              << std::endl;
+    std::cerr << "    Except 'all' and 'default', other format can be preceded "
+                 "by '-' to disable them"
+              << std::endl;
+    std::cerr << std::endl;
+    std::cerr << "{object_definition} might be a PROJ string, a WKT string, "
+                 " a AUTHORITY:CODE, or urn:ogc:def:OBJECT_TYPE:AUTHORITY::CODE"
+              << std::endl;
+    std::exit(1);
+}
+
+// ---------------------------------------------------------------------------
+
+static std::string un_c_ify_string(const std::string &str) {
+    std::string out(str);
+    out = out.substr(1, out.size() - 2);
+    out = replaceAll(out, "\\\"", "{ESCAPED_DOUBLE_QUOTE}");
+    out = replaceAll(out, "\\n\"", "");
+    out = replaceAll(out, "\"", "");
+    out = replaceAll(out, "{ESCAPED_DOUBLE_QUOTE}", "\"");
+    return out;
+}
+
+// ---------------------------------------------------------------------------
+
+static std::string c_ify_string(const std::string &str) {
+    std::string out(str);
+    out = replaceAll(out, "\"", "{DOUBLE_QUOTE}");
+    out = replaceAll(out, "\n", "\\n\"\n\"");
+    out = replaceAll(out, "{DOUBLE_QUOTE}", "\\\"");
+    return "\"" + out + "\"";
+}
+
+// ---------------------------------------------------------------------------
+
+static BaseObjectNNPtr buildObject(DatabaseContextPtr dbContext,
+                                   const std::string &user_string,
+                                   bool kindIsCRS, const std::string &context,
+                                   bool buildBoundCRSToWGS84, bool allowPivots,
+                                   bool quiet) {
+    BaseObjectPtr obj;
+
+    std::string l_user_string(user_string);
+    if (!user_string.empty() && user_string[0] == '@') {
+        std::ifstream fs;
+        auto filename = user_string.substr(1);
+        fs.open(filename, std::fstream::in | std::fstream::binary);
+        if (!fs.is_open()) {
+            std::cerr << context << ": cannot open " << filename << std::endl;
+            std::exit(1);
+        }
+        l_user_string.clear();
+        while (!fs.eof()) {
+            char buffer[256];
+            fs.read(buffer, sizeof(buffer));
+            l_user_string.append(buffer, static_cast<size_t>(fs.gcount()));
+            if (l_user_string.size() > 100 * 1000) {
+                fs.close();
+                std::cerr << context << ": too big file" << std::endl;
+                std::exit(1);
+            }
+        }
+        fs.close();
+    }
+    if (!l_user_string.empty() && l_user_string.back() == '\n') {
+        l_user_string.resize(l_user_string.size() - 1);
+    }
+    if (!l_user_string.empty() && l_user_string.back() == '\r') {
+        l_user_string.resize(l_user_string.size() - 1);
+    }
+
+    try {
+        auto tokens = split(l_user_string, ':');
+        if (!kindIsCRS && tokens.size() == 2) {
+            auto urn = "urn:ogc:def:coordinateOperation:" + tokens[0] + "::" +
+                       tokens[1];
+            obj = createFromUserInput(urn, dbContext).as_nullable();
+        } else {
+            // Convenience to be able to use C escaped strings...
+            if (l_user_string.size() > 2 && l_user_string[0] == '"' &&
+                l_user_string.back() == '"' &&
+                l_user_string.find("\\\"") != std::string::npos) {
+                l_user_string = un_c_ify_string(l_user_string);
+            }
+            WKTParser wktParser;
+            if (wktParser.guessDialect(l_user_string) !=
+                WKTParser::WKTGuessedDialect::NOT_WKT) {
+                wktParser.setStrict(false);
+                wktParser.attachDatabaseContext(dbContext);
+                obj = wktParser.createFromWKT(l_user_string).as_nullable();
+                if (!quiet) {
+                    auto warnings = wktParser.warningList();
+                    if (!warnings.empty()) {
+                        for (const auto &str : warnings) {
+                            std::cerr << "Warning: " << str << std::endl;
+                        }
+                    }
+                }
+            } else {
+                obj =
+                    createFromUserInput(l_user_string, dbContext).as_nullable();
+            }
+        }
+    } catch (const std::exception &e) {
+        std::cerr << context << ": parsing of user string failed: " << e.what()
+                  << std::endl;
+        std::exit(1);
+    }
+
+    if (buildBoundCRSToWGS84) {
+        auto crs = std::dynamic_pointer_cast<CRS>(obj);
+        if (crs) {
+            obj = crs->createBoundCRSToWGS84IfPossible(dbContext, allowPivots)
+                      .as_nullable();
+        }
+    }
+
+    return NN_NO_CHECK(obj);
+}
+
+// ---------------------------------------------------------------------------
+
+static void outputObject(DatabaseContextPtr dbContext, BaseObjectNNPtr obj,
+                         bool allowPivots, const OutputOptions &outputOpt) {
+
+    auto identified = dynamic_cast<const IdentifiedObject *>(obj.get());
+    if (!outputOpt.quiet && identified && identified->isDeprecated()) {
+        std::cout << "Warning: object is deprecated" << std::endl;
+        auto crs = dynamic_cast<const CRS *>(obj.get());
+        if (crs && dbContext) {
+            try {
+                auto list = crs->getNonDeprecated(NN_NO_CHECK(dbContext));
+                if (!list.empty()) {
+                    std::cout << "Alternative non-deprecated CRS:" << std::endl;
+                }
+                for (const auto &altCRS : list) {
+                    const auto &ids = altCRS->identifiers();
+                    if (!ids.empty()) {
+                        std::cout << "  " << *(ids[0]->codeSpace()) << ":"
+                                  << ids[0]->code() << std::endl;
+                    }
+                }
+            } catch (const std::exception &) {
+            }
+        }
+        std::cout << std::endl;
+    }
+
+    auto projStringExportable =
+        nn_dynamic_pointer_cast<IPROJStringExportable>(obj);
+    bool alreadyOutputed = false;
+    if (projStringExportable) {
+        if (outputOpt.PROJ5) {
+            try {
+                if (!outputOpt.quiet) {
+                    std::cout << "PROJ string:" << std::endl;
+                }
+                std::cout << projStringExportable->exportToPROJString(
+                                 PROJStringFormatter::create(
+                                     PROJStringFormatter::Convention::PROJ_5,
+                                     dbContext)
+                                     .get())
+                          << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to PROJ string: " << e.what()
+                          << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.PROJ4) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "PROJ.4 string:" << std::endl;
+                }
+
+                auto crs = nn_dynamic_pointer_cast<CRS>(obj);
+                std::shared_ptr<IPROJStringExportable> objToExport;
+                if (crs) {
+                    objToExport =
+                        nn_dynamic_pointer_cast<IPROJStringExportable>(
+                            crs->createBoundCRSToWGS84IfPossible(dbContext,
+                                                                 allowPivots));
+                }
+                if (!objToExport) {
+                    objToExport = projStringExportable;
+                }
+
+                std::cout << objToExport->exportToPROJString(
+                                 PROJStringFormatter::create(
+                                     PROJStringFormatter::Convention::PROJ_4,
+                                     dbContext)
+                                     .get())
+                          << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to PROJ string: " << e.what()
+                          << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+    }
+
+    auto wktExportable = nn_dynamic_pointer_cast<IWKTExportable>(obj);
+    if (wktExportable) {
+        if (outputOpt.WKT2_2015) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT2_2015 string:" << std::endl;
+                }
+                auto formatter =
+                    WKTFormatter::create(WKTFormatter::Convention::WKT2_2015);
+                if (outputOpt.singleLine) {
+                    formatter->setMultiLine(false);
+                }
+                auto wkt = wktExportable->exportToWKT(formatter.get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT2_2015: " << e.what()
+                          << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.WKT2_2015_SIMPLIFIED) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT2_2015_SIMPLIFIED string:" << std::endl;
+                }
+                auto formatter = WKTFormatter::create(
+                    WKTFormatter::Convention::WKT2_2015_SIMPLIFIED);
+                if (outputOpt.singleLine) {
+                    formatter->setMultiLine(false);
+                }
+                auto wkt = wktExportable->exportToWKT(formatter.get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT2_2015_SIMPLIFIED: "
+                          << e.what() << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.WKT2_2018) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT2_2018 string:" << std::endl;
+                }
+                auto formatter =
+                    WKTFormatter::create(WKTFormatter::Convention::WKT2_2018);
+                if (outputOpt.singleLine) {
+                    formatter->setMultiLine(false);
+                }
+                auto wkt = wktExportable->exportToWKT(formatter.get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT2_2018: " << e.what()
+                          << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.WKT2_2018_SIMPLIFIED) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT2_2018_SIMPLIFIED string:" << std::endl;
+                }
+                auto formatter = WKTFormatter::create(
+                    WKTFormatter::Convention::WKT2_2018_SIMPLIFIED);
+                if (outputOpt.singleLine) {
+                    formatter->setMultiLine(false);
+                }
+                auto wkt = wktExportable->exportToWKT(formatter.get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT2_2018_SIMPLIFIED: "
+                          << e.what() << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.WKT1_GDAL && !nn_dynamic_pointer_cast<Conversion>(obj)) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT1_GDAL:" << std::endl;
+                }
+
+                auto crs = nn_dynamic_pointer_cast<CRS>(obj);
+                std::shared_ptr<IWKTExportable> objToExport;
+                if (crs) {
+                    objToExport = nn_dynamic_pointer_cast<IWKTExportable>(
+                        crs->createBoundCRSToWGS84IfPossible(dbContext,
+                                                             allowPivots));
+                }
+                if (!objToExport) {
+                    objToExport = wktExportable;
+                }
+
+                auto formatter =
+                    WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL);
+                if (outputOpt.singleLine) {
+                    formatter->setMultiLine(false);
+                }
+                auto wkt = objToExport->exportToWKT(formatter.get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+                std::cout << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT1_GDAL: " << e.what()
+                          << std::endl;
+            }
+            alreadyOutputed = true;
+        }
+
+        if (outputOpt.WKT1_ESRI && !nn_dynamic_pointer_cast<Conversion>(obj)) {
+            try {
+                if (alreadyOutputed) {
+                    std::cout << std::endl;
+                }
+                if (!outputOpt.quiet) {
+                    std::cout << "WKT1_ESRI:" << std::endl;
+                }
+
+                auto wkt = wktExportable->exportToWKT(
+                    WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI,
+                                         dbContext)
+                        .get());
+                if (outputOpt.c_ify) {
+                    wkt = c_ify_string(wkt);
+                }
+                std::cout << wkt << std::endl;
+                std::cout << std::endl;
+            } catch (const std::exception &e) {
+                std::cerr << "Error when exporting to WKT1_ESRI: " << e.what()
+                          << std::endl;
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+static void outputOperationSummary(const CoordinateOperationNNPtr &op) {
+    auto ids = op->identifiers();
+    if (!ids.empty()) {
+        std::cout << *(ids[0]->codeSpace()) << ":" << ids[0]->code();
+    } else {
+        std::cout << "unknown id";
+    }
+
+    std::cout << ", ";
+
+    auto name = op->nameStr();
+    if (!name.empty()) {
+        std::cout << name;
+    } else {
+        std::cout << "unknown name";
+    }
+
+    std::cout << ", ";
+
+    auto accuracies = op->coordinateOperationAccuracies();
+    if (!accuracies.empty()) {
+        std::cout << accuracies[0]->value() << " m";
+    } else {
+        if (std::dynamic_pointer_cast<Conversion>(op.as_nullable())) {
+            std::cout << "0 m";
+        } else {
+            std::cout << "unknown accuracy";
+        }
+    }
+
+    std::cout << ", ";
+
+    auto domains = op->domains();
+    if (!domains.empty() && domains[0]->domainOfValidity() &&
+        domains[0]->domainOfValidity()->description().has_value()) {
+        std::cout << *(domains[0]->domainOfValidity()->description());
+    } else {
+        std::cout << "unknown domain of validity";
+    }
+
+    std::cout << std::endl;
+}
+
+// ---------------------------------------------------------------------------
+
+static void outputOperations(
+    DatabaseContextPtr dbContext, const std::string &sourceCRSStr,
+    const std::string &targetCRSStr, const ExtentPtr &bboxFilter,
+    CoordinateOperationContext::SpatialCriterion spatialCriterion,
+    CoordinateOperationContext::SourceTargetCRSExtentUse crsExtentUse,
+    CoordinateOperationContext::GridAvailabilityUse gridAvailabilityUse,
+    bool allowPivots,
+    const std::vector<std::pair<std::string, std::string>> &pivots,
+    const std::string &authority, bool usePROJGridAlternatives,
+    bool showSuperseded, const OutputOptions &outputOpt, bool summary) {
+    auto sourceObj = buildObject(dbContext, sourceCRSStr, true, "source CRS",
+                                 false, false, outputOpt.quiet);
+    auto sourceCRS = nn_dynamic_pointer_cast<CRS>(sourceObj);
+    if (!sourceCRS) {
+        std::cerr << "source CRS string is not a CRS" << std::endl;
+        std::exit(1);
+    }
+
+    auto targetObj = buildObject(dbContext, targetCRSStr, true, "target CRS",
+                                 false, false, outputOpt.quiet);
+    auto targetCRS = nn_dynamic_pointer_cast<CRS>(targetObj);
+    if (!targetCRS) {
+        std::cerr << "target CRS string is not a CRS" << std::endl;
+        std::exit(1);
+    }
+
+    std::vector<CoordinateOperationNNPtr> list;
+    try {
+        auto authFactory =
+            dbContext
+                ? AuthorityFactory::create(NN_NO_CHECK(dbContext), authority)
+                      .as_nullable()
+                : nullptr;
+        auto ctxt =
+            CoordinateOperationContext::create(authFactory, bboxFilter, 0);
+        ctxt->setSpatialCriterion(spatialCriterion);
+        ctxt->setSourceAndTargetCRSExtentUse(crsExtentUse);
+        ctxt->setGridAvailabilityUse(gridAvailabilityUse);
+        ctxt->setAllowUseIntermediateCRS(allowPivots);
+        ctxt->setIntermediateCRS(pivots);
+        ctxt->setUsePROJAlternativeGridNames(usePROJGridAlternatives);
+        ctxt->setDiscardSuperseded(!showSuperseded);
+        list = CoordinateOperationFactory::create()->createOperations(
+            NN_NO_CHECK(sourceCRS), NN_NO_CHECK(targetCRS), ctxt);
+    } catch (const std::exception &e) {
+        std::cerr << "createOperations() failed with: " << e.what()
+                  << std::endl;
+        std::exit(1);
+    }
+    if (outputOpt.quiet && !list.empty()) {
+        outputObject(dbContext, list[0], allowPivots, outputOpt);
+        return;
+    }
+    if (summary) {
+        std::cout << "Candidate operations found: " << list.size() << std::endl;
+        for (const auto &op : list) {
+            outputOperationSummary(op);
+        }
+    } else {
+        bool first = true;
+        for (size_t i = 0; i < list.size(); ++i) {
+            const auto &op = list[i];
+            if (list.size() > 1) {
+                if (!first) {
+                    std::cout << std::endl;
+                }
+                first = false;
+                std::cout << "-------------------------------------"
+                          << std::endl;
+                std::cout << "Operation n"
+                             "\xC2\xB0"
+                          << (i + 1) << ":" << std::endl
+                          << std::endl;
+            }
+            outputOperationSummary(op);
+            std::cout << std::endl;
+            outputObject(dbContext, op, allowPivots, outputOpt);
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+int main(int argc, char **argv) {
+
+    if (argc == 1) {
+        std::cerr << pj_get_release() << std::endl;
+        usage();
+    }
+
+    std::string user_string;
+    bool user_string_specified = false;
+    std::string sourceCRSStr;
+    std::string targetCRSStr;
+    bool outputSwithSpecified = false;
+    OutputOptions outputOpt;
+    bool kindIsCRS = true;
+    bool summary = false;
+    ExtentPtr bboxFilter = nullptr;
+    std::string area;
+    CoordinateOperationContext::SpatialCriterion spatialCriterion =
+        CoordinateOperationContext::SpatialCriterion::STRICT_CONTAINMENT;
+    CoordinateOperationContext::SourceTargetCRSExtentUse crsExtentUse =
+        CoordinateOperationContext::SourceTargetCRSExtentUse::SMALLEST;
+    bool buildBoundCRSToWGS84 = false;
+    CoordinateOperationContext::GridAvailabilityUse gridAvailabilityUse =
+        CoordinateOperationContext::GridAvailabilityUse::USE_FOR_SORTING;
+    bool allowPivots = true;
+    std::vector<std::pair<std::string, std::string>> pivots;
+    bool usePROJGridAlternatives = true;
+    std::string mainDBPath;
+    std::vector<std::string> auxDBPath;
+    bool guessDialect = false;
+    std::string authority;
+    bool identify = false;
+    bool showSuperseded = false;
+
+    for (int i = 1; i < argc; i++) {
+        std::string arg(argv[i]);
+        if (arg == "-o" && i + 1 < argc) {
+            outputSwithSpecified = true;
+            i++;
+            auto formats(split(argv[i], ','));
+            for (auto format : formats) {
+                if (ci_equal(format, "all")) {
+                    outputOpt.PROJ5 = true;
+                    outputOpt.PROJ4 = true;
+                    outputOpt.WKT2_2018 = true;
+                    outputOpt.WKT2_2015 = true;
+                    outputOpt.WKT1_GDAL = true;
+                    outputOpt.WKT1_ESRI = true;
+                } else if (ci_equal(format, "default")) {
+                    outputOpt.PROJ5 = true;
+                    outputOpt.PROJ4 = false;
+                    outputOpt.WKT2_2018 = false;
+                    outputOpt.WKT2_2015 = true;
+                    outputOpt.WKT1_GDAL = false;
+                } else if (ci_equal(format, "PROJ4") ||
+                           ci_equal(format, "PROJ.4")) {
+                    outputOpt.PROJ4 = true;
+                } else if (ci_equal(format, "-PROJ4") ||
+                           ci_equal(format, "-PROJ.4")) {
+                    outputOpt.PROJ4 = false;
+                } else if (ci_equal(format, "PROJ")) {
+                    outputOpt.PROJ5 = true;
+                } else if (ci_equal(format, "-PROJ")) {
+                    outputOpt.PROJ5 = false;
+                } else if (ci_equal(format, "WKT_ALL") ||
+                           ci_equal(format, "WKT-ALL")) {
+                    outputOpt.WKT2_2018 = true;
+                    outputOpt.WKT2_2015 = true;
+                    outputOpt.WKT1_GDAL = true;
+                } else if (ci_equal(format, "-WKT_ALL") ||
+                           ci_equal(format, "-WKT-ALL")) {
+                    outputOpt.WKT2_2018 = false;
+                    outputOpt.WKT2_2015 = false;
+                    outputOpt.WKT1_GDAL = false;
+                } else if (ci_equal(format, "WKT2_2018") ||
+                           ci_equal(format, "WKT2-2018") ||
+                           ci_equal(format, "WKT2:2018")) {
+                    outputOpt.WKT2_2018 = true;
+                } else if (ci_equal(format, "WKT2_2018_SIMPLIFIED") ||
+                           ci_equal(format, "WKT2-2018_SIMPLIFIED") ||
+                           ci_equal(format, "WKT2:2018_SIMPLIFIED")) {
+                    outputOpt.WKT2_2018_SIMPLIFIED = true;
+                } else if (ci_equal(format, "-WKT2_2018") ||
+                           ci_equal(format, "-WKT2-2018") ||
+                           ci_equal(format, "-WKT2:2018")) {
+                    outputOpt.WKT2_2018 = false;
+                } else if (ci_equal(format, "WKT2_2015") ||
+                           ci_equal(format, "WKT2-2015") ||
+                           ci_equal(format, "WKT2:2015")) {
+                    outputOpt.WKT2_2015 = true;
+                } else if (ci_equal(format, "WKT2_2015_SIMPLIFIED") ||
+                           ci_equal(format, "WKT2-2015_SIMPLIFIED") ||
+                           ci_equal(format, "WKT2:2015_SIMPLIFIED")) {
+                    outputOpt.WKT2_2015_SIMPLIFIED = true;
+                } else if (ci_equal(format, "-WKT2_2015") ||
+                           ci_equal(format, "-WKT2-2015") ||
+                           ci_equal(format, "-WKT2:2015")) {
+                    outputOpt.WKT2_2015 = false;
+                } else if (ci_equal(format, "WKT1_GDAL") ||
+                           ci_equal(format, "WKT1-GDAL") ||
+                           ci_equal(format, "WKT1:GDAL")) {
+                    outputOpt.WKT1_GDAL = true;
+                } else if (ci_equal(format, "-WKT1_GDAL") ||
+                           ci_equal(format, "-WKT1-GDAL") ||
+                           ci_equal(format, "-WKT1:GDAL")) {
+                    outputOpt.WKT1_GDAL = false;
+                } else if (ci_equal(format, "WKT1_ESRI") ||
+                           ci_equal(format, "WKT1-ESRI") ||
+                           ci_equal(format, "WKT1:ESRI")) {
+                    outputOpt.WKT1_ESRI = true;
+                } else if (ci_equal(format, "-WKT1_ESRI") ||
+                           ci_equal(format, "-WKT1-ESRI") ||
+                           ci_equal(format, "-WKT1:ESRI")) {
+                    outputOpt.WKT1_ESRI = false;
+                } else {
+                    std::cerr << "Unrecognized value for option -o: " << format
+                              << std::endl;
+                    usage();
+                }
+            }
+        } else if (arg == "--bbox" && i + 1 < argc) {
+            i++;
+            auto bboxStr(argv[i]);
+            auto bbox(split(bboxStr, ','));
+            if (bbox.size() != 4) {
+                std::cerr << "Incorrect number of values for option --bbox: "
+                          << bboxStr << std::endl;
+                usage();
+            }
+            try {
+                bboxFilter = Extent::createFromBBOX(
+                                 c_locale_stod(bbox[0]), c_locale_stod(bbox[1]),
+                                 c_locale_stod(bbox[2]), c_locale_stod(bbox[3]))
+                                 .as_nullable();
+            } catch (const std::exception &e) {
+                std::cerr << "Invalid value for option --bbox: " << bboxStr
+                          << ", " << e.what() << std::endl;
+                usage();
+            }
+        } else if (arg == "--area" && i + 1 < argc) {
+            i++;
+            area = argv[i];
+        } else if (arg == "-k" && i + 1 < argc) {
+            i++;
+            std::string kind(argv[i]);
+            if (ci_equal(kind, "crs") || ci_equal(kind, "srs")) {
+                kindIsCRS = true;
+            } else if (ci_equal(kind, "operation")) {
+                kindIsCRS = false;
+            } else {
+                std::cerr << "Unrecognized value for option -k: " << kind
+                          << std::endl;
+                usage();
+            }
+        } else if ((arg == "-s" || arg == "--ssource-crs") && i + 1 < argc) {
+            i++;
+            sourceCRSStr = argv[i];
+        } else if ((arg == "-t" || arg == "--target-crs") && i + 1 < argc) {
+            i++;
+            targetCRSStr = argv[i];
+        } else if (arg == "-q" || arg == "--quiet") {
+            outputOpt.quiet = true;
+        } else if (arg == "--c-ify") {
+            outputOpt.c_ify = true;
+        } else if (arg == "--single-line") {
+            outputOpt.singleLine = true;
+        } else if (arg == "--summary") {
+            summary = true;
+        } else if (ci_equal(arg, "--boundcrs-to-wgs84")) {
+            buildBoundCRSToWGS84 = true;
+
+            // undocumented: only for debugging purposes
+        } else if (ci_equal(arg, "--no-proj-grid-alternatives")) {
+            usePROJGridAlternatives = false;
+
+        } else if (arg == "--spatial-test" && i + 1 < argc) {
+            i++;
+            std::string value(argv[i]);
+            if (ci_equal(value, "contains")) {
+                spatialCriterion = CoordinateOperationContext::
+                    SpatialCriterion::STRICT_CONTAINMENT;
+            } else if (ci_equal(value, "intersects")) {
+                spatialCriterion = CoordinateOperationContext::
+                    SpatialCriterion::PARTIAL_INTERSECTION;
+            } else {
+                std::cerr << "Unrecognized value for option --spatial-test: "
+                          << value << std::endl;
+                usage();
+            }
+        } else if (arg == "--crs-extent-use" && i + 1 < argc) {
+            i++;
+            std::string value(argv[i]);
+            if (ci_equal(value, "NONE")) {
+                crsExtentUse =
+                    CoordinateOperationContext::SourceTargetCRSExtentUse::NONE;
+            } else if (ci_equal(value, "BOTH")) {
+                crsExtentUse =
+                    CoordinateOperationContext::SourceTargetCRSExtentUse::BOTH;
+            } else if (ci_equal(value, "INTERSECTION")) {
+                crsExtentUse = CoordinateOperationContext::
+                    SourceTargetCRSExtentUse::INTERSECTION;
+            } else if (ci_equal(value, "SMALLEST")) {
+                crsExtentUse = CoordinateOperationContext::
+                    SourceTargetCRSExtentUse::SMALLEST;
+            } else {
+                std::cerr << "Unrecognized value for option --crs-extent-use: "
+                          << value << std::endl;
+                usage();
+            }
+        } else if (arg == "--grid-check" && i + 1 < argc) {
+            i++;
+            std::string value(argv[i]);
+            if (ci_equal(value, "none")) {
+                gridAvailabilityUse = CoordinateOperationContext::
+                    GridAvailabilityUse::IGNORE_GRID_AVAILABILITY;
+            } else if (ci_equal(value, "discard_missing")) {
+                gridAvailabilityUse = CoordinateOperationContext::
+                    GridAvailabilityUse::DISCARD_OPERATION_IF_MISSING_GRID;
+            } else if (ci_equal(value, "sort")) {
+                gridAvailabilityUse = CoordinateOperationContext::
+                    GridAvailabilityUse::USE_FOR_SORTING;
+            } else {
+                std::cerr << "Unrecognized value for option --grid-check: "
+                          << value << std::endl;
+                usage();
+            }
+        } else if (arg == "--pivot-crs" && i + 1 < argc) {
+            i++;
+            auto value(argv[i]);
+            if (ci_equal(std::string(value), "none")) {
+                allowPivots = false;
+            } else {
+                auto splitValue(split(value, ','));
+                for (const auto &v : splitValue) {
+                    auto auth_code = split(v, ':');
+                    if (auth_code.size() != 2) {
+                        std::cerr
+                            << "Unrecognized value for option --grid-check: "
+                            << value << std::endl;
+                        usage();
+                    }
+                    pivots.emplace_back(
+                        std::make_pair(auth_code[0], auth_code[1]));
+                }
+            }
+        } else if (arg == "--main-db-path" && i + 1 < argc) {
+            i++;
+            mainDBPath = argv[i];
+        } else if (arg == "--aux-db-path" && i + 1 < argc) {
+            i++;
+            auxDBPath.push_back(argv[i]);
+        } else if (arg == "--guess-dialect") {
+            guessDialect = true;
+        } else if (arg == "--authority" && i + 1 < argc) {
+            i++;
+            authority = argv[i];
+        } else if (arg == "--identify") {
+            identify = true;
+        } else if (arg == "--show-superseded") {
+            showSuperseded = true;
+        } else if (arg == "-?" || arg == "--help") {
+            usage();
+        } else if (arg[0] == '-') {
+            std::cerr << "Unrecognized option: " << arg << std::endl;
+            usage();
+        } else {
+            if (!user_string_specified) {
+                user_string_specified = true;
+                user_string = arg;
+            } else {
+                std::cerr << "Too many parameters: " << arg << std::endl;
+                usage();
+            }
+        }
+    }
+
+    if (bboxFilter && !area.empty()) {
+        std::cerr << "ERROR: --bbox and --area are exclusive" << std::endl;
+        std::exit(1);
+    }
+
+    DatabaseContextPtr dbContext;
+    try {
+        dbContext =
+            DatabaseContext::create(mainDBPath, auxDBPath).as_nullable();
+    } catch (const std::exception &e) {
+        if (!mainDBPath.empty() || !auxDBPath.empty() || !area.empty()) {
+            std::cerr << "ERROR: Cannot create database connection: "
+                      << e.what() << std::endl;
+            std::exit(1);
+        }
+        std::cerr << "WARNING: Cannot create database connection: " << e.what()
+                  << std::endl;
+    }
+
+    if (!sourceCRSStr.empty() && targetCRSStr.empty()) {
+        std::cerr << "Source CRS specified, but missing target CRS"
+                  << std::endl;
+        usage();
+    } else if (sourceCRSStr.empty() && !targetCRSStr.empty()) {
+        std::cerr << "Target CRS specified, but missing source CRS"
+                  << std::endl;
+        usage();
+    } else if (!sourceCRSStr.empty() && !targetCRSStr.empty()) {
+        if (user_string_specified) {
+            std::cerr << "Unused extra value" << std::endl;
+            usage();
+        }
+    } else if (!user_string_specified) {
+        std::cerr << "Missing user string" << std::endl;
+        usage();
+    }
+
+    if (!outputSwithSpecified) {
+        outputOpt.PROJ5 = true;
+        outputOpt.WKT2_2015 = true;
+    }
+
+    if (outputOpt.quiet &&
+        (outputOpt.PROJ5 + outputOpt.PROJ4 + outputOpt.WKT2_2018 +
+         outputOpt.WKT2_2015 + outputOpt.WKT1_GDAL) != 1) {
+        std::cerr << "-q can only be used with a single output format"
+                  << std::endl;
+        usage();
+    }
+
+    if (!user_string.empty()) {
+        auto obj(buildObject(dbContext, user_string, kindIsCRS, "input string",
+                             buildBoundCRSToWGS84, allowPivots,
+                             outputOpt.quiet));
+        if (guessDialect) {
+            auto dialect = WKTParser().guessDialect(user_string);
+            std::cout << "Guessed WKT dialect: ";
+            if (dialect == WKTParser::WKTGuessedDialect::WKT2_2018) {
+                std::cout << "WKT2_2018";
+            } else if (dialect == WKTParser::WKTGuessedDialect::WKT2_2015) {
+                std::cout << "WKT2_2015";
+            } else if (dialect == WKTParser::WKTGuessedDialect::WKT1_GDAL) {
+                std::cout << "WKT1_GDAL";
+            } else if (dialect == WKTParser::WKTGuessedDialect::WKT1_ESRI) {
+                std::cout << "WKT1_ESRI";
+            } else {
+                std::cout << "Not WKT / unknown";
+            }
+            std::cout << std::endl;
+        }
+        outputObject(dbContext, obj, allowPivots, outputOpt);
+        if (identify) {
+            auto crs = dynamic_cast<CRS *>(obj.get());
+            if (crs) {
+                try {
+                    auto res = crs->identify(
+                        dbContext
+                            ? AuthorityFactory::create(NN_NO_CHECK(dbContext),
+                                                       authority)
+                                  .as_nullable()
+                            : nullptr);
+                    std::cout << std::endl;
+                    std::cout << "Identification match count: " << res.size()
+                              << std::endl;
+                    for (const auto &pair : res) {
+                        const auto &identifiedCRS = pair.first;
+                        const auto &ids = identifiedCRS->identifiers();
+                        if (!ids.empty()) {
+                            std::cout << *ids[0]->codeSpace() << ":"
+                                      << ids[0]->code() << ": " << pair.second
+                                      << " %" << std::endl;
+                        } else {
+                            auto boundCRS =
+                                dynamic_cast<BoundCRS *>(identifiedCRS.get());
+                            if (boundCRS &&
+                                !boundCRS->baseCRS()->identifiers().empty()) {
+                                const auto &idsBase =
+                                    boundCRS->baseCRS()->identifiers();
+                                std::cout << "BoundCRS of "
+                                          << *idsBase[0]->codeSpace() << ":"
+                                          << idsBase[0]->code() << ": "
+                                          << pair.second << " %" << std::endl;
+                            } else {
+                                std::cout
+                                    << "un-identifier CRS: " << pair.second
+                                    << " %" << std::endl;
+                            }
+                        }
+                    }
+                } catch (const std::exception &e) {
+                    std::cerr << "Identification failed: " << e.what()
+                              << std::endl;
+                }
+            }
+        }
+    } else {
+
+        if (!area.empty()) {
+            assert(dbContext);
+            try {
+                if (area.find(' ') == std::string::npos &&
+                    area.find(':') != std::string::npos) {
+                    auto tokens = split(area, ':');
+                    if (tokens.size() == 2) {
+                        const std::string &areaAuth = tokens[0];
+                        const std::string &areaCode = tokens[1];
+                        bboxFilter = AuthorityFactory::create(
+                                         NN_NO_CHECK(dbContext), areaAuth)
+                                         ->createExtent(areaCode)
+                                         .as_nullable();
+                    }
+                }
+                if (!bboxFilter) {
+                    auto authFactory = AuthorityFactory::create(
+                        NN_NO_CHECK(dbContext), std::string());
+                    auto res = authFactory->listAreaOfUseFromName(area, false);
+                    if (res.size() == 1) {
+                        bboxFilter =
+                            AuthorityFactory::create(NN_NO_CHECK(dbContext),
+                                                     res.front().first)
+                                ->createExtent(res.front().second)
+                                .as_nullable();
+                    } else {
+                        res = authFactory->listAreaOfUseFromName(area, true);
+                        if (res.size() == 1) {
+                            bboxFilter =
+                                AuthorityFactory::create(NN_NO_CHECK(dbContext),
+                                                         res.front().first)
+                                    ->createExtent(res.front().second)
+                                    .as_nullable();
+                        } else if (res.empty()) {
+                            std::cerr << "No area of use matching provided name"
+                                      << std::endl;
+                            std::exit(1);
+                        } else {
+                            std::cerr << "Several candidates area of use "
+                                         "matching provided name :"
+                                      << std::endl;
+                            for (const auto &candidate : res) {
+                                auto obj =
+                                    AuthorityFactory::create(
+                                        NN_NO_CHECK(dbContext), candidate.first)
+                                        ->createExtent(candidate.second);
+                                std::cerr << "  " << candidate.first << ":"
+                                          << candidate.second << " : "
+                                          << *obj->description() << std::endl;
+                            }
+                            std::exit(1);
+                        }
+                    }
+                }
+            } catch (const std::exception &e) {
+                std::cerr << "Area of use retrieval failed: " << e.what()
+                          << std::endl;
+                std::exit(1);
+            }
+        }
+
+        outputOperations(
+            dbContext, sourceCRSStr, targetCRSStr, bboxFilter, spatialCriterion,
+            crsExtentUse, gridAvailabilityUse, allowPivots, pivots, authority,
+            usePROJGridAlternatives, showSuperseded, outputOpt, summary);
+    }
+
+    return 0;
+}
+
+//! @endcond