svn-gvsig-desktop / tags / v2_0_Build_1217 / libraries / libjni-proj4 / src / proj.c @ 44178
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 |
} |