From 33e1c8e78f7bb826617f08f29e182b530d2ea153 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 28 Sep 2017 16:40:27 +0200 Subject: Introducing the cct 'Coordinate Conversion and Transformation' program (#574) * Introducing the cct 'Coordinate Conversion and Transformation' program * cct: Add some rudimentary documentation * Removed documentation again, moving to a separate doc PR * Minor corrections in response to a review by @kbevers --- src/CMakeLists.txt | 6 + src/Makefile.am | 6 +- src/bin_cct.cmake | 9 + src/cct.c | 324 +++++++++++++++++++++++++++++ src/makefile.vc | 8 +- src/optargpm.h | 593 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/proj_strtod.c | 251 +++++++++++++++++++++++ 7 files changed, 1194 insertions(+), 3 deletions(-) create mode 100644 src/bin_cct.cmake create mode 100644 src/cct.c create mode 100644 src/optargpm.h create mode 100644 src/proj_strtod.c (limited to 'src') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eec7ddec..788273a9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,6 +3,7 @@ include(lib_proj.cmake) # configure executable build +option(BUILD_CCT "Build cct (coordinate conversion and transformation tool)" ON) option(BUILD_CS2CS "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON) option(BUILD_PROJ "Build proj (cartographic projection tool : latlong <-> projected coordinates" ON) option(BUILD_GEOD "Build geod (computation of geodesic lines)" ON) @@ -22,6 +23,11 @@ if(NOT MSVC) endif () endif () +if(BUILD_CCT) + include(bin_cct.cmake) + set(BIN_TARGETS ${BIN_TARGETS} cct) +endif(BUILD_CCT) + if(BUILD_CS2CS) include(bin_cs2cs.cmake) set(BIN_TARGETS ${BIN_TARGETS} cs2cs) diff --git a/src/Makefile.am b/src/Makefile.am index 12d11c5f..de8fff4b 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ AM_CFLAGS = @C_WFLAGS@ -bin_PROGRAMS = proj nad2bin geod cs2cs +bin_PROGRAMS = proj nad2bin geod cs2cs cct EXTRA_PROGRAMS = multistresstest test228 TESTS = geodtest @@ -12,12 +12,13 @@ AM_CPPFLAGS = -DPROJ_LIB=\"$(pkgdatadir)\" \ include_HEADERS = proj.h proj_api.h projects.h geodesic.h \ org_proj4_Projections.h org_proj4_PJ.h -EXTRA_DIST = makefile.vc proj.def bin_cs2cs.cmake \ +EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_cs2cs.cmake \ bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake \ lib_proj.cmake CMakeLists.txt bin_geodtest.cmake geodtest.c proj_SOURCES = proj.c gen_cheb.c p_series.c cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c +cct_SOURCES = cct.c proj_strtod.c optargpm.h nad2bin_SOURCES = nad2bin.c geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h multistresstest_SOURCES = multistresstest.c @@ -26,6 +27,7 @@ geodtest_SOURCES = geodtest.c proj_LDADD = libproj.la cs2cs_LDADD = libproj.la +cct_LDADD = libproj.la nad2bin_LDADD = libproj.la geod_LDADD = libproj.la multistresstest_LDADD = libproj.la @THREAD_LIB@ diff --git a/src/bin_cct.cmake b/src/bin_cct.cmake new file mode 100644 index 00000000..a204e7e7 --- /dev/null +++ b/src/bin_cct.cmake @@ -0,0 +1,9 @@ +set(CCT_SRC cct.c proj_strtod.c) +set(CCT_INCLUDE optargpm.h) + +source_group("Source Files\\Bin" FILES ${CCT_SRC}) + +add_executable(cct ${CCT_SRC} ${CCT_INCLUDE}) +target_link_libraries(cct ${PROJ_LIBRARIES}) +install(TARGETS cct + RUNTIME DESTINATION ${BINDIR}) diff --git a/src/cct.c b/src/cct.c new file mode 100644 index 00000000..83a0b0a3 --- /dev/null +++ b/src/cct.c @@ -0,0 +1,324 @@ +/*********************************************************************** + + The cct 4D Transformation program + +************************************************************************ + +cct is a 4D equivalent to the "proj" projection program. + +cct is an acronym meaning "Coordinate Conversion and Transformation". + +The acronym refers to definitions given in the OGC 08-015r2/ISO-19111 +standard "Geographical Information -- Spatial Referencing by Coordinates", +which defines two different classes of coordinate operations: + +*Coordinate Conversions*, which are coordinate operations where input +and output datum are identical (e.g. conversion from geographical to +cartesian coordinates) and + +*Coordinate Transformations*, which are coordinate operations where +input and output datums differ (e.g. change of reference frame). + +cct, however, also refers to Carl Christian Tscherning (1942--2014), +professor of Geodesy at the University of Copenhagen, mentor and advisor +for a generation of Danish geodesists, colleague and collaborator for +two generations of global geodesists, Secretary General for the +International Association of Geodesy, IAG (1995--2007), fellow of the +Amercan Geophysical Union (1991), recipient of the IAG Levallois Medal +(2007), the European Geosciences Union Vening Meinesz Medal (2008), and +of numerous other honours. + +cct, or Christian, as he was known to most of us, was recognized for his +good mood, his sharp wit, his tireless work, and his great commitment to +the development of geodesy - both through his scientific contributions, +comprising more than 250 publications, and by his mentoring and teaching +of the next generations of geodesists. + +As Christian was an avid Fortran programmer, and a keen Unix connoiseur, +he would have enjoyed to know that his initials would be used to name a +modest Unix style transformation filter, hinting at the tireless aspect +of his personality, which was certainly one of the reasons he accomplished +so much, and meant so much to so many people. + +Hence, in honour of cct (the geodesist) this is cct (the program). + +************************************************************************ + +Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-09-19 + +************************************************************************ + +* Copyright (c) 2016, 2017 Thomas Knudsen +* Copyright (c) 2017, SDFE +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +* DEALINGS IN THE SOFTWARE. + +***********************************************************************/ + +#include "optargpm.h" +#include +#include +#include +#include +#include +#include +#include + +double proj_strtod(const char *str, char **endptr); +double proj_atof(const char *str); + +char *column (char *buf, int n); +PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time); +int print_output_line (FILE *fout, char *buf, PJ_COORD point); +int main(int argc, char **argv); + + + +static const char usage[] = { + "--------------------------------------------------------------------------------\n" + "Usage: %s [-options]... [+operator_specs]... infile...\n" + "--------------------------------------------------------------------------------\n" + "Options:\n" + "--------------------------------------------------------------------------------\n" + " -o /path/to/file Specify output file name\n" + " -c x,y,z,t Specify input columns for (up to) 4 input parameters.\n" + " Defaults to 1,2,3,4\n" + " -z value Provide a fixed z value for all input data (e.g. -z 0)\n" + " -t value Provide a fixed t value for all input data (e.g. -t 0)\n" + " -v Verbose: Provide non-essential informational output.\n" + " Repeat -v for more verbosity (e.g. -vv)\n" + "--------------------------------------------------------------------------------\n" + "Long Options:\n" + "--------------------------------------------------------------------------------\n" + " --output Alias for -o\n" + " --columns Alias for -c\n" + " --height Alias for -z\n" + " --time Alias for -t\n" + " --verbose Alias for -v\n" + " --help Alias for -h\n" + "--------------------------------------------------------------------------------\n" + "Operator Specs:\n" + "--------------------------------------------------------------------------------\n" + "The operator specs describe the action to be performed by cct, e.g:\n" + "\n" + " +proj=utm +ellps=GRS80 +zone=32\n" + "\n" + "instructs cct to convert input data to Universal Transverse Mercator, zone 32\n" + "coordinates, based on the GRS80 ellipsoid.\n" + "\n" + "Hence, the command\n" + "\n" + " echo 12 55 | cct -z0 -t0 +proj=utm +zone=32 +ellps=GRS80\n" + "\n" + "Should give results comparable to the classic proj command\n" + "\n" + " echo 12 55 | proj +proj=utm +zone=32 +ellps=GRS80\n" + "--------------------------------------------------------------------------------\n" + "Examples:\n" + "--------------------------------------------------------------------------------\n" + "1. convert geographical input to UTM zone 32 on the GRS80 ellipsoid:\n" + " cct +proj=utm +ellps=GRS80 +zone=32\n" + "2. roundtrip accuracy check for the case above:\n" + " cct +proj=pipeline +proj=utm +ellps=GRS80 +zone=32 +step +step +inv\n" + "3. as (1) but specify input columns for longitude, latitude, height and time:\n" + " cct -c 5,2,1,4 +proj=utm +ellps=GRS80 +zone=32\n" + "4. as (1) but specify fixed height and time, hence needing only 2 cols in input:\n" + " cct -t 0 -z 0 +proj=utm +ellps=GRS80 +zone=32\n" + "--------------------------------------------------------------------------------\n" +}; + +int main(int argc, char **argv) { + PJ *P; + PJ_COORD point; + OPTARGS *o; + FILE *fout = stdout; + char *buf; + int input_unit, output_unit, nfields = 4, direction = 1, verbose; + double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL; + int columns_xyzt[] = {1, 2, 3, 4}; + const char *longflags[] = {"v=verbose", "h=help", 0}; + const char *longkeys[] = {"o=output", "c=columns", "z=height", "t=time", 0}; + + o = opt_parse (argc, argv, "hv", "cozt", longflags, longkeys); + if (0==o) + return 0; + + if (opt_given (o, "h")) { + printf (usage, o->progname); + return 0; + } + + + direction = opt_given (o, "I")? -1: 1; + verbose = opt_given (o, "v"); + + if (opt_given (o, "o")) + fout = fopen (opt_arg (o, "output"), "rt"); + if (0==fout) { + fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output")); + free (o); + return 1; + } + if (verbose > 3) + fprintf (fout, "%s: Running in very verbose mode\n", o->progname); + + + + if (opt_given (o, "z")) { + fixed_z = proj_atof (opt_arg (o, "z")); + nfields--; + } + + if (opt_given (o, "t")) { + fixed_time = proj_atof (opt_arg (o, "t")); + nfields--; + } + + if (opt_given (o, "c")) { + int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+3, columns_xyzt+3); + if (ncols != nfields) { + fprintf (stderr, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c")); + free (o); + if (stdout != fout) + fclose (fout); + return 1; + } + } + + /* Setup transformation */ + P = proj_create_argv (0, o->pargc, o->pargv); + if ((0==P) || (0==o->pargc)) { + fprintf (stderr, "%s: Bad transformation arguments. '%s -h' for help\n", o->progname, o->progname); + free (o); + if (stdout != fout) + fclose (fout); + return 1; + } + + input_unit = P->left; + output_unit = P->right; + if (PJ_IO_UNITS_CLASSIC==P->left) + input_unit = PJ_IO_UNITS_RADIANS; + if (PJ_IO_UNITS_CLASSIC==P->right) + output_unit = PJ_IO_UNITS_METERS; + if (direction==-1) { + enum pj_io_units swap = input_unit; + input_unit = output_unit; + output_unit = swap; + } + + /* Allocate input buffer */ + buf = calloc (1, 10000); + if (0==buf) { + fprintf (stderr, "%s: Out of memory\n", o->progname); + pj_free (P); + free (o); + if (stdout != fout) + fclose (fout); + return 1; + } + + + /* Loop over all lines of all input files */ + while (opt_input_loop (o, optargs_file_format_text)) { + void *ret = fgets (buf, 10000, o->input); + int res; + opt_eof_handler (o); + if (0==ret) { + fprintf (stderr, "Read error in record %d\n", (int) o->record_index); + continue; + } + point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time); + if (PJ_IO_UNITS_RADIANS==input_unit) { + point.lpzt.lam = proj_torad (point.lpzt.lam); + point.lpzt.phi = proj_torad (point.lpzt.phi); + } + point = proj_trans_coord (P, direction, point); + if (PJ_IO_UNITS_RADIANS==output_unit) { + point.lpzt.lam = proj_todeg (point.lpzt.lam); + point.lpzt.phi = proj_todeg (point.lpzt.phi); + } + res = print_output_line (fout, buf, point); + if (0==res) { + fprintf (fout, "# UNREADABLE: %s", buf); + if (verbose) + fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o)); + } + } + if (stdout != fout) + fclose (fout); + free (o); + return 0; +} + + + + + +/* return a pointer to the n'th column of buf */ +char *column (char *buf, int n) { + int i; + if (n <= 0) + return buf; + for (i = 0; i < n; i++) { + while (isspace(*buf)) + buf++; + if (i == n - 1) + break; + while ((0 != *buf) && !isspace(*buf)) + buf++; + } + return buf; +} + + +PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time) { + PJ_COORD err = proj_coord (HUGE_VAL, HUGE_VAL, HUGE_VAL, HUGE_VAL); + PJ_COORD result = err; + int prev_errno = errno; + char *endptr = 0; + errno = 0; + + result.xyzt.z = fixed_height; + result.xyzt.t = fixed_time; + result.xyzt.x = proj_strtod (column (buf, columns[0]), &endptr); + result.xyzt.y = proj_strtod (column (buf, columns[1]), &endptr); + if (result.xyzt.z==HUGE_VAL) + result.xyzt.z = proj_strtod (column (buf, columns[2]), &endptr); + if (result.xyzt.t==HUGE_VAL) + result.xyzt.t = proj_strtod (column (buf, columns[3]), &endptr); + + if (0!=errno) + return err; + + errno = prev_errno; + return result; +} + + +int print_output_line (FILE *fout, char *buf, PJ_COORD point) { + char *c; + if (HUGE_VAL!=point.xyzt.x) + return fprintf (fout, "%20.15f %20.15f %20.15f %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); + c = column (buf, 1); + /* reflect comments and blanks */ + if (c && ((*c=='\0') || (*c=='#'))) + return fprintf (fout, "%s\n", buf); + return 0; +} diff --git a/src/makefile.vc b/src/makefile.vc index cc57d806..7ab0dfc7 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -70,11 +70,13 @@ LIBOBJ = $(support) $(pseudo) $(azimuthal) $(conic) $(cylinder) $(misc) \ PROJEXE_OBJ = proj.obj gen_cheb.obj p_series.obj emess.obj CS2CSEXE_OBJ = cs2cs.obj gen_cheb.obj p_series.obj emess.obj GEODEXE_OBJ = geod.obj geod_set.obj geod_interface.obj emess.obj +CCTEXE_OBJ = cct.obj proj_strtod.obj MULTISTRESSTEST_OBJ = multistresstest.obj PROJ_DLL = proj$(VERSION).dll PROJ_EXE = proj.exe CS2CS_EXE = cs2cs.exe GEOD_EXE = geod.exe +CCT_EXE = cct.exe NAD2BIN_EXE = nad2bin.exe MULTISTRESSTEST_EXE = multistresstest.exe @@ -83,7 +85,7 @@ CFLAGS = /nologo -I. -DPROJ_LIB=\"$(PROJ_LIB_DIR)\" \ default: all -all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(NAD2BIN_EXE) +all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(CCT_EXE) $(NAD2BIN_EXE) proj.lib: $(LIBOBJ) if exist proj.lib del proj.lib @@ -108,6 +110,10 @@ $(GEOD_EXE): $(GEODEXE_OBJ) $(EXE_PROJ) cl $(GEODEXE_OBJ) $(EXE_PROJ) if exist $(GEOD_EXE).manifest mt -manifest $(GEOD_EXE).manifest -outputresource:$(GEOD_EXE);1 +$(CCT_EXE): $(CCTEXE_OBJ) $(EXE_PROJ) + cl $(CCTEXE_OBJ) $(EXE_PROJ) + if exist $(CCT_EXE).manifest mt -manifest $(CCT_EXE).manifest -outputresource:$(CCT_EXE);1 + $(NAD2BIN_EXE): nad2bin.obj emess.obj $(EXE_PROJ) cl nad2bin.obj emess.obj $(EXE_PROJ) diff --git a/src/optargpm.h b/src/optargpm.h new file mode 100644 index 00000000..6be2c9ef --- /dev/null +++ b/src/optargpm.h @@ -0,0 +1,593 @@ +/*********************************************************************** + + OPTARGPM - a header-only library for decoding + PROJ.4 style command line options + + Thomas Knudsen, 2017-09-10 + +************************************************************************ + +For PROJ.4 command line programs, we have a somewhat complex option +decoding situation, since we have to navigate in a cocktail of classic +single letter style options, prefixed by "-", GNU style long options +prefixwd by "--", transformation specification elements prefixed by "+", +and input file names prefixed by "" nothing. + +Hence, classic getopt.h style decoding does not cut the mustard, so +this is an attempt to catch up and chop the ketchup. + +Since optargpm (for "optarg plus minus") does not belong, in any +obvious way, in any systems development library, it is provided as +a "header only" library. + +While this is conventional in C++, it is frowned at in plain C. +But frown away - "header only" has its places, and this is one of +them. + +By convention, we expect a command line to consist of the following +elements: + + + [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 claiming neither to save the world, nor to hint at the "shape of +jazz to come", at least optargpm has shown useful in constructing cs2cs +style transformation filters. + +Supporting a wide range of option syntax, the getoptpm API is somewhat +quirky, but also compact, consisting of one data type, 3(+2) functions, +and one enumeration: + +OPTARGS + Housekeeping data type. An instance of OPTARGS is conventionally + called o or opt +opt_parse (opt, argc, argv ...): + The work horse: Define supported options; Split (argc, argv) + into groups (options, op specs, operands); Parse option + arguments. +opt_given (o, option): + The number of times