From e6de172371ea203f6393d745641d66c82b5b13e2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 19 Dec 2018 12:25:33 +0100 Subject: cpp conversion: move source files in apps/ iso19111/ conversions/ projections/ transformations/ tests/ subdirectories --- src/apps/cct.cpp | 472 ++++++++++++++ src/apps/cs2cs.cpp | 640 +++++++++++++++++++ src/apps/emess.cpp | 68 ++ src/apps/emess.h | 29 + src/apps/gen_cheb.cpp | 106 ++++ src/apps/geod.cpp | 241 +++++++ src/apps/geod_interface.cpp | 33 + src/apps/geod_interface.h | 45 ++ src/apps/geod_set.cpp | 75 +++ src/apps/gie.cpp | 1458 +++++++++++++++++++++++++++++++++++++++++++ src/apps/nad2bin.cpp | 382 ++++++++++++ src/apps/optargpm.h | 635 +++++++++++++++++++ src/apps/p_series.cpp | 42 ++ src/apps/proj.cpp | 582 +++++++++++++++++ src/apps/proj_strtod.cpp | 440 +++++++++++++ src/apps/proj_strtod.h | 4 + src/apps/projinfo.cpp | 1067 +++++++++++++++++++++++++++++++ 17 files changed, 6319 insertions(+) create mode 100644 src/apps/cct.cpp create mode 100644 src/apps/cs2cs.cpp create mode 100644 src/apps/emess.cpp create mode 100644 src/apps/emess.h create mode 100644 src/apps/gen_cheb.cpp create mode 100644 src/apps/geod.cpp create mode 100644 src/apps/geod_interface.cpp create mode 100644 src/apps/geod_interface.h create mode 100644 src/apps/geod_set.cpp create mode 100644 src/apps/gie.cpp create mode 100644 src/apps/nad2bin.cpp create mode 100644 src/apps/optargpm.h create mode 100644 src/apps/p_series.cpp create mode 100644 src/apps/proj.cpp create mode 100644 src/apps/proj_strtod.cpp create mode 100644 src/apps/proj_strtod.h create mode 100644 src/apps/projinfo.cpp (limited to 'src/apps') 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 +#include +#include +#include +#include +#include + +#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(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(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 +#include +#include +#include +#include +#include + +#include +#include + +#include + +// 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]* " 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("-"); + 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("-"); + + /* + * 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(""); + + } 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 +#include +#include +#include +#include + +#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 = ""; +#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 +#include +#include +#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(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(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(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 +# include +# include + +# 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("-"); + 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("-"); + for ( ; eargc-- ; ++eargv) { + if (**eargv == '-') { + fid = stdin; + emess_dat.File_name = const_cast(""); + } 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 +#include +#include + +#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 +#include +#include +#include +#include +#include +#include +#include + +#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[] = { + "", "operation", "use_proj4_init_rules", + "accept", "expect", "roundtrip", "banner", "verbose", + "direction", "tolerance", "ignore", "require_grid", "echo", "skip", "" +}; + +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 '' cmnd - bye!\n", fname); + if (F->level && F->level%2) + return errmsg (-4, "File '%s':Missing '' 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 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 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(calloc (1, sizeof (ffio))); + if (nullptr==G) + return nullptr; + + if (0==max_record_size) + max_record_size = 1000; + + G->args = static_cast(calloc (1, 5*max_record_size)); + if (nullptr==G->args) { + free (G); + return nullptr; + } + + G->next_args = static_cast(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 , 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 -block. Return 1 on success, 0 otherwise. +****************************************************************************************/ + /* Already inside */ + if (G->level % 2) + return 1; + + if (0==locate_tag (G, "")) + return 0; + + while (0!=strncmp ("", 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 block, locate the next and retry */ + if (0==strcmp (c, "")) { + 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(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 +#include + +#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: + + + [short ("-")/long ("--") options} + [operator ("+") specs] + [operands/input files] + +or less verbose: + + [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