/Proj4/proj.c
C | 502 lines | 472 code | 19 blank | 11 comment | 192 complexity | 9af58c2e6db7ca82ee3006e321e3b20c MD5 | raw file
1/* <<<< Cartographic projection filter program >>>> */ 2#ifndef lint 3static const char SCCSID[]="@(#)proj.c 4.12 95/09/23 GIE REL"; 4#endif 5#include "projects.h" 6#include <stdio.h> 7#include <stdlib.h> 8#include <ctype.h> 9#include <string.h> 10#include <math.h> 11#include "emess.h" 12 13/* TK 1999-02-13 */ 14#if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__WIN32__) 15# include <fcntl.h> 16# include <io.h> 17# define SET_BINARY_MODE(file) setmode(fileno(file), O_BINARY) 18#else 19# define SET_BINARY_MODE(file) 20#endif 21/* ! TK 1999-02-13 */ 22 23#define MAX_LINE 1000 24#define MAX_PARGS 100 25#define PJ_INVERS(P) (P->inv ? 1 : 0) 26 static PJ 27*Proj; 28 static projUV 29(*proj)(projUV, PJ *); 30 static int 31reversein = 0, /* != 0 reverse input arguments */ 32reverseout = 0, /* != 0 reverse output arguments */ 33bin_in = 0, /* != 0 then binary input */ 34bin_out = 0, /* != 0 then binary output */ 35echoin = 0, /* echo input data to output line */ 36tag = '#', /* beginning of line tag character */ 37inverse = 0, /* != 0 then inverse projection */ 38prescale = 0, /* != 0 apply cartesian scale factor */ 39dofactors = 0, /* determine scale factors */ 40facs_bad = 0, /* return condition from pj_factors */ 41very_verby = 0, /* very verbose mode */ 42postscale = 0; 43 static char 44*cheby_str, /* string controlling Chebychev evaluation */ 45*oform = (char *)0, /* output format for x-y or decimal degrees */ 46*oterr = "*\t*", /* output line for unprojectable input */ 47*usage = 48"%s\nusage: %s [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]\n"; 49 static struct FACTORS 50facs; 51 static double 52(*informat)(const char *, char **), /* input data deformatter function */ 53fscale = 0.; /* cartesian scale factor */ 54 static projUV 55int_proj(projUV data) { 56 if (prescale) { data.u *= fscale; data.v *= fscale; } 57 data = (*proj)(data, Proj); 58 if (postscale && data.u != HUGE_VAL) 59 { data.u *= fscale; data.v *= fscale; } 60 return(data); 61} 62 static void /* file processing function */ 63process(FILE *fid) { 64 char line[MAX_LINE+3], *s, pline[40]; 65 projUV data; 66 67 for (;;) { 68 ++emess_dat.File_line; 69 if (bin_in) { /* binary input */ 70 if (fread(&data, sizeof(projUV), 1, fid) != 1) 71 break; 72 } else { /* ascii input */ 73 if (!(s = fgets(line, MAX_LINE, fid))) 74 break; 75 if (!strchr(s, '\n')) { /* overlong line */ 76 int c; 77 (void)strcat(s, "\n"); 78 /* gobble up to newline */ 79 while ((c = fgetc(fid)) != EOF && c != '\n') ; 80 } 81 if (*s == tag) { 82 if (!bin_out) 83 (void)fputs(line, stdout); 84 continue; 85 } 86 if (reversein) { 87 data.v = (*informat)(s, &s); 88 data.u = (*informat)(s, &s); 89 } else { 90 data.u = (*informat)(s, &s); 91 data.v = (*informat)(s, &s); 92 } 93 if (data.v == HUGE_VAL) 94 data.u = HUGE_VAL; 95 if (!*s && (s > line)) --s; /* assumed we gobbled \n */ 96 if (!bin_out && echoin) { 97 int t; 98 t = *s; 99 *s = '\0'; 100 (void)fputs(line, stdout); 101 *s = t; 102 putchar('\t'); 103 } 104 } 105 if (data.u != HUGE_VAL) { 106 if (prescale) { data.u *= fscale; data.v *= fscale; } 107 if (dofactors && !inverse) 108 facs_bad = pj_factors(data, Proj, 0., &facs); 109 data = (*proj)(data, Proj); 110 if (dofactors && inverse) 111 facs_bad = pj_factors(data, Proj, 0., &facs); 112 if (postscale && data.u != HUGE_VAL) 113 { data.u *= fscale; data.v *= fscale; } 114 } 115 if (bin_out) { /* binary output */ 116 (void)fwrite(&data, sizeof(projUV), 1, stdout); 117 continue; 118 } else if (data.u == HUGE_VAL) /* error output */ 119 (void)fputs(oterr, stdout); 120 else if (inverse && !oform) { /*ascii DMS output */ 121 if (reverseout) { 122 (void)fputs(rtodms(pline, data.v, 'N', 'S'), stdout); 123 putchar('\t'); 124 (void)fputs(rtodms(pline, data.u, 'E', 'W'), stdout); 125 } else { 126 (void)fputs(rtodms(pline, data.u, 'E', 'W'), stdout); 127 putchar('\t'); 128 (void)fputs(rtodms(pline, data.v, 'N', 'S'), stdout); 129 } 130 } else { /* x-y or decimal degree ascii output */ 131 if (inverse) { 132 data.v *= RAD_TO_DEG; 133 data.u *= RAD_TO_DEG; 134 } 135 if (reverseout) { 136 (void)printf(oform,data.v); putchar('\t'); 137 (void)printf(oform,data.u); 138 } else { 139 (void)printf(oform,data.u); putchar('\t'); 140 (void)printf(oform,data.v); 141 } 142 } 143 if (dofactors) /* print scale factor data */ 144 if (!facs_bad) 145 (void)printf("\t<%g %g %g %g %g %g>", 146 facs.h, facs.k, facs.s, 147 facs.omega * RAD_TO_DEG, facs.a, facs.b); 148 else 149 (void)fputs("\t<* * * * * *>", stdout); 150 (void)fputs(bin_in ? "\n" : s, stdout); 151 } 152} 153 static void /* file processing function --- verbosely */ 154vprocess(FILE *fid) { 155 char line[MAX_LINE+3], *s, pline[40]; 156 projUV dat_ll, dat_xy; 157 int linvers; 158 159 if (!oform) 160 oform = "%.3f"; 161 if (bin_in || bin_out) 162 emess(1,"binary I/O not available in -V option"); 163 for (;;) { 164 ++emess_dat.File_line; 165 if (!(s = fgets(line, MAX_LINE, fid))) 166 break; 167 if (!strchr(s, '\n')) { /* overlong line */ 168 int c; 169 (void)strcat(s, "\n"); 170 /* gobble up to newline */ 171 while ((c = fgetc(fid)) != EOF && c != '\n') ; 172 } 173 if (*s == tag) { /* pass on data */ 174 (void)fputs(s, stdout); 175 continue; 176 } 177 /* check to override default input mode */ 178 if (*s == 'I' || *s == 'i') { 179 linvers = 1; 180 ++s; 181 } else if (*s == 'I' || *s == 'i') { 182 linvers = 0; 183 ++s; 184 } else 185 linvers = inverse; 186 if (linvers) { 187 if (!PJ_INVERS(Proj)) { 188 emess(-1,"inverse for this projection not avail.\n"); 189 continue; 190 } 191 dat_xy.u = strtod(s, &s); 192 dat_xy.v = strtod(s, &s); 193 if (dat_xy.u == HUGE_VAL || dat_xy.v == HUGE_VAL) { 194 emess(-1,"lon-lat input conversion failure\n"); 195 continue; 196 } 197 if (prescale) { dat_xy.u *= fscale; dat_xy.v *= fscale; } 198 dat_ll = pj_inv(dat_xy, Proj); 199 } else { 200 dat_ll.u = dmstor(s, &s); 201 dat_ll.v = dmstor(s, &s); 202 if (dat_ll.u == HUGE_VAL || dat_ll.v == HUGE_VAL) { 203 emess(-1,"lon-lat input conversion failure\n"); 204 continue; 205 } 206 dat_xy = pj_fwd(dat_ll, Proj); 207 if (postscale) { dat_xy.u *= fscale; dat_xy.v *= fscale; } 208 } 209 if (pj_errno) { 210 emess(-1, pj_strerrno(pj_errno)); 211 continue; 212 } 213 if (!*s && (s > line)) --s; /* assumed we gobbled \n */ 214 if (pj_factors(dat_ll, Proj, 0., &facs)) { 215 emess(-1,"failed to conpute factors\n\n"); 216 continue; 217 } 218 if (*s != '\n') 219 (void)fputs(s, stdout); 220 (void)fputs("Longitude: ", stdout); 221 (void)fputs(rtodms(pline, dat_ll.u, 'E', 'W'), stdout); 222 (void)printf(" [ %.11g ]\n", dat_ll.u * RAD_TO_DEG); 223 (void)fputs("Latitude: ", stdout); 224 (void)fputs(rtodms(pline, dat_ll.v, 'N', 'S'), stdout); 225 (void)printf(" [ %.11g ]\n", dat_ll.v * RAD_TO_DEG); 226 (void)fputs("Easting (x): ", stdout); 227 (void)printf(oform, dat_xy.u); putchar('\n'); 228 (void)fputs("Northing (y): ", stdout); 229 (void)printf(oform, dat_xy.v); putchar('\n'); 230 (void)printf("Meridian scale (h)%c: %.8f ( %.4g %% error )\n", 231 facs.code & IS_ANAL_HK ? '*' : ' ', facs.h, (facs.h-1.)*100.); 232 (void)printf("Parallel scale (k)%c: %.8f ( %.4g %% error )\n", 233 facs.code & IS_ANAL_HK ? '*' : ' ', facs.k, (facs.k-1.)*100.); 234 (void)printf("Areal scale (s): %.8f ( %.4g %% error )\n", 235 facs.s, (facs.s-1.)*100.); 236 (void)printf("Angular distortion (w): %.3f\n", facs.omega * 237 RAD_TO_DEG); 238 (void)printf("Meridian/Parallel angle: %.5f\n", 239 facs.thetap * RAD_TO_DEG); 240 (void)printf("Convergence%c: ",facs.code & IS_ANAL_CONV ? '*' : ' '); 241 (void)fputs(rtodms(pline, facs.conv, 0, 0), stdout); 242 (void)printf(" [ %.8f ]\n", facs.conv * RAD_TO_DEG); 243 (void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n", 244 facs.a, facs.b); 245 } 246} 247 248int main(int argc, char **argv) { 249 char *arg, **eargv = argv, *pargv[MAX_PARGS], **iargv = argv; 250 FILE *fid; 251 int pargc = 0, iargc = argc, eargc = 0, c, mon = 0; 252 253 if (emess_dat.Prog_name = strrchr(*argv,DIR_CHAR)) 254 ++emess_dat.Prog_name; 255 else emess_dat.Prog_name = *argv; 256 inverse = ! strncmp(emess_dat.Prog_name, "inv", 3); 257 if (argc <= 1 ) { 258 (void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name); 259 exit (0); 260 } 261 /* process run line arguments */ 262 while (--argc > 0) { /* collect run line arguments */ 263 if(**++argv == '-') for(arg = *argv;;) { 264 switch(*++arg) { 265 case '\0': /* position of "stdin" */ 266 if (arg[-1] == '-') eargv[eargc++] = "-"; 267 break; 268 case 'b': /* binary I/O */ 269 bin_in = bin_out = 1; 270 continue; 271 case 'v': /* monitor dump of initialization */ 272 mon = 1; 273 continue; 274 case 'i': /* input binary */ 275 bin_in = 1; 276 continue; 277 case 'o': /* output binary */ 278 bin_out = 1; 279 continue; 280 case 'I': /* alt. method to spec inverse */ 281 inverse = 1; 282 continue; 283 case 'E': /* echo ascii input to ascii output */ 284 echoin = 1; 285 continue; 286 case 'V': /* very verbose processing mode */ 287 very_verby = 1; 288 mon = 1; 289 case 'S': /* compute scale factors */ 290 dofactors = 1; 291 continue; 292 case 't': /* set col. one char */ 293 if (arg[1]) tag = *++arg; 294 else emess(1,"missing -t col. 1 tag"); 295 continue; 296 case 'l': /* list projections, ellipses or units */ 297 if (!arg[1] || arg[1] == 'p' || arg[1] == 'P') { 298 /* list projections */ 299 struct PJ_LIST *lp; 300 int do_long = arg[1] == 'P', c; 301 char *str; 302 303 for (lp = pj_get_list_ref() ; lp->id ; ++lp) { 304 if( strcmp(lp->id,"latlong") == 0 305 || strcmp(lp->id,"longlat") == 0 306 || strcmp(lp->id,"geocent") == 0 ) 307 continue; 308 309 (void)printf("%s : ", lp->id); 310 if (do_long) /* possibly multiline description */ 311 (void)puts(*lp->descr); 312 else { /* first line, only */ 313 str = *lp->descr; 314 while ((c = *str++) && c != '\n') 315 putchar(c); 316 putchar('\n'); 317 } 318 } 319 } else if (arg[1] == '=') { /* list projection 'descr' */ 320 struct PJ_LIST *lp; 321 322 arg += 2; 323 for (lp = pj_get_list_ref(); lp->id ; ++lp) 324 if (!strcmp(lp->id, arg)) { 325 (void)printf("%9s : %s\n", lp->id, *lp->descr); 326 break; 327 } 328 } else if (arg[1] == 'e') { /* list ellipses */ 329 struct PJ_ELLPS *le; 330 331 for (le = pj_get_ellps_ref(); le->id ; ++le) 332 (void)printf("%9s %-16s %-16s %s\n", 333 le->id, le->major, le->ell, le->name); 334 } else if (arg[1] == 'u') { /* list units */ 335 struct PJ_UNITS *lu; 336 337 for (lu = pj_get_units_ref(); lu->id ; ++lu) 338 (void)printf("%12s %-20s %s\n", 339 lu->id, lu->to_meter, lu->name); 340 } else if (arg[1] == 'd') { /* list datums */ 341 struct PJ_DATUMS *ld; 342 343 printf("__datum_id__ __ellipse___ __definition/comments______________________________\n" ); 344 for (ld = pj_get_datums_ref(); ld->id ; ++ld) 345 { 346 printf("%12s %-12s %-30s\n", 347 ld->id, ld->ellipse_id, ld->defn); 348 if( ld->comments != NULL && strlen(ld->comments) > 0 ) 349 printf( "%25s %s\n", " ", ld->comments ); 350 } 351 } else 352 emess(1,"invalid list option: l%c",arg[1]); 353 exit(0); 354 continue; /* artificial */ 355 case 'e': /* error line alternative */ 356 if (--argc <= 0) 357 noargument: 358 emess(1,"missing argument for -%c",*arg); 359 oterr = *++argv; 360 continue; 361 case 'T': /* generate Chebyshev coefficients */ 362 if (--argc <= 0) goto noargument; 363 cheby_str = *++argv; 364 continue; 365 case 'm': /* cartesian multiplier */ 366 if (--argc <= 0) goto noargument; 367 postscale = 1; 368 if (!strncmp("1/",*++argv,2) || 369 !strncmp("1:",*argv,2)) { 370 if((fscale = atof((*argv)+2)) == 0.) 371 goto badscale; 372 fscale = 1. / fscale; 373 } else 374 if ((fscale = atof(*argv)) == 0.) { 375 badscale: 376 emess(1,"invalid scale argument"); 377 } 378 continue; 379 case 'W': /* specify seconds precision */ 380 case 'w': /* -W for constant field width */ 381 if ((c = arg[1]) != 0 && isdigit(c)) { 382 set_rtodms(c - '0', *arg == 'W'); 383 ++arg; 384 } else 385 emess(1,"-W argument missing or non-digit"); 386 continue; 387 case 'f': /* alternate output format degrees or xy */ 388 if (--argc <= 0) goto noargument; 389 oform = *++argv; 390 continue; 391 case 'r': /* reverse input */ 392 reversein = 1; 393 continue; 394 case 's': /* reverse output */ 395 reverseout = 1; 396 continue; 397 default: 398 emess(1, "invalid option: -%c",*arg); 399 break; 400 } 401 break; 402 } else if (**argv == '+') { /* + argument */ 403 if (pargc < MAX_PARGS) 404 pargv[pargc++] = *argv + 1; 405 else 406 emess(1,"overflowed + argument table"); 407 } else /* assumed to be input file name(s) */ 408 eargv[eargc++] = *argv; 409 } 410 if (eargc == 0 && !cheby_str) /* if no specific files force sysin */ 411 eargv[eargc++] = "-"; 412 else if (eargc > 0 && cheby_str) /* warning */ 413 emess(4, "data files when generating Chebychev prohibited"); 414 /* done with parameter and control input */ 415 if (inverse && postscale) { 416 prescale = 1; 417 postscale = 0; 418 fscale = 1./fscale; 419 } 420 if (!(Proj = pj_init(pargc, pargv))) 421 emess(3,"projection initialization failure\ncause: %s", 422 pj_strerrno(pj_errno)); 423 424 if( pj_is_latlong( Proj ) ) 425 { 426 emess( 3, "+proj=latlong unsuitable for use with proj program." ); 427 exit( 0 ); 428 } 429 430 if (inverse) { 431 if (!Proj->inv) 432 emess(3,"inverse projection not available"); 433 proj = pj_inv; 434 } else 435 proj = pj_fwd; 436 if (cheby_str) { 437 extern void gen_cheb(int, projUV(*)(projUV), char *, PJ *, int, char **); 438 439 gen_cheb(inverse, int_proj, cheby_str, Proj, iargc, iargv); 440 exit(0); 441 } 442 /* set input formating control */ 443 if (mon) { 444 pj_pr_list(Proj); 445 if (very_verby) { 446 (void)printf("#Final Earth figure: "); 447 if (Proj->es) { 448 (void)printf("ellipsoid\n# Major axis (a): "); 449 (void)printf(oform ? oform : "%.3f", Proj->a); 450 (void)printf("\n# 1/flattening: %.6f\n", 451 1./(1. - sqrt(1. - Proj->es))); 452 (void)printf("# squared eccentricity: %.12f\n", Proj->es); 453 } else { 454 (void)printf("sphere\n# Radius: "); 455 (void)printf(oform ? oform : "%.3f", Proj->a); 456 (void)putchar('\n'); 457 } 458 } 459 } 460 if (inverse) 461 informat = strtod; 462 else { 463 informat = dmstor; 464 if (!oform) 465 oform = "%.2f"; 466 } 467 468 if (bin_out) 469 { 470 SET_BINARY_MODE(stdout); 471 } 472 473 /* process input file list */ 474 for ( ; eargc-- ; ++eargv) { 475 if (**eargv == '-') { 476 fid = stdin; 477 emess_dat.File_name = "<stdin>"; 478 479 if (bin_in) 480 { 481 SET_BINARY_MODE(stdin); 482 } 483 484 } else { 485 if ((fid = fopen(*eargv, "rb")) == NULL) { 486 emess(-2, *eargv, "input file"); 487 continue; 488 } 489 emess_dat.File_name = *eargv; 490 } 491 emess_dat.File_line = 0; 492 if (very_verby) 493 vprocess(fid); 494 else 495 process(fid); 496 (void)fclose(fid); 497 emess_dat.File_name = 0; 498 } 499 if( Proj ) 500 pj_free(Proj); 501 exit(0); /* normal completion */ 502}