Introducing the cct 'Coordinate Conversion and Transformation' program (#574)

* Introducing the cct 'Coordinate Conversion and Transformation' program
* cct: Add some rudimentary documentation
* Removed documentation again, moving to a separate doc PR
* Minor corrections in response to a review by @kbevers

7 files changed, 1194 insertions, 3 deletions

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index eec7ddec..788273a9 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,6 +3,7 @@
 include(lib_proj.cmake)
 
 # configure executable build
+option(BUILD_CCT   "Build cct (coordinate conversion and transformation tool)" ON)
 option(BUILD_CS2CS "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON)
 option(BUILD_PROJ "Build proj (cartographic projection tool : latlong <-> projected coordinates" ON)
 option(BUILD_GEOD "Build geod (computation of geodesic lines)" ON)
@@ -22,6 +23,11 @@ if(NOT MSVC)
   endif ()
 endif ()
 
+if(BUILD_CCT)
+ include(bin_cct.cmake)
+ set(BIN_TARGETS ${BIN_TARGETS} cct)
+endif(BUILD_CCT)
+
 if(BUILD_CS2CS)
  include(bin_cs2cs.cmake)
  set(BIN_TARGETS ${BIN_TARGETS} cs2cs)
diff --git a/src/Makefile.am b/src/Makefile.am
index 12d11c5f..de8fff4b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,6 +1,6 @@
 AM_CFLAGS = @C_WFLAGS@
 
-bin_PROGRAMS =	proj nad2bin geod cs2cs
+bin_PROGRAMS =	proj nad2bin geod cs2cs cct
 EXTRA_PROGRAMS = multistresstest test228
 
 TESTS = geodtest
@@ -12,12 +12,13 @@ AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 include_HEADERS = proj.h proj_api.h projects.h geodesic.h \
 	org_proj4_Projections.h org_proj4_PJ.h
 
-EXTRA_DIST = makefile.vc proj.def bin_cs2cs.cmake \
+EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_cs2cs.cmake \
 			 bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake \
 			 lib_proj.cmake CMakeLists.txt bin_geodtest.cmake geodtest.c
 
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
+cct_SOURCES = cct.c proj_strtod.c optargpm.h
 nad2bin_SOURCES = nad2bin.c
 geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
 multistresstest_SOURCES = multistresstest.c
@@ -26,6 +27,7 @@ geodtest_SOURCES = geodtest.c
 
 proj_LDADD = libproj.la
 cs2cs_LDADD = libproj.la
+cct_LDADD = libproj.la
 nad2bin_LDADD = libproj.la
 geod_LDADD = libproj.la
 multistresstest_LDADD = libproj.la @THREAD_LIB@
diff --git a/src/bin_cct.cmake b/src/bin_cct.cmake
new file mode 100644
index 00000000..a204e7e7
--- /dev/null
+++ b/src/bin_cct.cmake
@@ -0,0 +1,9 @@
+set(CCT_SRC cct.c proj_strtod.c) 
+set(CCT_INCLUDE optargpm.h)
+
+source_group("Source Files\\Bin" FILES ${CCT_SRC})
+
+add_executable(cct ${CCT_SRC} ${CCT_INCLUDE})
+target_link_libraries(cct ${PROJ_LIBRARIES})
+install(TARGETS cct 
+        RUNTIME DESTINATION ${BINDIR})
diff --git a/src/cct.c b/src/cct.c
new file mode 100644
index 00000000..83a0b0a3
--- /dev/null
+++ b/src/cct.c
@@ -0,0 +1,324 @@
+/***********************************************************************
+
+                 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
+Amercan 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 connoiseur,
+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-09-19
+
+************************************************************************
+
+* 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 "optargpm.h"
+#include <proj.h>
+#include <projects.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <math.h>
+
+double proj_strtod(const char *str, char **endptr);
+double proj_atof(const char *str);
+
+char *column (char *buf, int n);
+PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time);
+int print_output_line (FILE *fout, char *buf, PJ_COORD point);
+int main(int argc, char **argv);
+    
+
+
+static const char usage[] = {
+    "--------------------------------------------------------------------------------\n"
+    "Usage: %s [-options]... [+operator_specs]... infile...\n"
+    "--------------------------------------------------------------------------------\n"
+    "Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    -o /path/to/file  Specify output file name\n"
+    "    -c x,y,z,t        Specify input columns for (up to) 4 input parameters.\n"
+    "                      Defaults to 1,2,3,4\n"
+    "    -z value          Provide a fixed z value for all input data (e.g. -z 0)\n"
+    "    -t value          Provide a fixed t value for all input data (e.g. -t 0)\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"
+    "    --height          Alias for -z\n"
+    "    --time            Alias for -t\n"
+    "    --verbose         Alias for -v\n"
+    "    --help            Alias for -h\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"    
+};
+
+int main(int argc, char **argv) {
+    PJ *P;
+    PJ_COORD point;
+    OPTARGS *o;
+    FILE *fout = stdout;
+    char *buf;
+    int input_unit, output_unit, nfields = 4, direction = 1, verbose;
+    double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL;
+    int columns_xyzt[] = {1, 2, 3, 4};
+    const char *longflags[]  = {"v=verbose", "h=help", 0};
+    const char *longkeys[]   = {"o=output",  "c=columns", "z=height", "t=time", 0};
+    
+    o = opt_parse (argc, argv, "hv", "cozt", longflags, longkeys);
+    if (0==o)
+        return 0;
+
+    if (opt_given (o, "h")) {
+        printf (usage, o->progname);
+        return 0;
+    }
+
+
+    direction = opt_given (o, "I")? -1: 1;
+    verbose   = opt_given (o, "v");
+
+    if (opt_given (o, "o"))
+        fout = fopen (opt_arg (o, "output"), "rt");
+    if (0==fout) {
+        fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
+        free (o);
+        return 1;
+    }
+    if (verbose > 3)
+        fprintf (fout, "%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, "c")) {
+        int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+3, columns_xyzt+3);
+        if (ncols != nfields) {
+            fprintf (stderr, "%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 (0, o->pargc, o->pargv);
+    if ((0==P) || (0==o->pargc)) {
+        fprintf (stderr, "%s: Bad transformation arguments. '%s -h' for help\n", o->progname, o->progname);
+        free (o);
+        if (stdout != fout)
+            fclose (fout);
+        return 1;
+    }
+
+    input_unit   =  P->left;
+    output_unit  =  P->right;
+    if (PJ_IO_UNITS_CLASSIC==P->left)
+        input_unit = PJ_IO_UNITS_RADIANS;
+    if (PJ_IO_UNITS_CLASSIC==P->right)
+        output_unit = PJ_IO_UNITS_METERS;
+    if (direction==-1) {
+        enum pj_io_units swap = input_unit;
+        input_unit = output_unit;
+        output_unit = swap;
+    }
+
+    /* Allocate input buffer */
+    buf = calloc (1, 10000);
+    if (0==buf) {
+        fprintf (stderr, "%s: Out of memory\n", o->progname);
+        pj_free (P);
+        free (o);
+        if (stdout != fout)
+            fclose (fout);
+        return 1;
+    }
+
+
+    /* Loop over all lines of all input files */
+    while (opt_input_loop (o, optargs_file_format_text)) {
+        void *ret = fgets (buf, 10000, o->input);
+        int res;
+        opt_eof_handler (o);
+        if (0==ret) {
+            fprintf (stderr, "Read error in record %d\n", (int) o->record_index);
+            continue;
+        }
+        point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time);
+        if (PJ_IO_UNITS_RADIANS==input_unit) {
+            point.lpzt.lam = proj_torad (point.lpzt.lam);
+            point.lpzt.phi = proj_torad (point.lpzt.phi);
+        }
+        point = proj_trans_coord (P, direction, point);
+        if (PJ_IO_UNITS_RADIANS==output_unit) {
+            point.lpzt.lam = proj_todeg (point.lpzt.lam);
+            point.lpzt.phi = proj_todeg (point.lpzt.phi);
+        }
+        res = print_output_line (fout, buf, point);
+        if (0==res) {
+            fprintf (fout, "# UNREADABLE: %s", buf);
+            if (verbose)
+                fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o));
+        }
+    }
+    if (stdout != fout)
+        fclose (fout);
+    free (o);
+    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;
+}
+
+
+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;
+    char *endptr = 0;
+    errno = 0;
+
+    result.xyzt.z = fixed_height;
+    result.xyzt.t = fixed_time;
+    result.xyzt.x = proj_strtod (column (buf, columns[0]), &endptr);
+    result.xyzt.y = proj_strtod (column (buf, columns[1]), &endptr);
+    if (result.xyzt.z==HUGE_VAL)
+        result.xyzt.z = proj_strtod (column (buf, columns[2]), &endptr);
+    if (result.xyzt.t==HUGE_VAL)
+        result.xyzt.t = proj_strtod (column (buf, columns[3]), &endptr);
+
+    if (0!=errno)
+        return err;
+
+    errno = prev_errno;
+    return result;
+}
+
+
+int print_output_line (FILE *fout, char *buf, PJ_COORD point) {
+    char *c;
+    if (HUGE_VAL!=point.xyzt.x)
+        return fprintf (fout, "%20.15f  %20.15f  %20.15f  %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
+    c = column (buf, 1);
+    /* reflect comments and blanks */
+    if (c && ((*c=='\0') || (*c=='#')))
+        return fprintf (fout, "%s\n", buf);
+    return 0;
+}
diff --git a/src/makefile.vc b/src/makefile.vc
index cc57d806..7ab0dfc7 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -70,11 +70,13 @@ LIBOBJ	=	$(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \
 PROJEXE_OBJ	= proj.obj gen_cheb.obj p_series.obj emess.obj
 CS2CSEXE_OBJ	= cs2cs.obj gen_cheb.obj p_series.obj emess.obj
 GEODEXE_OBJ	= geod.obj geod_set.obj geod_interface.obj emess.obj
+CCTEXE_OBJ = cct.obj proj_strtod.obj
 MULTISTRESSTEST_OBJ = multistresstest.obj
 PROJ_DLL 	= proj$(VERSION).dll
 PROJ_EXE    = proj.exe
 CS2CS_EXE   = cs2cs.exe
 GEOD_EXE    = geod.exe
+CCT_EXE     = cct.exe
 NAD2BIN_EXE = nad2bin.exe
 MULTISTRESSTEST_EXE = multistresstest.exe
 
@@ -83,7 +85,7 @@ CFLAGS	=	/nologo -I. -DPROJ_LIB=\"$(PROJ_LIB_DIR)\" \
 
 default:	all
 
-all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(NAD2BIN_EXE)
+all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(CCT_EXE) $(NAD2BIN_EXE)
 
 proj.lib:	$(LIBOBJ)
 	if exist proj.lib del proj.lib
@@ -108,6 +110,10 @@ $(GEOD_EXE):	$(GEODEXE_OBJ) $(EXE_PROJ)
 	cl $(GEODEXE_OBJ) $(EXE_PROJ)
 	if exist $(GEOD_EXE).manifest mt -manifest $(GEOD_EXE).manifest -outputresource:$(GEOD_EXE);1
 
+$(CCT_EXE):	$(CCTEXE_OBJ) $(EXE_PROJ)
+	cl $(CCTEXE_OBJ) $(EXE_PROJ)
+	if exist $(CCT_EXE).manifest mt -manifest $(CCT_EXE).manifest -outputresource:$(CCT_EXE);1
+
 $(NAD2BIN_EXE):	nad2bin.obj emess.obj $(EXE_PROJ)
 	cl nad2bin.obj emess.obj $(EXE_PROJ)
 
diff --git a/src/optargpm.h b/src/optargpm.h
new file mode 100644
index 00000000..6be2c9ef
--- /dev/null
+++ b/src/optargpm.h
@@ -0,0 +1,593 @@
+/***********************************************************************
+
+          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
+prefixwd by "--", transformation specification elements prefixed by "+",
+and input file names prefixed by "" 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 claiming neither 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 conventinal "-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 (0==o)
+        return 0;
+
+
+    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 (0==ret) {
+            fprintf (stderr, "Read error in record %d\n", (int) o->record_index);
+            continue;
+        }
+        do_what_needs_to_be_done (buf);
+    }
+
+    return 0;
+}
+
+
+-------------------------------------------------------------------------------
+
+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.
+
+***********************************************************************/
+
+#define PJ_LIB__
+#include <proj.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include <math.h>
+#include <errno.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, char *option);
+int opt_given (OPTARGS *opt, char *option);
+char *opt_arg (OPTARGS *opt, char *option);
+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;
+    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 (0==opt)
+        return 0;
+    if (0==opt->fargc)
+        return opt->flaglevel;
+    return opt->fargv[opt->input_index];
+}
+
+static int opt_eof (OPTARGS *opt) {
+    if (0==opt)
+        return 1;
+    return feof (opt->input);
+}
+
+/* record number of most recently read record */
+int opt_record (OPTARGS *opt) {
+    if (0==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 (0==opt)
+        return 0;
+
+    /* most common case: increment record index and read on */
+    if ( (opt->input!=0) && !feof (opt->input) ) {
+        opt->record_index++;
+        return 1;
+    }
+
+    opt->record_index = 0;
+        
+    /* no input files specified - read from stdin */
+    if ((0==opt->fargc) && (0==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 (0!=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");
+
+    /* ignore non-existing files - go on! */
+    if (0==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, char *option) {
+    int i;
+    if (0==opt)
+        return 0;
+    if (0==option)
+        return 0;
+    if (0==option[0])
+        return 0;
+    /* An ordinary -o style short option */
+    if (strlen (option)==1) {
+        /* Undefined option? */
+        if (0==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 (0==f)
+            break;
+        if (0==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 (0==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 (0==v)
+            return 0;
+        if (0==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 (0==opt->optarg[(int) v[i][0]])
+                return 0;
+            return (int) v[i][0];
+        }
+
+    }
+
+    return 0;
+}
+
+
+/* Returns 0 if option was not given on command line, non-0 otherwise */
+int opt_given (OPTARGS *opt, 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, char *option) {
+    int ordinal = opt_ordinal (opt, option);
+    if (0==ordinal)
+        return 0;
+    return opt->optarg[ordinal];
+}
+
+
+/* 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;
+    OPTARGS *o;
+    char *last_path_delim;
+    
+    o = (OPTARGS *) calloc (1, sizeof(OPTARGS));
+    if (0==o)
+        return 0;
+
+    o->argc = argc;
+    o->argv = argv;
+    o->progname = argv[0];
+    last_path_delim = strrchr (argv[0], '\\');
+    if (last_path_delim > o->progname)
+        o->progname = last_path_delim;
+    last_path_delim = strrchr (argv[0], '/');
+    if (last_path_delim > o->progname)
+        o->progname = last_path_delim;
+        
+
+    /* 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 (0==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 0;            
+        }
+    }
+    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 (0==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 0;            
+        }
+    }
+
+    /* 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 != 0) && (longflags[i - 128] != 0); i++) {
+        if (i==192) {
+            free (o);
+            fprintf (stderr, "Too many flag style long options\n");
+            return 0;
+        }
+        o->optarg[i] = o->flaglevel;
+    }
+
+    for (i = 192;  (longkeys != 0) && (longkeys[i - 192] != 0);  i++) {
+        if (i==256) {
+            free (o);
+            fprintf (stderr, "Too many value style long options\n");
+            return 0;
+        }
+        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 (0==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;
+
+                /* need to maniplulate a bit to support gnu style --pap=pop syntax */
+                equals = strchr (crepr, '=');
+                if (equals)
+                    *equals = 0;
+                c = opt_ordinal (o, crepr);
+                if (0==c)
+                    return fprintf (stderr, "Invalid option \"%s\"\n", crepr), (OPTARGS *) 0;
+
+                /* inline (gnu) --foo=bar style arg */
+                if (equals) {
+                    *equals = '=';
+                    if (opt_is_flag (o, c))
+                        return fprintf (stderr, "Option \"%s\" takes no arguments\n", crepr), (OPTARGS *) 0;
+                    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]))
+                        return fprintf (stderr, "Missing argument for option \"%s\"\n", crepr), (OPTARGS *) 0;
+                    o->optarg[c] = argv[i + 1];
+                    i++; /* eat the arg */
+                    break;
+                }
+
+                if (!opt_is_flag (o, c))
+                    return fprintf (stderr, "Expected flag style long option here, but got \"%s\"\n", crepr), (OPTARGS *) 0;
+
+                /* Flag style option, i.e. taking no arguments */
+                opt_raise_flag (o, c);
+                break;
+            }
+
+            /* classic short options */
+            if (0==o->optarg[c])
+                return fprintf (stderr, "Invalid option \"%s\"\n", crepr), (OPTARGS *) 0;
+
+            /* 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]))
+                    return fprintf (stderr, "Bad or missing arg for option \"%s\"\n", crepr), (OPTARGS *) 0;
+                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;
+    for (/* empty */; i < argc; i++) {
+        if ('-' == argv[i][0]) {
+            free (o);
+            fprintf (stderr, "+ and - style options must not be mixed\n");
+            return 0;
+        }
+
+        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/proj_strtod.c b/src/proj_strtod.c
new file mode 100644
index 00000000..e9942e5c
--- /dev/null
+++ b/src/proj_strtod.c
@@ -0,0 +1,251 @@
+/***********************************************************************
+
+ 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 <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 proj_atof(const char *str);
+    
+
+double proj_strtod(const char *str, char **endptr) {
+    double number = 0;
+    int exponent = 0;
+    int sign = 0;
+    char *p = (char *) str;
+    int n = 0;
+    int num_digits_total = 0;
+    int num_digits_after_comma = 0;
+
+    if (0==str) {
+        errno = EFAULT;
+        if (endptr)
+            *endptr = p;
+        return HUGE_VAL;
+    }
+
+    /* First skip leading whitespace */
+    while (isspace(*p))
+        p++;
+
+    /* Empty string? */
+    if (0==*p) {
+        errno = EINVAL;
+        if (endptr)
+            *endptr = p;
+        return HUGE_VAL;
+    }
+
+    /* Then handle optional prefixed sign */
+    switch (*p) {
+        case '-':
+            sign = -1,  p++;  break;
+        case '+':
+            sign =  1,  p++;  break;
+        default:
+            if (isdigit(*p) || '_'==*p || '.'==*p)
+                break;
+            if (endptr)
+                *endptr = p;
+            errno = EINVAL;
+            return HUGE_VAL;
+    }
+
+    /* 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++;
+    }
+        
+    /* Do we have a fractional part? */
+    if ('.'==*p) {
+        p++;
+
+        while (isdigit(*p) || '_'==*p) {
+            if ('_'==*p) {
+                p++;
+                continue;
+            }
+            number = number * 10. + (*p - '0');
+            p++;
+            num_digits_total++;
+            num_digits_after_comma++;
+        }
+
+        exponent = -num_digits_after_comma;
+    }
+
+    /* 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? */
+    if (*p == 'e' || *p == 'E') {
+        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;
+    }
+
+    if ((exponent < DBL_MIN_EXP) || (exponent > DBL_MAX_EXP)) {
+      errno = ERANGE;
+      if (endptr)
+          *endptr = p;
+      return HUGE_VAL;
+    }
+
+    number *= pow (10, exponent);
+
+    if (fabs(number) > DBL_MAX)
+        errno = ERANGE;
+
+    if (endptr)
+        *endptr = p;
+
+    return number;
+}
+
+double proj_atof(const char *str) {
+  return proj_strtod(str, (void *) 0);
+}
+
+#ifdef TEST
+#include <string.h>
+
+int main (int argc, char **argv) {
+    double res;
+    char *endptr;
+    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