Statistics
| Revision:

root / tags / v2_0_Build_1210 / libraries / libjni-proj4 / src / proj.c @ 38161

History | View | Annotate | Download (16.5 KB)

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 <stdio.h>
6
#include <stdlib.h>
7
#include <ctype.h>
8
#include <string.h>
9
#include <math.h>
10
#include "projects.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 200
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
31
reversein = 0,        /* != 0 reverse input arguments */
32
reverseout = 0,        /* != 0 reverse output arguments */
33
bin_in = 0,        /* != 0 then binary input */
34
bin_out = 0,        /* != 0 then binary output */
35
echoin = 0,        /* echo input data to output line */
36
tag = '#',        /* beginning of line tag character */
37
inverse = 0,        /* != 0 then inverse projection */
38
prescale = 0,        /* != 0 apply cartesian scale factor */
39
dofactors = 0,        /* determine scale factors */
40
facs_bad = 0,        /* return condition from pj_factors */
41
very_verby = 0, /* very verbose mode */
42
postscale = 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
50
facs;
51
        static double
52
(*informat)(const char *, char **),        /* input data deformatter function */
53
fscale = 0.;        /* cartesian scale factor */
54
        static projUV
55
int_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 */
63
process(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 */
154
vprocess(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

    
248
int 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
    exit(0); /* normal completion */
500
}