/Proj4/proj.c

http://github.com/route-me/route-me · C · 502 lines · 472 code · 19 blank · 11 comment · 193 complexity · 9af58c2e6db7ca82ee3006e321e3b20c MD5 · raw file

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