PageRenderTime 99ms CodeModel.GetById 60ms app.highlight 34ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/proj.c

http://github.com/route-me/route-me
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}