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