1 /*
2 ** lproj: projgram to exercise library of cartographic projections
3 **
4 ** Copyright (c) 2003, 2005, 2006 Gerald I. Evenden
5 */
6 static const char
7 RCS_ID[] = "Id";
8 /*
9 ** Permission is hereby granted, free of charge, to any person obtaining
10 ** a copy of this software and associated documentation files (the
11 ** "Software"), to deal in the Software without restriction, including
12 ** without limitation the rights to use, copy, modify, merge, publish,
13 ** distribute, sublicense, and/or sell copies of the Software, and to
14 ** permit persons to whom the Software is furnished to do so, subject to
15 ** the following conditions:
16 **
17 ** The above copyright notice and this permission notice shall be
18 ** included in all copies or substantial portions of the Software.
19 **
20 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21 ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
22 ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
23 ** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
24 ** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
25 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
26 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
27 */
28 #define PROJ_UV_TYPE
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <ctype.h>
32 #include <string.h>
33 #include <errno.h>
34 #include <stdarg.h>
35 #include "proj_config.h"
36 #ifdef PROJ_LIST_EXTERNAL
37 /* we must define our own proj_list since libproj4 does not */
38 # define PROJ_LIST_H "proj_list.h"
39 #endif /* PROJ_LIST_EXTERNAL */
40 #include <lib_proj.h>
41
42 #define MAX_LINE 200
43 #define MAX_PARGS 100
44 #define PROJ_INVERS(P) (P->inv ? 1 : 0)
45 static PROJ_UV
46 (*proj)(PROJ_UV, PROJ *);
47 static PROJ *Proj;
48 static int
49 reversein = 0, /* != 0 reverse input arguments */
50 reverseout = 0, /* != 0 reverse output arguments */
51 bin_in = 0, /* != 0 then binary input */
52 bin_out = 0, /* != 0 then binary output */
53 echoin = 0, /* echo input data to output line */
54 tag = '#', /* beginning of line tag character */
55 inverse = 0, /* != 0 then inverse projection */
56 prescale = 0, /* != 0 apply cartesian scale factor */
57 dofactors = 0, /* determine scale factors */
58 facs_bad = 0, /* return condition from proj_factors */
59 very_verby = 0, /* very verbose mode */
60 postscale = 0;
61 static char
62 *oform = (char *)0, /* output format for x-y or decimal degrees */
63 *oterr = "*\t*", /* output line for unprojectable input */
64 *usage =
65 "usage: %s [ -beEfiIlormsStvVwW [args] ] [ +opts[=arg] ] [ files ]\n";
66 static struct PROJ_FACTORS
67 facs;
68 static double
69 (*informat)(const char *, char **), /* input data deformatter function */
70 fscale = 0.; /* cartesian scale factor */
71 static struct {
72 char *File_name, /* input file name */
73 *Prog_name; /* name of program */
74 int File_line; /* approximate line read where error occured */
75 } emess_dat = { (char *)0, (char *)0, 0 };
76 static void
emess(int code,char * fmt,...)77 emess(int code, char *fmt, ...) {
78 va_list args;
79
80 va_start(args, fmt);
81 /* prefix program name, if given */
82 if (fmt != NULL)
83 (void)fprintf(stderr,"libproj4 system\n<%s>: ",emess_dat.Prog_name);
84 /* print file name and line, if given */
85 if (emess_dat.File_name != NULL && *emess_dat.File_name) {
86 (void)fprintf(stderr,"while processing file: %s", emess_dat.File_name);
87 if (emess_dat.File_line > 0)
88 (void)fprintf(stderr,", line %d\n", emess_dat.File_line);
89 else
90 (void)fputc('\n', stderr);
91 } else
92 putc('\n', stderr);
93 /* if |code|==2, print errno code data */
94 if (code == 2 || code == -2)
95 (void)fprintf(stderr, "Sys errno: %d: %s\n", errno, strerror(errno));
96 /* post remainder of call data */
97 (void)vfprintf(stderr,fmt,args);
98 va_end(args);
99 /* die if code positive */
100 if (code > 0) {
101 (void)fputs("\nprogram abnormally terminated\n", stderr);
102 exit(code);
103 }
104 else
105 putc('\n', stderr);
106 }
107 static void /* file processing function */
process(FILE * fid)108 process(FILE *fid) {
109 char line[MAX_LINE+3], *s = 0, pline[40];
110 PROJ_UV data;
111
112 for (;;) {
113 ++emess_dat.File_line;
114 if (bin_in) { /* binary input */
115 if (fread(&data, sizeof(PROJ_UV), 1, fid) != 1)
116 break;
117 } else { /* ascii input */
118 if (!(s = fgets(line, MAX_LINE, fid)))
119 break;
120 if (!strchr(s, '\n')) { /* overlong line */
121 int c;
122 (void)strcat(s, "\n");
123 /* gobble up to newline */
124 while ((c = fgetc(fid)) != EOF && c != '\n') ;
125 }
126 if (*s == tag) {
127 if (!bin_out)
128 (void)fputs(line, stdout);
129 continue;
130 }
131 if (reversein) {
132 data.v = (*informat)(s, &s);
133 data.u = (*informat)(s, &s);
134 } else {
135 data.u = (*informat)(s, &s);
136 data.v = (*informat)(s, &s);
137 }
138 if (data.v == HUGE_VAL)
139 data.u = HUGE_VAL;
140 if (!*s && (s > line)) --s; /* assumed we gobbled \n */
141 if (!bin_out && echoin) {
142 int t;
143 t = *s;
144 *s = '\0';
145 (void)fputs(line, stdout);
146 *s = (char)t;
147 putchar('\t');
148 }
149 }
150 if (data.u != HUGE_VAL) {
151 if (prescale) { data.u *= fscale; data.v *= fscale; }
152 if (dofactors && !inverse)
153 facs_bad = proj_factors(data, Proj, 0., &facs);
154 data = (*proj)(data, Proj);
155 if (dofactors && inverse)
156 facs_bad = proj_factors(data, Proj, 0., &facs);
157 if (postscale && data.u != HUGE_VAL)
158 { data.u *= fscale; data.v *= fscale; }
159 }
160 if (bin_out) { /* binary output */
161 (void)fwrite(&data, sizeof(PROJ_UV), 1, stdout);
162 continue;
163 } else if (data.u == HUGE_VAL) /* error output */
164 (void)fputs(oterr, stdout);
165 else if (inverse && !oform) { /*ascii DMS output */
166 if (reverseout) {
167 (void)fputs(proj_rtodms(pline, data.v, "NS"), stdout);
168 putchar('\t');
169 (void)fputs(proj_rtodms(pline, data.u, "EW"), stdout);
170 } else {
171 (void)fputs(proj_rtodms(pline, data.u, "EW"), stdout);
172 putchar('\t');
173 (void)fputs(proj_rtodms(pline, data.v, "NS"), stdout);
174 }
175 } else { /* x-y or decimal degree ascii output */
176 if (inverse) {
177 data.v *= RAD_TO_DEG;
178 data.u *= RAD_TO_DEG;
179 }
180 if (reverseout) {
181 (void)printf(oform,data.v); putchar('\t');
182 (void)printf(oform,data.u);
183 } else {
184 (void)printf(oform,data.u); putchar('\t');
185 (void)printf(oform,data.v);
186 }
187 }
188 if (dofactors) { /* print scale factor data */
189 if (!facs_bad)
190 (void)printf("\t<%.7g %.7g %g %g %g %g>",
191 facs.h, facs.k, facs.s,
192 facs.omega * RAD_TO_DEG, facs.a, facs.b);
193 else
194 (void)fputs("\t<* * * * * *>", stdout);
195 }
196 (void)fputs(bin_in ? "\n" : s, stdout);
197 }
198 }
199 static void /* file processing function --- verbosely */
vprocess(FILE * fid)200 vprocess(FILE *fid) {
201 char line[MAX_LINE+3], *s, pline[40];
202 PROJ_UV dat_ll, dat_xy;
203 int linvers;
204
205 if (!oform)
206 oform = "%.3f";
207 if (bin_in || bin_out)
208 emess(1,"binary I/O not available in -V option");
209 for (;;) {
210 ++emess_dat.File_line;
211 if (!(s = fgets(line, MAX_LINE, fid)))
212 break;
213 if (!strchr(s, '\n')) { /* overlong line */
214 int c;
215 (void)strcat(s, "\n");
216 /* gobble up to newline */
217 while ((c = fgetc(fid)) != EOF && c != '\n') ;
218 }
219 if (*s == tag) { /* pass on data */
220 (void)fputs(s, stdout);
221 continue;
222 }
223 /* check to override default input mode */
224 if (*s == 'I' || *s == 'i') {
225 linvers = 1;
226 ++s;
227 } else if (*s == 'I' || *s == 'i') {
228 linvers = 0;
229 ++s;
230 } else
231 linvers = inverse;
232 if (linvers) {
233 if (!PROJ_INVERS(Proj)) {
234 emess(-1,"inverse for this projection not avail.\n");
235 continue;
236 }
237 dat_xy.u = strtod(s, &s);
238 dat_xy.v = strtod(s, &s);
239 if (dat_xy.u == HUGE_VAL || dat_xy.v == HUGE_VAL) {
240 emess(-1,"lon-lat input conversion failure\n");
241 continue;
242 }
243 if (prescale) { dat_xy.u *= fscale; dat_xy.v *= fscale; }
244 dat_ll = proj_inv(dat_xy, Proj);
245 } else {
246 dat_ll.u = proj_dmstor(s, &s);
247 dat_ll.v = proj_dmstor(s, &s);
248 if (dat_ll.u == HUGE_VAL || dat_ll.v == HUGE_VAL) {
249 emess(-1,"lon-lat input conversion failure\n");
250 continue;
251 }
252 dat_xy = proj_fwd(dat_ll, Proj);
253 if (postscale) { dat_xy.u *= fscale; dat_xy.v *= fscale; }
254 }
255 if (proj_errno) {
256 emess(-1, proj_strerrno(proj_errno));
257 continue;
258 }
259 if (!*s && (s > line)) --s; /* assumed we gobbled \n */
260 if (proj_factors(dat_ll, Proj, 0., &facs)) {
261 emess(-1,"failed to conpute factors\n\n");
262 continue;
263 }
264 if (*s != '\n')
265 (void)fputs(s, stdout);
266 (void)fputs("Longitude: ", stdout);
267 (void)fputs(proj_rtodms(pline, dat_ll.u, "EW"), stdout);
268 (void)printf(" [ %.11g ]\n", dat_ll.u * RAD_TO_DEG);
269 (void)fputs("Latitude: ", stdout);
270 (void)fputs(proj_rtodms(pline, dat_ll.v, "NS"), stdout);
271 (void)printf(" [ %.11g ]\n", dat_ll.v * RAD_TO_DEG);
272 (void)fputs("Easting (x): ", stdout);
273 (void)printf(oform, dat_xy.u); putchar('\n');
274 (void)fputs("Northing (y): ", stdout);
275 (void)printf(oform, dat_xy.v); putchar('\n');
276 (void)printf("Meridian scale (h)%c: %.8f ( %.4g %% error )\n",
277 facs.code & IS_ANAL_HK ? '*' : ' ', facs.h, (facs.h-1.)*100.);
278 (void)printf("Parallel scale (k)%c: %.8f ( %.4g %% error )\n",
279 facs.code & IS_ANAL_HK ? '*' : ' ', facs.k, (facs.k-1.)*100.);
280 (void)printf("Areal scale (s): %.8f ( %.4g %% error )\n",
281 facs.s, (facs.s-1.)*100.);
282 (void)printf("Angular distortion (w): %.3f\n", facs.omega *
283 RAD_TO_DEG);
284 (void)printf("Meridian/Parallel angle: %.5f\n",
285 facs.thetap * RAD_TO_DEG);
286 (void)printf("Convergence%c: ",facs.code & IS_ANAL_CONV ? '*' : ' ');
287 (void)fputs(proj_rtodms(pline, facs.conv, 0), stdout);
288 (void)printf(" [ %.8f ]\n", facs.conv * RAD_TO_DEG);
289 (void)printf("Max-min (Tissot axis a-b) scale error: %.5f %.5f\n\n",
290 facs.a, facs.b);
291 }
292 }
293 int
main(int argc,char ** argv)294 main(int argc, char **argv) {
295 char *arg, **eargv = argv, *pargv[MAX_PARGS];
296 FILE *fid;
297 int pargc = 0, eargc = 0, c, mon = 0;
298
299 emess_dat.Prog_name = *argv;
300 inverse = ! strncmp(emess_dat.Prog_name, "inv", 3);
301 if (argc <= 1 ) {
302 (void)fprintf(stderr, usage, emess_dat.Prog_name);
303 exit (0);
304 }
305 /* process run line arguments */
306 while (--argc > 0) { /* collect run line arguments */
307 if(**++argv == '-') for(arg = *argv;;) {
308 switch(*++arg) {
309 case '\0': /* position of "stdin" */
310 if (arg[-1] == '-') eargv[eargc++] = "-";
311 break;
312 case 'b': /* binary I/O */
313 bin_in = bin_out = 1;
314 continue;
315 case 'v': /* monitor dump of initialization */
316 mon = 1;
317 continue;
318 case 'i': /* input binary */
319 bin_in = 1;
320 continue;
321 case 'o': /* output binary */
322 bin_out = 1;
323 continue;
324 case 'I': /* alt. method to spec inverse */
325 inverse = 1;
326 continue;
327 case 'E': /* echo ascii input to ascii output */
328 echoin = 1;
329 continue;
330 case 'V': /* very verbose processing mode */
331 very_verby = 1;
332 mon = 1;
333 case 'S': /* compute scale factors */
334 dofactors = 1;
335 continue;
336 case 't': /* set col. one char */
337 if (arg[1]) tag = *++arg;
338 else emess(1,"missing -t col. 1 tag");
339 continue;
340 case 'l': /* list projections, ellipses or units */
341 if (!arg[1] || arg[1] == 'p' || arg[1] == 'P') {
342 /* list projections */
343 const struct PROJ_LIST *lp;
344 int do_long = arg[1] == 'P';
345 char *str;
346
347 for (lp = proj_list ; lp->id ; ++lp) {
348 (void)printf("%s : ", lp->id);
349 if (do_long) /* possibly multiline description */
350 (void)puts(*lp->descr);
351 else { /* first line, only */
352 str = *lp->descr;
353 while ((c = *str++) && c != '\n')
354 putchar(c);
355 putchar('\n');
356 }
357 }
358 } else if (arg[1] == '=') { /* list projection 'descr' */
359 const struct PROJ_LIST *lp;
360
361 arg += 2;
362 for (lp = proj_list ; lp->id ; ++lp)
363 if (!strcmp(lp->id, arg)) {
364 (void)printf("%9s : %s\n", lp->id, *lp->descr);
365 break;
366 }
367 } else if (arg[1] == 'e') { /* list ellipses */
368 const struct PROJ_ELLPS *le;
369
370 for (le = proj_ellps; le->id ; ++le)
371 (void)printf("%9s %-16s %-16s %s\n",
372 le->id, le->major, le->ell, le->name);
373 } else if (arg[1] == 'u') { /* list units */
374 const struct PROJ_UNITS *lu;
375
376 for (lu = proj_units; lu->id ; ++lu)
377 (void)printf("%12s %-20s %s\n",
378 lu->id, lu->to_meter, lu->name);
379 } else
380 emess(1,"invalid list option: l%c",arg[1]);
381 exit(0);
382 continue; /* artificial */
383 case 'e': /* error line alternative */
384 if (--argc <= 0)
385 noargument: emess(1,"missing argument for -%c",*arg);
386 oterr = *++argv;
387 continue;
388 case 'm': /* cartesian multiplier */
389 if (--argc <= 0) goto noargument;
390 postscale = 1;
391 if (!strncmp("1/",*++argv,2) ||
392 !strncmp("1:",*argv,2)) {
393 if((fscale = atof((*argv)+2)) == 0.)
394 goto badscale;
395 fscale = 1. / fscale;
396 } else
397 if ((fscale = atof(*argv)) == 0.) {
398 badscale:
399 emess(1,"invalid scale argument");
400 }
401 continue;
402 case 'W': /* specify seconds precision */
403 case 'w': /* -W for constant field width */
404 if ((c = arg[1]) != 0 && isdigit(c)) {
405 proj_set_rtodms(c - '0', *arg == 'W');
406 ++arg;
407 } else
408 emess(1,"-W argument missing or non-digit");
409 continue;
410 case 'f': /* alternate output format degrees or xy */
411 if (--argc <= 0) goto noargument;
412 oform = *++argv;
413 continue;
414 case 'r': /* reverse input */
415 reversein = 1;
416 continue;
417 case 's': /* reverse output */
418 reverseout = 1;
419 continue;
420 default:
421 emess(1, "invalid option: -%c",*arg);
422 break;
423 }
424 break;
425 } else if (**argv == '+') { /* + argument */
426 if (pargc < MAX_PARGS)
427 pargv[pargc++] = *argv + 1;
428 else
429 emess(1,"overflowed + argument table");
430 } else /* assumed to be input file name(s) */
431 eargv[eargc++] = *argv;
432 }
433 if (eargc == 0) /* if no specific files force sysin */
434 eargv[eargc++] = "-";
435 /* done with parameter and control input */
436 if (inverse && postscale) {
437 prescale = 1;
438 postscale = 0;
439 fscale = 1./fscale;
440 }
441 if (!(Proj = proj_init(pargc, pargv)))
442 emess(3,"projection initialization failure\ncause: %s",
443 proj_strerrno(proj_errno));
444 if (inverse) {
445 if (!Proj->inv)
446 emess(3,"inverse projection not available");
447 proj = proj_inv;
448 } else
449 proj = proj_fwd;
450 /* set input formating control */
451 if (mon) {
452 proj_pr_list(Proj);
453 if (very_verby) {
454 (void)printf("#Final Earth figure: ");
455 if (Proj->es) {
456 (void)printf("ellipsoid\n# Major axis (a): ");
457 (void)printf(oform ? oform : "%.3f", Proj->a);
458 (void)printf("\n# 1/flattening: %.6f\n",
459 1./(1. - sqrt(1. - Proj->es)));
460 (void)printf("# squared eccentricity: %.12f\n", Proj->es);
461 } else {
462 (void)printf("sphere\n# Radius: ");
463 (void)printf(oform ? oform : "%.3f", Proj->a);
464 (void)putchar('\n');
465 }
466 }
467 }
468 if (inverse)
469 informat = strtod;
470 else {
471 informat = proj_dmstor;
472 if (!oform)
473 oform = "%.2f";
474 }
475 /* process input file list */
476 for ( ; eargc-- ; ++eargv) {
477 if (**eargv == '-') {
478 fid = stdin;
479 emess_dat.File_name = "<stdin>";
480 } else {
481 if ((fid = fopen(*eargv, "r")) == NULL) {
482 emess(-2, *eargv, "input file");
483 continue;
484 }
485 emess_dat.File_name = *eargv;
486 }
487 emess_dat.File_line = 0;
488 if (very_verby)
489 vprocess(fid);
490 else
491 process(fid);
492 (void)fclose(fid);
493 emess_dat.File_name = 0;
494 }
495 return 0; /* normal completion */
496 }
497 /*
498 ** Log: lproj.c
499 ** Revision 1.2 2008-11-07 21:40:43 jeff
500 ** ENH: Fixing some proj.4 warnings.
501 **
502 ** Revision 1.1 2008-11-07 16:41:13 jeff
503 ** ENH: Adding a 2D geoview. Adding the geographic projection library libproj4
504 ** to Utilities. Updating the architecture of the geospatial views. All
505 ** multi-resolution sources are now subclasses of vtkGeoSource. Each source
506 ** has its own worker thread for fetching refined images or geometry.
507 ** On the 3D side, vtkGeoGlobeSource is an appropriate source for vtkGeoTerrain,
508 ** and vtkGeoAlignedImageSource is an appropriate source for
509 ** vtkGeoAlignedImageRepresentation. On the 2D side, vtkGeoProjectionSource is an
510 ** appropriate source for vtkGeoTerrain2D, and the image source is the same.
511 **
512 ** Revision 3.2 2006/01/24 02:00:12 gie
513 ** additional corrections
514 **
515 ** Revision 3.1 2006/01/11 02:41:14 gie
516 ** Initial
517 **
518 */
519