1 /***********************************************************************
2
3 The cct 4D Transformation program
4
5 ************************************************************************
6
7 cct is a 4D equivalent to the "proj" projection program.
8
9 cct is an acronym meaning "Coordinate Conversion and Transformation".
10
11 The acronym refers to definitions given in the OGC 08-015r2/ISO-19111
12 standard "Geographical Information -- Spatial Referencing by Coordinates",
13 which defines two different classes of coordinate operations:
14
15 *Coordinate Conversions*, which are coordinate operations where input
16 and output datum are identical (e.g. conversion from geographical to
17 cartesian coordinates) and
18
19 *Coordinate Transformations*, which are coordinate operations where
20 input and output datums differ (e.g. change of reference frame).
21
22 cct, however, also refers to Carl Christian Tscherning (1942--2014),
23 professor of Geodesy at the University of Copenhagen, mentor and advisor
24 for a generation of Danish geodesists, colleague and collaborator for
25 two generations of global geodesists, Secretary General for the
26 International Association of Geodesy, IAG (1995--2007), fellow of the
27 American Geophysical Union (1991), recipient of the IAG Levallois Medal
28 (2007), the European Geosciences Union Vening Meinesz Medal (2008), and
29 of numerous other honours.
30
31 cct, or Christian, as he was known to most of us, was recognized for his
32 good mood, his sharp wit, his tireless work, and his great commitment to
33 the development of geodesy - both through his scientific contributions,
34 comprising more than 250 publications, and by his mentoring and teaching
35 of the next generations of geodesists.
36
37 As Christian was an avid Fortran programmer, and a keen Unix connoisseur,
38 he would have enjoyed to know that his initials would be used to name a
39 modest Unix style transformation filter, hinting at the tireless aspect
40 of his personality, which was certainly one of the reasons he accomplished
41 so much, and meant so much to so many people.
42
43 Hence, in honour of cct (the geodesist) this is cct (the program).
44
45 ************************************************************************
46
47 Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26
48
49 ************************************************************************
50
51 * Copyright (c) 2016, 2017 Thomas Knudsen
52 * Copyright (c) 2017, SDFE
53 *
54 * Permission is hereby granted, free of charge, to any person obtaining a
55 * copy of this software and associated documentation files (the "Software"),
56 * to deal in the Software without restriction, including without limitation
57 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
58 * and/or sell copies of the Software, and to permit persons to whom the
59 * Software is furnished to do so, subject to the following conditions:
60 *
61 * The above copyright notice and this permission notice shall be included
62 * in all copies or substantial portions of the Software.
63 *
64 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
65 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
66 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
67 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
68 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
69 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
70 * DEALINGS IN THE SOFTWARE.
71
72 ***********************************************************************/
73
74 #include <ctype.h>
75 #include <math.h>
76 #include <stdio.h>
77 #include <stdlib.h>
78 #include <string.h>
79 #include <stdarg.h>
80
81 #include "proj.h"
82 #include "proj_internal.h"
83 #include "proj_strtod.h"
84 #include "optargpm.h"
85
86
87 static void logger(void *data, int level, const char *msg);
88 static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...);
89
90 /* Prototypes from functions in this file */
91 char *column (char *buf, int n);
92 PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time);
93
94
95 static const char usage[] = {
96 "--------------------------------------------------------------------------------\n"
97 "Usage: %s [-options]... [+operator_specs]... infile...\n"
98 "--------------------------------------------------------------------------------\n"
99 "Options:\n"
100 "--------------------------------------------------------------------------------\n"
101 " -c x,y,z,t Specify input columns for (up to) 4 input parameters.\n"
102 " Defaults to 1,2,3,4\n"
103 " -d n Specify number of decimals in output.\n"
104 " -I Do the inverse transformation\n"
105 " -o /path/to/file Specify output file name\n"
106 " -t value Provide a fixed t value for all input data (e.g. -t 0)\n"
107 " -z value Provide a fixed z value for all input data (e.g. -z 0)\n"
108 " -s n Skip n first lines of a infile\n"
109 " -v Verbose: Provide non-essential informational output.\n"
110 " Repeat -v for more verbosity (e.g. -vv)\n"
111 "--------------------------------------------------------------------------------\n"
112 "Long Options:\n"
113 "--------------------------------------------------------------------------------\n"
114 " --output Alias for -o\n"
115 " --columns Alias for -c\n"
116 " --decimals Alias for -d\n"
117 " --height Alias for -z\n"
118 " --time Alias for -t\n"
119 " --verbose Alias for -v\n"
120 " --inverse Alias for -I\n"
121 " --skip-lines Alias for -s\n"
122 " --help Alias for -h\n"
123 " --version Print version number\n"
124 "--------------------------------------------------------------------------------\n"
125 "Operator Specs:\n"
126 "--------------------------------------------------------------------------------\n"
127 "The operator specs describe the action to be performed by cct, e.g:\n"
128 "\n"
129 " +proj=utm +ellps=GRS80 +zone=32\n"
130 "\n"
131 "instructs cct to convert input data to Universal Transverse Mercator, zone 32\n"
132 "coordinates, based on the GRS80 ellipsoid.\n"
133 "\n"
134 "Hence, the command\n"
135 "\n"
136 " echo 12 55 | cct -z0 -t0 +proj=utm +zone=32 +ellps=GRS80\n"
137 "\n"
138 "Should give results comparable to the classic proj command\n"
139 "\n"
140 " echo 12 55 | proj +proj=utm +zone=32 +ellps=GRS80\n"
141 "--------------------------------------------------------------------------------\n"
142 "Examples:\n"
143 "--------------------------------------------------------------------------------\n"
144 "1. convert geographical input to UTM zone 32 on the GRS80 ellipsoid:\n"
145 " cct +proj=utm +ellps=GRS80 +zone=32\n"
146 "2. roundtrip accuracy check for the case above:\n"
147 " cct +proj=pipeline +proj=utm +ellps=GRS80 +zone=32 +step +step +inv\n"
148 "3. as (1) but specify input columns for longitude, latitude, height and time:\n"
149 " cct -c 5,2,1,4 +proj=utm +ellps=GRS80 +zone=32\n"
150 "4. as (1) but specify fixed height and time, hence needing only 2 cols in input:\n"
151 " cct -t 0 -z 0 +proj=utm +ellps=GRS80 +zone=32\n"
152 "--------------------------------------------------------------------------------\n"
153 };
154
155
logger(void * data,int level,const char * msg)156 static void logger(void *data, int level, const char *msg) {
157 FILE *stream;
158 int log_tell = proj_log_level(PJ_DEFAULT_CTX, PJ_LOG_TELL);
159
160 stream = (FILE *) data;
161
162 /* if we use PJ_LOG_NONE we always want to print stuff to stream */
163 if (level == PJ_LOG_NONE) {
164 fprintf(stream, "%s\n", msg);
165 return;
166 }
167
168 /* otherwise only print if log level set by user is high enough or error */
169 if (level <= log_tell || level == PJ_LOG_ERROR)
170 fprintf(stderr, "%s\n", msg);
171 }
172
173 FILE *fout;
174
print(PJ_LOG_LEVEL log_level,const char * fmt,...)175 static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...) {
176
177 va_list args;
178 char *msg_buf;
179
180 va_start( args, fmt );
181
182 msg_buf = (char *) malloc(100000);
183 if( msg_buf == nullptr ) {
184 va_end( args );
185 return;
186 }
187
188 vsprintf( msg_buf, fmt, args );
189
190 logger((void *) fout, log_level, msg_buf);
191
192 va_end( args );
193 free( msg_buf );
194 }
195
196
main(int argc,char ** argv)197 int main(int argc, char **argv) {
198 PJ *P;
199 PJ_COORD point;
200 PJ_PROJ_INFO info;
201 OPTARGS *o;
202 char blank_comment[] = "";
203 char whitespace[] = " ";
204 char *comment;
205 char *comment_delimiter;
206 char *buf;
207 int i, nfields = 4, skip_lines = 0, verbose;
208 double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL;
209 int decimals_angles = 10;
210 int decimals_distances = 4;
211 int columns_xyzt[] = {1, 2, 3, 4};
212 const char *longflags[] = {"v=verbose", "h=help", "I=inverse", "version", nullptr};
213 const char *longkeys[] = {
214 "o=output",
215 "c=columns",
216 "d=decimals",
217 "z=height",
218 "t=time",
219 "s=skip-lines",
220 nullptr};
221
222 fout = stdout;
223
224 /* coverity[tainted_data] */
225 o = opt_parse (argc, argv, "hvI", "cdozts", longflags, longkeys);
226 if (nullptr==o)
227 return 0;
228
229 if (opt_given (o, "h") || argc==1) {
230 printf (usage, o->progname);
231 return 0;
232 }
233
234 PJ_DIRECTION direction = opt_given (o, "I")? PJ_INV: PJ_FWD;
235
236 verbose = MIN(opt_given (o, "v"), 3); /* log level can't be larger than 3 */
237 if( verbose > 0 ) {
238 proj_log_level (PJ_DEFAULT_CTX, static_cast<PJ_LOG_LEVEL>(verbose));
239 }
240 proj_log_func (PJ_DEFAULT_CTX, (void *) fout, logger);
241
242 if (opt_given (o, "version")) {
243 print (PJ_LOG_NONE, "%s: %s", o->progname, pj_get_release ());
244 return 0;
245 }
246
247 if (opt_given (o, "o"))
248 fout = fopen (opt_arg (o, "output"), "wt");
249 if (nullptr==fout) {
250 print (PJ_LOG_ERROR, "%s: Cannot open '%s' for output", o->progname, opt_arg (o, "output"));
251 free (o);
252 return 1;
253 }
254
255 print (PJ_LOG_TRACE, "%s: Running in very verbose mode", o->progname);
256
257 if (opt_given (o, "z")) {
258 fixed_z = proj_atof (opt_arg (o, "z"));
259 nfields--;
260 }
261
262 if (opt_given (o, "t")) {
263 fixed_time = proj_atof (opt_arg (o, "t"));
264 nfields--;
265 }
266
267 if (opt_given (o, "d")) {
268 int dec = atoi (opt_arg (o, "d"));
269 decimals_angles = dec;
270 decimals_distances = dec;
271 }
272
273 if (opt_given (o, "s")) {
274 skip_lines = atoi (opt_arg(o, "s"));
275 }
276
277 if (opt_given (o, "c")) {
278 int ncols;
279 /* reset column numbers to ease comment output later on */
280 for (i=0; i<4; i++)
281 columns_xyzt[i] = 0;
282
283 /* cppcheck-suppress invalidscanf */
284 ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3);
285 if (ncols != nfields) {
286 print (PJ_LOG_ERROR, "%s: Too few input columns given: '%s'", o->progname, opt_arg (o, "c"));
287 free (o);
288 if (stdout != fout)
289 fclose (fout);
290 return 1;
291 }
292 }
293
294 /* Setup transformation */
295 P = proj_create_argv (nullptr, o->pargc, o->pargv);
296 if ((nullptr==P) || (0==o->pargc)) {
297 print (PJ_LOG_ERROR, "%s: Bad transformation arguments - (%s)\n '%s -h' for help",
298 o->progname, pj_strerrno (proj_errno(P)), o->progname);
299 free (o);
300 if (stdout != fout)
301 fclose (fout);
302 return 1;
303 }
304
305 info = proj_pj_info (P);
306 print (PJ_LOG_TRACE, "Final: %s argc=%d pargc=%d", info.definition, argc, o->pargc);
307
308 if (direction== PJ_INV) {
309 /* fail if an inverse operation is not available */
310 if (!info.has_inverse) {
311 print (PJ_LOG_ERROR, "Inverse operation not available");
312 if (stdout != fout)
313 fclose (fout);
314 return 1;
315 }
316 /* We have no API call for inverting an operation, so we brute force it. */
317 P->inverted = !(P->inverted);
318 }
319 direction = PJ_FWD;
320
321 /* Allocate input buffer */
322 buf = static_cast<char*>(calloc (1, 10000));
323 if (nullptr==buf) {
324 print (PJ_LOG_ERROR, "%s: Out of memory", o->progname);
325 pj_free (P);
326 free (o);
327 if (stdout != fout)
328 fclose (fout);
329 return 1;
330 }
331
332
333 /* Loop over all records of all input files */
334 while (opt_input_loop (o, optargs_file_format_text)) {
335 int err;
336 void *ret = fgets (buf, 10000, o->input);
337 char *c = column (buf, 1);
338 opt_eof_handler (o);
339 if (nullptr==ret) {
340 print (PJ_LOG_ERROR, "Read error in record %d", (int) o->record_index);
341 continue;
342 }
343 point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time);
344 if (skip_lines > 0) {
345 skip_lines--;
346 continue;
347 }
348
349 /* if it's a comment or blank line, we reflect it */
350 if (c && ((*c=='\0') || (*c=='#'))) {
351 fprintf (fout, "%s", buf);
352 continue;
353 }
354
355 if (HUGE_VAL==point.xyzt.x) {
356 /* otherwise, it must be a syntax error */
357 print (PJ_LOG_NONE, "# Record %d UNREADABLE: %s", (int) o->record_index, buf);
358 print (PJ_LOG_ERROR, "%s: Could not parse file '%s' line %d", o->progname, opt_filename (o), opt_record (o));
359 continue;
360 }
361
362 if (proj_angular_input (P, direction)) {
363 point.lpzt.lam = proj_torad (point.lpzt.lam);
364 point.lpzt.phi = proj_torad (point.lpzt.phi);
365 }
366 err = proj_errno_reset (P);
367 /* coverity[returned_value] */
368 point = proj_trans (P, direction, point);
369
370 if (HUGE_VAL==point.xyzt.x) {
371 /* transformation error */
372 print (PJ_LOG_NONE, "# Record %d TRANSFORMATION ERROR: %s (%s)",
373 (int) o->record_index, buf, pj_strerrno (proj_errno(P)));
374 proj_errno_restore (P, err);
375 continue;
376 }
377 proj_errno_restore (P, err);
378
379 /* handle comment string */
380 comment = column(buf, nfields+1);
381 if (opt_given(o, "c")) {
382 /* what number is the last coordinate column in the input data? */
383 int colmax = 0;
384 for (i=0; i<4; i++)
385 colmax = MAX(colmax, columns_xyzt[i]);
386 comment = column(buf, colmax+1);
387 }
388 /* remove the line feed from comment, as logger() above, invoked
389 by print() below (output), will add one */
390 size_t len = strlen(comment);
391 if (len >= 1)
392 comment[len - 1] = '\0';
393 comment_delimiter = (comment && *comment) ? whitespace : blank_comment;
394
395 /* Time to print the result */
396 /* use same arguments to printf format string for both radians and
397 degrees; convert radians to degrees before printing */
398 if (proj_angular_output (P, direction) ||
399 proj_degree_output (P, direction)) {
400 if (proj_angular_output (P, direction)) {
401 point.lpzt.lam = proj_todeg (point.lpzt.lam);
402 point.lpzt.phi = proj_todeg (point.lpzt.phi);
403 }
404 print (PJ_LOG_NONE, "%14.*f %14.*f %12.*f %12.4f%s%s",
405 decimals_angles, point.xyzt.x,
406 decimals_angles, point.xyzt.y,
407 decimals_distances, point.xyzt.z,
408 point.xyzt.t, comment_delimiter, comment
409 );
410 }
411 else
412 print (PJ_LOG_NONE, "%13.*f %13.*f %12.*f %12.4f%s%s",
413 decimals_distances, point.xyzt.x,
414 decimals_distances, point.xyzt.y,
415 decimals_distances, point.xyzt.z,
416 point.xyzt.t, comment_delimiter, comment
417 );
418 if( fout == stdout )
419 fflush(stdout);
420 }
421
422 proj_destroy(P);
423
424 if (stdout != fout)
425 fclose (fout);
426 free (o);
427 free (buf);
428 return 0;
429 }
430
431
432
433
434
435 /* return a pointer to the n'th column of buf */
column(char * buf,int n)436 char *column (char *buf, int n) {
437 int i;
438 if (n <= 0)
439 return buf;
440 for (i = 0; i < n; i++) {
441 while (isspace(*buf))
442 buf++;
443 if (i == n - 1)
444 break;
445 while ((0 != *buf) && !isspace(*buf))
446 buf++;
447 }
448 return buf;
449 }
450
451 /* column to double */
cold(char * args,int col)452 static double cold (char *args, int col) {
453 char *endp;
454 char *target;
455 double d;
456 target = column (args, col);
457 d = proj_strtod (target, &endp);
458 if (endp==target)
459 return HUGE_VAL;
460 return d;
461 }
462
parse_input_line(char * buf,int * columns,double fixed_height,double fixed_time)463 PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time) {
464 PJ_COORD err = proj_coord (HUGE_VAL, HUGE_VAL, HUGE_VAL, HUGE_VAL);
465 PJ_COORD result = err;
466 int prev_errno = errno;
467 errno = 0;
468
469 result.xyzt.z = fixed_height;
470 result.xyzt.t = fixed_time;
471 result.xyzt.x = cold (buf, columns[0]);
472 result.xyzt.y = cold (buf, columns[1]);
473 if (result.xyzt.z==HUGE_VAL)
474 result.xyzt.z = cold (buf, columns[2]);
475 if (result.xyzt.t==HUGE_VAL)
476 result.xyzt.t = cold (buf, columns[3]);
477
478 if (0!=errno)
479 return err;
480
481 errno = prev_errno;
482 return result;
483 }
484