1 /* Copyright (C) 1988-2005 by Brian Doty and the Institute
2 of Global Environment and Society (IGES).
3
4 See file COPYRIGHT for more information. */
5
6 /* Authored by B. Doty */
7
8 #ifdef HAVE_CONFIG_H
9 #include <config.h>
10
11 /* If autoconfed, only include malloc.h when it's presen */
12 #ifdef HAVE_MALLOC_H
13 #include <malloc.h>
14 #endif
15
16 #else /* undef HAVE_CONFIG_H */
17
18 #include <malloc.h>
19
20 #endif /* HAVE_CONFIG_H */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <math.h>
25 #include <float.h>
26 #include "gx.h"
27 #include "gxmap.h"
28
29
30 static int sizes[9];
31 static float *mapx, *mapy;
32 static int ifirst=0;
33 static FILE *mapin, *imap;
34
35 static int adjtyp = 0; /* Direction adjustment class */
36 static float lonref; /* Reference longitude for adjustment */
37
gxdmap(struct mapopt * mopt)38 void gxdmap (struct mapopt *mopt) {
39 float lon[255],lat[255],xx,yy,lnmin,lnmax,ltmin,ltmax,lnfact;
40 int num,i,ipen,rc,type,ilon,ilat,rnum,flag;
41 int st1,st2,spos;
42 float sln1,sln2,slt1,slt2;
43 float lnsav,ltsav,lndif,ltdif,lntmp,lttmp,llinc,llsum,lldist;
44 char *fname;
45 char hdr[3],rec[1530];
46
47 llinc = hypot(mopt->lnmax-mopt->lnmin,mopt->ltmax-mopt->ltmin);
48 llinc = llinc/200;
49 if (llinc<0.0001) llinc=0.0001;
50
51 /* Open the map data set */
52
53 if (*(mopt->mpdset)=='/' || *(mopt->mpdset)=='\\') {
54 imap = fopen(mopt->mpdset,"rb");
55 if (imap==NULL) {
56 printf ("Open Error on Map Data Set: %s\n",mopt->mpdset);
57 return;
58 }
59 } else {
60 fname = gxgnam(mopt->mpdset);
61 imap = fopen(fname,"rb");
62 if (imap==NULL) {
63 imap = fopen(mopt->mpdset,"rb");
64 if (imap==NULL) {
65 printf ("Open Error on Map Data Set: %s\n",fname);
66 free (fname);
67 return;
68 }
69 }
70 free (fname);
71 }
72
73 /* Read and process each record */
74
75 rnum = 0;
76 while (1) {
77 rc = fread(hdr,1,3,imap);
78 if (rc!=3) break;
79 rnum++;
80 i = gagby (hdr,0,1);
81 if (i<1 || i>3) {
82 printf ("Map file format error: Invalid rec type %i rec num %i\n",i,rnum);
83 return;
84 }
85 if (i==2) {
86 st1 = gagby(hdr,1,1);
87 st2 = gagby(hdr,2,1);
88 fread(rec,1,16,imap);
89 spos = gagby(rec,0,4);
90 ilon = gagby(rec,4,3);
91 sln1 = ((float)ilon)/1e4;
92 ilon = gagby(rec,7,3);
93 sln2 = ((float)ilon)/1e4;
94 ilat = gagby(rec,10,3);
95 slt1 = ((float)ilat)/1e4 - 90.0;
96 ilat = gagby(rec,13,3);
97 slt2 = ((float)ilat)/1e4 - 90.0;
98 flag = 0;
99 for (i=0; i<256; i++) {
100 if (*(mopt->mcol+i)!=-9 && i>=st1 && i<=st2) flag = 1;
101 }
102 if (flag==0) {
103 if (spos==0) {
104 fclose(imap);
105 return;
106 }
107 fseek(imap,spos,0);
108 continue;
109 }
110 flag = 0;
111 if (sln1>360.0) flag = 1;
112 else {
113 if (slt2 <= mopt->ltmin || slt1 >= mopt->ltmax) flag = 0;
114 else {
115 lnfact = 0.0;
116 while (sln2+lnfact > mopt->lnmin) lnfact -= 360.0;
117 lnfact += 360.0;
118 if (sln1+lnfact >= mopt->lnmax) flag = 0;
119 else flag = 1;
120 }
121 }
122 if (flag==0) {
123 if (spos==0) {
124 fclose(imap);
125 return;
126 }
127 fseek(imap,spos,0);
128 }
129 continue;
130 }
131 type = gagby(hdr,1,1);
132 num = gagby(hdr,2,1);
133
134 /* Read the next record; convert the data points;
135 and get the lat/lon bounds for this line segment */
136
137 fread(rec,1,num*6,imap);
138 if (*(mopt->mcol+type) == -9) continue;
139 if (*(mopt->mcol+type) == -1) {
140 gxcolr(mopt->dcol);
141 gxstyl(mopt->dstl);
142 gxwide(mopt->dthk);
143 } else {
144 gxcolr(*(mopt->mcol+type));
145 gxstyl(*(mopt->mstl+type));
146 gxwide(*(mopt->mthk+type));
147 }
148 lnmin = 9999.9; lnmax = -9999.9; ltmin = 9999.9; ltmax = -9999.9;
149 for (i=0; i<num; i++) {
150 ilon = gagby(rec,i*6,3);
151 ilat = gagby(rec,i*6+3,3);
152 lat[i] = ((float)ilat)/1e4 - 90.0;
153 lon[i] = ((float)ilon)/1e4;
154 if (lat[i]<ltmin) ltmin=lat[i]; if (lat[i]>ltmax) ltmax=lat[i];
155 if (lon[i]<lnmin) lnmin=lon[i]; if (lon[i]>lnmax) lnmax=lon[i];
156 }
157
158 /* Plot this line segment if it falls within the
159 appropriate lat/lon bounds */
160
161 if (ltmax < mopt->ltmin) continue;
162 if (ltmin > mopt->ltmax) continue;
163
164 lnfact = 0.0;
165 while (lnmax+lnfact > mopt->lnmin) lnfact -= 360.0;
166 lnfact += 360.0;
167
168 while (lnmin+lnfact < mopt->lnmax) {
169 if (lnmax+lnfact < mopt->lnmin) {
170 lnfact += 360.0;
171 continue;
172 }
173
174 /* Split long lines into shorter segments and limit
175 drawing at lat-lon bounds */
176
177 ipen = 3;
178 lnsav = lon[0]; ltsav = lat[0];
179 for (i=1; i<num; i++) {
180 lndif = fabs(lon[i] - lon[i-1]);
181 ltdif = fabs(lat[i] - lat[i-1]);
182 if (lndif>ltdif) lldist = lndif;
183 else lldist = ltdif;
184 llsum = llinc;
185 lntmp = lnsav; lttmp = ltsav;
186 while (llsum<lldist+llinc) {
187 if (llsum>=lldist-llinc/4.0) {
188 lntmp = lon[i]; lttmp = lat[i];
189 llsum += llinc; /* Insure loop dropout */
190 } else {
191 if (lndif>ltdif) {
192 if (lon[i-1]<lon[i]) {
193 lntmp += llinc;
194 lttmp += llinc * (lat[i]-lat[i-1])/(lon[i]-lon[i-1]);
195 } else {
196 lntmp -= llinc;
197 lttmp -= llinc * (lat[i]-lat[i-1])/(lon[i]-lon[i-1]);
198 }
199 } else {
200 if (lat[i-1]<lat[i]) {
201 lttmp += llinc;
202 lntmp += llinc * (lon[i]-lon[i-1])/(lat[i]-lat[i-1]);
203 } else {
204 lttmp -= llinc;
205 lntmp -= llinc * (lon[i]-lon[i-1])/(lat[i]-lat[i-1]);
206 }
207 }
208 }
209 if (lntmp+lnfact<mopt->lnmin ||
210 lntmp+lnfact>mopt->lnmax ||
211 lttmp<mopt->ltmin || lttmp>mopt->ltmax) {
212 if (ipen==2) {
213 gxconv (lntmp+lnfact,lttmp,&xx,&yy,2);
214 gxplot (xx,yy,ipen);
215 }
216 ipen = 3;
217 } else {
218 if (ipen==3) {
219 gxconv (lnsav+lnfact,ltsav,&xx,&yy,2);
220 gxplot (xx,yy,ipen);
221 }
222 ipen = 2;
223 gxconv (lntmp+lnfact,lttmp,&xx,&yy,2);
224 gxplot (xx,yy,ipen);
225 }
226 lnsav = lntmp; ltsav = lttmp;
227 llsum += llinc;
228 }
229 }
230 lnfact += 360.0;
231 }
232 }
233 fclose (imap);
234 }
235
236 /* Routine to set up scaling for lat-lon projection. The aspect
237 ratio is *not* maintained. */
238
gxscld(struct mapprj * mpj,int xflip,int yflip)239 int gxscld (struct mapprj *mpj, int xflip, int yflip) {
240 float x1,x2,y1,y2;
241
242 if (mpj->lnmn>=mpj->lnmx) return(1);
243 if (mpj->ltmn>=mpj->ltmx) return(1);
244 if (mpj->xmn>=mpj->xmx) return(1);
245 if (mpj->ymn>=mpj->ymx) return(1);
246 mpj->axmn = mpj->xmn;
247 mpj->axmx = mpj->xmx;
248 mpj->aymn = mpj->ymn;
249 mpj->aymx = mpj->ymx;
250 x1 = mpj->lnmn; x2 = mpj->lnmx; y1 = mpj->ltmn; y2 = mpj->ltmx;
251 if (xflip) { x1 = mpj->lnmx; x2 = mpj->lnmn; }
252 if (yflip) { y1 = mpj->ltmx; y2 = mpj->ltmn; }
253 gxscal (mpj->axmn, mpj->axmx, mpj->aymn, mpj->aymx, x1, x2, y1, y2);
254 gxproj (NULL);
255 adjtyp = 0;
256 return (0);
257 }
258
259 /* Routine to set up scaling for lat-lon projection. Aspect
260 ratio of the projection is maintained as a constant, and it
261 fills the plotting area as much as possible. */
262
gxltln(struct mapprj * mpj)263 int gxltln (struct mapprj *mpj) {
264 float lndif,ltdif,aspect,aspect2,xdif,xlo,xhi,ydif,ylo,yhi;
265
266 if (mpj->lnmn>=mpj->lnmx) return(1);
267 if (mpj->ltmn>=mpj->ltmx) return(1);
268 if (mpj->xmn>=mpj->xmx) return(1);
269 if (mpj->ymn>=mpj->ymx) return(1);
270
271 lndif = mpj->lnmx - mpj->lnmn;
272 ltdif = mpj->ltmx - mpj->ltmn;
273 aspect = 1.2*ltdif/lndif;
274 aspect2 = (mpj->ymx - mpj->ymn) / (mpj->xmx - mpj->xmn);
275 if (aspect > aspect2) {
276 xdif = (mpj->xmx - mpj->xmn) * aspect2/aspect;
277 xlo = ((mpj->xmx - mpj->xmn)/2.0)-(xdif*0.5);
278 xhi = ((mpj->xmx - mpj->xmn)/2.0)+(xdif*0.5);
279 mpj->axmx = mpj->xmn + xhi;
280 mpj->axmn = mpj->xmn + xlo;
281 mpj->aymn = mpj->ymn;
282 mpj->aymx = mpj->ymx;
283 } else {
284 ydif = (mpj->ymx - mpj->ymn) * aspect/aspect2;
285 ylo = ((mpj->ymx - mpj->ymn)/2.0)-(ydif*0.5);
286 yhi = ((mpj->ymx - mpj->ymn)/2.0)+(ydif*0.5);
287 mpj->aymx = mpj->ymn + yhi;
288 mpj->aymn = mpj->ymn + ylo;
289 mpj->axmn = mpj->xmn;
290 mpj->axmx = mpj->xmx;
291 }
292
293 gxscal (mpj->axmn, mpj->axmx, mpj->aymn, mpj->aymx,
294 mpj->lnmn, mpj->lnmx, mpj->ltmn, mpj->ltmx);
295 gxproj (NULL);
296 adjtyp = 0;
297 return (0);
298 }
299
300 /* Routine for north polar stereographic. Projection scaling
301 is set along with level 1 linear scaling. The only difficult
302 aspect to this is to set the level 1 linear scaling such that
303 the proper aspect ratio is maintained. */
304
305 static float londif;
306
gxnste(struct mapprj * mpj)307 int gxnste (struct mapprj *mpj) {
308 float x1,x2,y1,y2,dum,lonave;
309 float w1,xave,yave;
310 float lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;
311
312 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
313 latmn = mpj->ltmn; latmx = mpj->ltmx;
314 xmin = mpj->xmn; xmax = mpj->xmx;
315 ymin = mpj->ymn; ymax = mpj->ymx;
316
317 if ((lonmx-lonmn) > 360.0) {
318 return (1);
319 }
320 if (lonmn<-360.0 || lonmx>360.0) {
321 return (1);
322 }
323 if (latmn<-80.0 || latmx>90.0) {
324 return (1);
325 }
326 if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) {
327 return(1);
328 }
329
330 lonave = (lonmx+lonmn)/2.0; /* Longitude adjustment to put */
331 londif = -90.0 - lonave; /* central meridian at bottom.*/
332 lonref = lonave;
333
334 /* Plotting limits depend on how much of the hemisphere we are
335 actually plotting. */
336
337 if ( (lonmx-lonmn) < 180.0 ) {
338
339 gxnpst ( lonmn, latmn, &x1, &dum ); /* Left side coord */
340 gxnpst ( lonmx, latmn, &x2, &dum ); /* Right side coord */
341 gxnpst ( lonmn, latmx, &dum, &y2 ); /* Top coord */
342 gxnpst ( lonave, latmn, &dum, &y1 ); /* Bottom coord */
343
344 } else {
345
346 gxnpst ( lonave-90.0, latmn, &x1, &dum ); /* Left side coord */
347 gxnpst ( lonave+90.0, latmn, &x2, &dum ); /* Right side coord */
348 gxnpst ( lonmn, latmn, &dum, &y2 ); /* Top coord */
349 gxnpst ( lonave, latmn, &dum, &y1 ); /* Bottom coord */
350
351 }
352
353 /* Set up linear level scaling while maintaining aspect ratio. */
354
355 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
356 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
357 xave = (xmax+xmin)/2.0;
358 gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
359 mpj->axmn = xave-w1; mpj->axmx = xave+w1;
360 mpj->aymn = ymin; mpj->aymx = ymax;
361 } else {
362 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
363 yave = (ymax+ymin)/2.0;
364 gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
365 mpj->axmn = xmin; mpj->axmx = xmax;
366 mpj->aymn = yave-w1; mpj->aymx = yave+w1;
367 }
368
369 gxproj (gxnpst);
370 gxback (gxnrev);
371 adjtyp = 1;
372 return (0);
373
374 }
375
gxnpst(float rlon,float rlat,float * x,float * y)376 void gxnpst (float rlon, float rlat, float *x, float *y) {
377 float radius,theta;
378
379 radius = tan (0.785315-(0.00872572*rlat));
380 theta = (rlon+londif)*0.0174514;
381 *x = radius * cos(theta);
382 *y = radius * sin(theta);
383 }
384
385 /* Routine for back transform for npst */
386
gxnrev(float x,float y,float * rlon,float * rlat)387 void gxnrev (float x, float y, float *rlon, float *rlat) {
388 float rad,alpha,theta;
389
390 rad = hypot(x,y);
391 alpha = 180.0*atan(rad)/3.14159;
392 *rlat = 90.0 - 2.0*alpha;
393
394 if (x==0.0 && y==0.0) *rlon = 0.0;
395 else {
396 *rlon = (180.0*atan2(y,x)/3.14159)-londif;
397 while (*rlon < lonref-180.0) *rlon += 360.0;
398 while (*rlon > lonref+180.0) *rlon -= 360.0;
399 }
400 }
401
402 /* Routine for south polar stereographic. Projection scaling
403 is set along with level 1 linear scaling. The only difficult
404 aspect to this is to set the level 1 linear scaling such that
405 the proper aspect ratio is maintained. */
406
gxsste(struct mapprj * mpj)407 int gxsste (struct mapprj *mpj) {
408 float x1,x2,y1,y2,dum,lonave;
409 float w1,xave,yave;
410 float lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;
411
412 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
413 latmn = mpj->ltmn; latmx = mpj->ltmx;
414 xmin = mpj->xmn; xmax = mpj->xmx;
415 ymin = mpj->ymn; ymax = mpj->ymx;
416
417 if ((lonmx-lonmn) > 360.0) {
418 return (1);
419 }
420 if (lonmn<-360.0 || lonmx>360.0) {
421 return (1);
422 }
423 if (latmn<-90.0 || latmx>80.0) {
424 return (1);
425 }
426 if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) {
427 return(1);
428 }
429
430 lonave = (lonmx+lonmn)/2.0; /* Longitude adjustment to put */
431 londif = -90.0 - lonave; /* central meridian at bottom.*/
432 lonref = lonave;
433
434 /* Plotting limits depend on how much of the hemisphere we are
435 actually plotting. */
436
437 if ( (lonmx-lonmn) < 180.0 ) {
438
439 gxspst ( lonmn, latmx, &x1, &dum ); /* Left side coord */
440 gxspst ( lonmx, latmx, &x2, &dum ); /* Right side coord */
441 gxspst ( lonmn, latmn, &dum, &y1 ); /* Top coord */
442 gxspst ( lonave, latmx, &dum, &y2 ); /* Bottom coord */
443
444 } else {
445
446 gxspst ( lonave-90.0, latmx, &x1, &dum ); /* Left side coord */
447 gxspst ( lonave+90.0, latmx, &x2, &dum ); /* Right side coord */
448 gxspst ( lonmn, latmx, &dum, &y1 ); /* Top coord */
449 gxspst ( lonave, latmx, &dum, &y2 ); /* Bottom coord */
450
451 }
452
453 /* Set up linear level scaling while maintaining aspect ratio. */
454
455 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
456 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
457 xave = (xmax+xmin)/2.0;
458 gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
459 mpj->axmn = xave-w1; mpj->axmx = xave+w1;
460 mpj->aymn = ymin; mpj->aymx = ymax;
461 } else {
462 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
463 yave = (ymax+ymin)/2.0;
464 gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
465 mpj->axmn = xmin; mpj->axmx = xmax;
466 mpj->aymn = yave-w1; mpj->aymx = yave+w1;
467 }
468 gxproj (gxspst);
469 gxback (gxsrev);
470 adjtyp = 2;
471 return (0);
472
473 }
474
gxspst(float rlon,float rlat,float * x,float * y)475 void gxspst (float rlon, float rlat, float *x, float *y) {
476 float radius,theta;
477
478 radius = tan(0.785315+(0.00872572*rlat));
479 theta = (rlon+londif)*(-0.0174514);
480 *x = radius * cos(theta);
481 *y = radius * sin(theta);
482 }
483
484 /* Routine for back transform for spst */
485
gxsrev(float x,float y,float * rlon,float * rlat)486 void gxsrev (float x, float y, float *rlon, float *rlat) {
487 float rad,alpha,theta;
488
489 rad = hypot(x,y);
490 alpha = 180.0*atan(rad)/3.14159;
491 *rlat = 2.0*alpha - 90.0;
492
493 if (x==0.0 && y==0.0) *rlon = 0.0;
494 else {
495 *rlon = (-180.0*atan2(y,x)/3.14159)-londif;
496 while (*rlon < lonref-180.0) *rlon += 360.0;
497 while (*rlon > lonref+180.0) *rlon -= 360.0;
498 }
499 }
500
501 /* Return adjustment angle (in radians) to apply to a wind direction
502 to correct for current map projection and position. */
503
gxaarw(float lon,float lat)504 float gxaarw (float lon, float lat) {
505
506 if (adjtyp==0) return(0.0);
507 if (adjtyp==1) {
508 lon = (lon - lonref)*3.14159/180.0;
509 return (lon);
510 }
511 if (adjtyp==2) {
512 lon = (lonref - lon)*3.14159/180.0;
513 return (lon);
514 }
515 return (0.0);
516 }
517
518 /* Set up Robinson Projection */
519
gxrobi(struct mapprj * mpj)520 int gxrobi (struct mapprj *mpj) {
521 float lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;
522 float x1,x2,y1,y2,xd,yd,xave,yave,w1;
523
524 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
525 latmn = mpj->ltmn; latmx = mpj->ltmx;
526 xmin = mpj->xmn; xmax = mpj->xmx;
527 ymin = mpj->ymn; ymax = mpj->ymx;
528
529 /* Check for errors */
530
531 if (lonmn<-180.0 || lonmx>180.0 || latmn<-90.0 || latmx>90.0) {
532 return (1);
533 }
534 if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) {
535 return(1);
536 }
537
538 /* Get bounds of the map in linear units */
539
540 gxrobp ( lonmn, latmn, &x1, &y1); /* Lower Left */
541 gxrobp ( lonmn, latmx, &xd, &y2); /* Upper Left */
542 if (xd<x1) x1 = xd;
543 if (latmn<0.0 && latmx>0.0) {
544 gxrobp (lonmn, 0.0, &xd, &yd); /* Left Middle */
545 if (xd<x1) x1 = xd;
546 }
547 gxrobp ( lonmx, latmn, &x2, &y1); /* Lower Right */
548 gxrobp ( lonmx, latmx, &xd, &y2); /* Upper Right */
549 if (xd>x2) x2 = xd;
550 if (latmn<0.0 && latmx>0.0) {
551 gxrobp (lonmx, 0.0, &xd, &yd); /* Right Middle */
552 if (xd>x2) x2 = xd;
553 }
554
555 /* Set up linear level scaling while maintaining aspect ratio. */
556
557 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
558 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
559 xave = (xmax+xmin)/2.0;
560 gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
561 mpj->axmn = xave-w1; mpj->axmx = xave+w1;
562 mpj->aymn = ymin; mpj->aymx = ymax;
563 } else {
564 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
565 yave = (ymax+ymin)/2.0;
566 gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
567 mpj->axmn = xmin; mpj->axmx = xmax;
568 mpj->aymn = yave-w1; mpj->aymx = yave+w1;
569 }
570
571 gxproj (gxrobp);
572 gxback (gxrobb);
573 adjtyp = 3;
574 return (0);
575 }
576
577 /* Transform routine for Robinson Projection */
578
579 float rob1[37] = {-1.349,-1.317,-1.267,-1.206,-1.138,-1.066,-0.991,
580 -0.913,-0.833,-0.752,-0.669,-0.586,-0.502,-0.418,-0.334,-0.251,
581 -0.167,-0.084,0.000,0.084,0.167,0.251,0.334,0.418,0.502,0.586,
582 0.669,0.752,0.833,0.913,0.991,1.066,1.138,1.206,1.267,1.317,1.349};
583 float rob2[37] = {1.399,1.504,1.633,1.769,1.889,1.997,2.099,
584 2.195,2.281,2.356,2.422,2.478,2.532,2.557,2.582,2.602,2.616,
585 2.625,2.628,2.625,2.616,2.602,2.582,2.557,2.532,2.478,2.422,
586 2.356,2.281,2.195,2.099,1.997,1.889,1.769,1.633,1.504,1.399};
587
gxrobp(float rlon,float rlat,float * x,float * y)588 void gxrobp (float rlon, float rlat, float *x, float *y) {
589 int i;
590 rlat = (rlat+90.0)/5.0;
591 i = (int)rlat;
592 rlat = rlat - (float)i;
593 if (i<0) {
594 *y = -1.349;
595 *x = 1.399*rlon/180.0;
596 return;
597 }
598 if (i>=36) {
599 *y = 1.349;
600 *x = 1.399*rlon/180.0;
601 return;
602 }
603 *y = rob1[i] + rlat*(rob1[i+1]-rob1[i]);
604 *x = rob2[i] + rlat*(rob2[i+1]-rob2[i]);
605 *x = *x * rlon/180.0;
606 return;
607 }
608
609 /* Back Transform for Robinson Projection */
610
gxrobb(float x,float y,float * rlon,float * rlat)611 void gxrobb (float x, float y, float *rlon, float *rlat) {
612 }
613 /*------------------------------------------------------------------
614 DKRZ appends: Mollweide Projection
615 10.08.95 Karin Meier (karin.meier@dkrz.de)
616 ------------------------------------------------------------------*/
617
gxmoll(struct mapprj * mpj)618 int gxmoll (struct mapprj *mpj) {
619 float lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;
620 float x1,x2,y1,y2,xd,yd,xave,yave,w1;
621
622 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
623 latmn = mpj->ltmn; latmx = mpj->ltmx;
624 xmin = mpj->xmn; xmax = mpj->xmx;
625 ymin = mpj->ymn; ymax = mpj->ymx;
626 lomin = lonmn; lomax = lonmx;
627 lamin = latmn; lamax = latmx;
628
629 /* Check for errors */
630
631 if (latmn<-90.0 || latmx>90.0) {
632 return (1);
633 }
634 if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) {
635 return(1);
636 }
637
638 /* Get bounds of the map in linear units */
639
640 gxmollp ( lonmn, latmn, &x1, &y1); /* Lower Left */
641 gxmollp ( lonmn, latmx, &xd, &y2); /* Upper Left */
642 if (xd<x1) x1 = xd;
643 if (latmn<0.0 && latmx>0.0) {
644 gxmollp (lonmn, 0.0, &xd, &yd); /* Left Middle */
645 if (xd<x1) x1 = xd;
646 }
647 gxmollp ( lonmx, latmn, &x2, &y1); /* Lower Right */
648 gxmollp ( lonmx, latmx, &xd, &y2); /* Upper Right */
649 if (xd>x2) x2 = xd;
650 if (latmn<0.0 && latmx>0.0) {
651 gxmollp (lonmx, 0.0, &xd, &yd); /* Right Middle */
652 if (xd>x2) x2 = xd;
653 }
654
655 /* Set up linear level scaling while maintaining aspect ratio. */
656
657 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
658 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
659 xave = (xmax+xmin)/2.0;
660 gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
661 mpj->axmn = xave-w1; mpj->axmx = xave+w1;
662 mpj->aymn = ymin; mpj->aymx = ymax;
663 } else {
664 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
665 yave = (ymax+ymin)/2.0;
666 gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
667 mpj->axmn = xmin; mpj->axmx = xmax;
668 mpj->aymn = yave-w1; mpj->aymx = yave+w1;
669 }
670
671 gxproj (gxmollp);
672 gxback (gxmollb);
673 adjtyp = 3;
674 return (0);
675 }
676
gxmollp(float rlon,float rlat,float * x,float * y)677 void gxmollp (float rlon, float rlat, float *x, float *y) {
678 float diff, radlat, radlon;
679
680 if (lomin != -180.0) {
681 diff = -180.0 - lomin;
682 rlon = rlon + diff;
683 }
684 radlat = (pi*rlat)/180.0;
685 radlon = (pi*rlon)/180.0;
686
687 *x = cos(radlat);
688 *y = sin(radlat)/2.0;
689 *x = *x*rlon/180.0;
690
691 return;
692 }
693
694 /* Back Transform for Mollweide Projection */
695
gxmollb(float x,float y,float * rlon,float * rlat)696 void gxmollb (float x, float y, float *rlon, float *rlat) {
697 }
698
699
700 /*------------------------------------------------------------------
701 DKRZ appends: Orthographic Projection
702 15.08.95 Karin Meier (karin.meier@dkrz.d400.de)
703 ------------------------------------------------------------------*/
704
gxortg(struct mapprj * mpj)705 int gxortg (struct mapprj *mpj) {
706 float lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;
707 float lonave;
708 float x1,x2,y1,y2,xd,yd,xave,yave,w1;
709
710 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
711 latmn = mpj->ltmn; latmx = mpj->ltmx;
712 xmin = mpj->xmn; xmax = mpj->xmx;
713 ymin = mpj->ymn; ymax = mpj->ymx;
714 lomin = lonmn; lomax = lonmx;
715 lamin = latmn; lamax = latmx;
716 lonref = (lonmx+lonmn)/2.0;
717
718 /* Check boundaries */
719
720 if (latmn != -90.0 || latmx != 90.0) {
721 printf("Map Projection Error: Latitude must be in range -90 90\n");
722 return (1);
723 }
724 if ((lonmx - lonmn) > 180.0 ) {
725 printf("Map Projection Error: %.1f - %.1f > 180.0\n",
726 lonmx, lonmn);
727 return (1);
728 }
729 if ((lonmx - lonmn) < 180.0) {
730 printf("Map Projection Error: %.1f - %.1f < 180.0\n",
731 lonmx, lonmn);
732 return (1);
733 }
734 if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) return(1);
735
736 if (lonmn < -180.0) {
737 mpj->lnmn = lonmn + 360.0;
738 mpj->lnmx = lonmx + 360.0;
739 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
740 }
741 if (lonmx > 180.0 ) {
742 mpj->lnmn = lonmn - 360.0;
743 mpj->lnmx = lonmx - 360.0;
744 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
745 }
746
747 /* Get bounds of the map in linear units */
748
749 gxortgp ( lonmn, latmn, &x1, &y1);
750 gxortgp ( lonmn, latmx, &xd, &y2);
751 if (xd<x1) x1 = xd;
752 if (latmn<0.0 && latmx>0.0) {
753 gxortgp ( lonmn, 0.0, &xd, &yd);
754 if (xd<x1) x1 = xd;
755 }
756 gxortgp ( lonmx, latmn, &x2, &y1);
757 gxortgp ( lonmx, latmx, &xd, &y2);
758 if(xd>x2) x2 = xd;
759 if (latmn<0.0 && latmx>0.0) {
760 gxortgp ( lonmx, 0.0, &xd, &yd);
761 if (xd>x2) x2 = xd;
762 }
763
764 /* Set up linear level scaling while maintaining aspect ratio. */
765
766 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
767 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
768 xave = (xmax+xmin)/2.0;
769 gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
770 mpj->axmn = xave-w1; mpj->axmx = xave+w1;
771 mpj->aymn = ymin; mpj->aymx = ymax;
772 } else {
773 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
774 yave = (ymax+ymin)/2.0;
775 gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
776 mpj->axmn = xmin; mpj->axmx = xmax;
777 mpj->aymn = yave-w1; mpj->aymx = yave+w1;
778 }
779
780 gxproj (gxortgp);
781 gxback (gxortgb);
782 adjtyp = 3;
783 return (0);
784 }
785
gxortgp(float rlon,float rlat,float * x,float * y)786 void gxortgp (float rlon, float rlat, float *x, float *y) {
787 int i;
788 float radlat, radlon, diff;
789
790 if ( lomin != -90.0 ){
791 diff = -90.0 - lomin;
792 rlon = rlon + diff;
793 }
794
795 radlat = (pi*rlat)/180.0;
796 radlon = (pi*rlon)/180.0;
797
798 *x = cos(radlat);
799 *y = sin(radlat);
800 *x = *x * sin(radlon);
801
802 return;
803 }
804
805 /* Back Transform for Orthographic Projection */
806
gxortgb(float x,float y,float * rlon,float * rlat)807 void gxortgb (float x, float y, float *rlon, float *rlat) {
808 }
809
810 /*------------------------------------------------------------------
811 DKRZ appends: Lambert conformal conic Projection
812 15.03.96 Karin Meier (karin.meier@dkrz.de)
813 ------------------------------------------------------------------*/
814 static float hemi, r;
815
gxlamc(struct mapprj * mpj)816 int gxlamc (struct mapprj *mpj) {
817 float lonmn, lonmx, latmn, latmx, dlat, dlon, dx, dy;
818 float xave,yave, w1, lonave, xmin, xmax, ymin, ymax, x1, x2, y1, y2, xd, yd;
819 int flag=0;
820
821 lonmn = mpj->lnmn; lonmx = mpj->lnmx;
822 latmn = mpj->ltmn; latmx = mpj->ltmx;
823 xmin = mpj->xmn; xmax = mpj->xmx;
824 ymin = mpj->ymn; ymax = mpj->ymx;
825 lomin = lonmn; lomax = lonmx;
826 lamin = latmn; lamax = latmx;
827 lonave = (lonmx+lonmn)/2.0;
828 dlat = lamax - lamin; dlon = lomax - lomin;
829 dx = xmax - xmin; dy = ymax - ymin;
830
831 if ((lonmn>=lonmx)||(latmn>=latmx)||(xmin>=xmax)||(ymin>=ymax)) {
832 return(1);
833 }
834 if (((latmn > 0.0) && (latmx < 0.0)) || ((latmn < 0.0) && (latmx >0.0))) {
835 printf("Map Projection Error: Latitude must be in range -90 0 or 0 90\n");
836 return (1);
837 }
838
839 /*--- set constant for northern or southern hemisphere ---*/
840
841 if (latmn >= 0.0) {
842 hemi = 1.0; /** northern hemisphere **/
843 } else {
844 hemi = -1.0; /** southern hemisphere **/
845 }
846
847 /*--- reset 90.0/-90.0 degrees to 89.99/-89.99 because of tangent ---*/
848
849 if (latmn == -90.0) latmn = -89.99;
850 if (latmx == 90.0) latmx = 89.99;
851
852 /*--- get viewport coordinates x1, x2, y1, y2---*/
853
854 gxlamcp (lonmn, latmn, &x1, &y1);
855 gxlamcp (lonmn, latmx, &xd, &y2);
856 if (xd<x1) x1=xd;
857 if (y2<y1) {
858 yd = y2; y2 = y1; y1 = yd;
859 }
860 if (latmn>=0.0 && latmx>0.0) {
861 gxlamcp (lonmn,0.0,&xd,&yd);
862 if (xd<x1) x1=xd;
863 }
864 gxlamcp (lonmx, latmn, &x2, &y1);
865 gxlamcp (lonmx, latmx, &xd, &y2);
866 if (xd>x2) x2=xd;
867 if (y2<y1) {
868 yd = y2; y2 = y1; y1 = yd;
869 }
870 if (latmn<0.0 && latmx<=0.0) {
871 gxlamcp (lonmx,0.0,&xd,&yd);
872 if (xd>x2) x2=xd;
873 }
874
875 /*--- determining terms for scaling ---*/
876
877 xave = (xmin+xmax)/2.0;
878 yave = (ymin+ymax)/2.0;
879
880 if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) )
881 {
882 if (hemi==-1.0 && 180.0<(lomax-lomin) && (lomax-lomin)<=270.0)
883 yave -= 1.5;
884 else if (hemi==1.0 && 180.0<(lomax-lomin) && (lomax-lomin)<=270.0)
885 yave += 1.5;
886 else if (hemi==-1.0 && 270.0<=(lomax-lomin) && (lomax-lomin)<=360.0)
887 yave -= 1.2;
888 else if (hemi==1.0 && 270.0<=(lomax-lomin) && (lomax-lomin)<=360.0)
889 yave += 1.2;
890 else if (hemi==-1.0 && 90.0<(lomax-lomin) && (lomax-lomin)<=180.0)
891 yave -= 0.5;
892 else if (hemi==1.0 && 90.0<(lomax-lomin) && (lomax-lomin)<=180.0)
893 yave += 1.0;
894 else if (hemi==-1.0 && (lomax-lomin)<=90.0)
895 yave += 0.0;
896 else if (hemi==1.0 && (lomax-lomin)<=90.0)
897 yave += 1.0;
898
899 w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
900 if (w1 < 1.0)
901 gxscal (xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
902 else if (w1 < 2.0)
903 gxscal (xave-0.5*(w1), xave+0.5*w1, yave-w1, yave+w1,
904 x1, x2, y1, y2 );
905 else if (w1 < 3.0)
906 gxscal (xave-0.5*w1, xave+0.5*w1, yave-w1, yave+w1,
907 x1, x2, y1, y2 );
908 else if (w1 > 3.0)
909 gxscal (xave-0.75*w1, xave+0.75*w1, yave-0.75*w1,
910 yave+0.75*w1, x1, x2, y1, y2 );
911 else
912 gxscal (xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
913 }
914 else
915 {
916 if (hemi==-1.0 && 180.0<(lomax-lomin) && (lomax-lomin)<=270.0)
917 yave -= 1.0;
918 else if (hemi==1.0 && 180.0<(lomax-lomin) && (lomax-lomin)<=270.0)
919 yave += 1.5;
920 else if (hemi==-1.0 && 270.0<=(lomax-lomin) && (lomax-lomin)<=360.0)
921 yave -= 1.0;
922 else if (hemi==1.0 && 270.0<=(lomax-lomin) && (lomax-lomin)<=360.0)
923 yave += 1.0;
924 else if (hemi==-1.0 && 90.0<(lomax-lomin) && (lomax-lomin)<=180.0)
925 yave -= 0.5;
926 else if (hemi==1.0 && 90.0<(lomax-lomin) && (lomax-lomin)<=180.0)
927 yave += 1.0;
928 else if (hemi==-1.0 && (lomax-lomin)<=90.0)
929 yave += 0.0;
930 else if (hemi==1.0 && (lomax-lomin)<=90.0)
931 yave += 1.0;
932
933 w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
934 if (w1 < 1.0)
935 gxscal (xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
936 else if (w1 < 2.0)
937 gxscal (xmin+0.5*w1, xmax-0.5*w1, yave-1.25*w1,
938 yave+1.25*w1, x1, x2, y1, y2 );
939 else if (w1 < 3.0)
940 gxscal (xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
941 else if (w1 > 3.0)
942 gxscal (xave-0.5*w1, xave+0.5*w1, yave-0.5*w1,
943 yave+0.5*w1, x1, x2, y1, y2 );
944 else
945 gxscal (xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
946 }
947
948 mpj->axmn = xmin; mpj->axmx = xmax;
949 mpj->aymn = ymin; mpj->aymx = ymax;
950
951 gxproj (gxlamcp);
952 gxback (gxlamcb);
953 adjtyp = 0;
954 return (0);
955 }
956
957
958
959 /*--- transform routine for lambert conformal conic projection ---*/
960
gxlamcp(float rlon,float rlat,float * x,float * y)961 void gxlamcp (float rlon, float rlat, float *x, float *y)
962 {
963 float d2r, cone, phis, phin, clon, term1, term2;
964
965 d2r = pi/180.0;
966
967 /*--- standard latitudes: north - phin; south - phis ---*/
968 phis = lamin;
969 phin = lamax;
970
971 /*--- reset 90.0/-90.0 degrees to 89.99/-89.99 because of tangent ---*/
972 if(phis == -90.0) phis = -89.99;
973 if(phin == 90.0) phin = 89.99;
974
975 /*--- calculate the constant of the cone +++ radius, x, y ---*/
976 /*--- clon - central meridian; cone - cone constant ---*/
977 clon = floor((lomax + lomin)/2.0);
978 term1 = tan((45.0-hemi*phis/2.0)*d2r);
979 term2 = tan((45.0-hemi*phin/2.0)*d2r);
980
981 if(phis!=phin)
982 cone = (log10(cos(phis*d2r))-log10(cos(phin*d2r)))/
983 (log10(term1)-log10(term2));
984 else
985 cone = cos((90.0-hemi*phis)*d2r);
986
987 r = pow(tan((45.0-hemi*rlat/2.0)*d2r),cone);
988 *x = r*sin((rlon-clon)*d2r*cone);
989 *y = -hemi*r*cos((rlon-clon)*d2r*cone);
990
991 return;
992 }
993
994
995
996 /*--- Back Transform for Lambert conformal Projection ---*/
997
gxlamcb(float x,float y,float * rlon,float * rlat)998 void gxlamcb (float x, float y, float *rlon, float *rlat) {
999 }
1000
1001 /* Interpolate lat/lon boundaries, and convert to xy, on
1002 behalf of 'draw mappoly' . For most part, the same
1003 code as in gxdmap */
1004
gxmpoly(float * xy,int cnt,float llinc,int * newcnt)1005 float *gxmpoly(float *xy, int cnt, float llinc, int *newcnt) {
1006 float ln1, ln2, lt1, lt2, lnsav, ltsav, llsum;
1007 float lndif, ltdif, lldist, lntmp, lttmp, xx, yy;
1008 float *newxy;
1009 int i,j,ip,ncnt;
1010
1011 /* Determine total 'path' length */
1012
1013 llsum = 0.0;
1014 for (i=1; i<cnt; i++) {
1015 ip = (i-1)*2;
1016 lndif = fabs(*(xy+ip+2) - *(xy+ip));
1017 ltdif = fabs(*(xy+ip+3) - *(xy+ip+1));
1018 if (lndif>ltdif) lldist = lndif;
1019 else lldist = ltdif;
1020 llsum += lldist;
1021 }
1022
1023 /* Estimate number of output points, and allocate storage
1024 for them. */
1025
1026 ncnt = cnt + llsum/llinc;
1027 newxy = (float *)malloc(sizeof(float)*ncnt*2);
1028 if (newxy==NULL) return(NULL); /* caller issues error */
1029
1030 /* Now interpolate each point, convert to x,y, and put in list */
1031
1032 j = 0;
1033 lnsav = *xy; ltsav = *(xy+1);
1034 for (i=1; i<cnt; i++) {
1035 ip = (i-1)*2;
1036 ln1 = *(xy+ip); ln2 = *(xy+ip+2);
1037 lt1 = *(xy+ip+1); lt2 = *(xy+ip+3);
1038 lndif = fabs(ln2-ln1);
1039 ltdif = fabs(lt2-lt1);
1040 if (lndif>ltdif) lldist = lndif;
1041 else lldist = ltdif;
1042 llsum = llinc;
1043 lntmp = lnsav; lttmp = ltsav;
1044 while (llsum<lldist+llinc) {
1045 if (llsum>=lldist-llinc/4.0) {
1046 lntmp = ln2; lttmp = lt2;
1047 llsum += llinc; /* Insure loop dropout */
1048 } else {
1049 if (lndif>ltdif) {
1050 if (ln1<ln2) {
1051 lntmp += llinc;
1052 lttmp += llinc * (lt2-lt1)/(ln2-ln1);
1053 } else {
1054 lntmp -= llinc;
1055 lttmp -= llinc * (lt2-lt1)/(ln2-ln1);
1056 }
1057 } else {
1058 if (lt1<lt2) {
1059 lttmp += llinc;
1060 lntmp += llinc * (ln2-ln1)/(lt2-lt1);
1061 } else {
1062 lttmp -= llinc;
1063 lntmp -= llinc * (ln2-ln1)/(lt2-lt1);
1064 }
1065 }
1066 }
1067 gxconv (lntmp,lttmp,&xx,&yy,2);
1068 *(newxy+j*2) = xx;
1069 *(newxy+j*2+1) = yy;
1070 j++;
1071 if (j>=ncnt) {
1072 printf ("Logic Error in gxmpoly\n");
1073 free (newxy);
1074 return (NULL);
1075 }
1076 lnsav = lntmp; ltsav = lttmp;
1077 llsum += llinc;
1078 }
1079 }
1080 *newcnt = j;
1081 return (newxy);
1082 }
1083