7 #define NVERT 20 /* max number of vertices in a -v polygon */
8 #define HALFWIDTH 8192 /* output scaled to fit in -HALFWIDTH,HALFWIDTH */
9 #define LONGLINES (HALFWIDTH*4) /* permissible segment lengths */
10 #define SHORTLINES (HALFWIDTH/8)
11 #define SCALERATIO 10 /* of abs to rel data (see map(5)) */
12 #define RESOL 2. /* coarsest resolution for tracing grid (degrees) */
13 #define TWO_THRD 0.66666666666666667
15 int normproj(double, double, double *, double *);
16 int posproj(double, double, double *, double *);
17 int picut(struct place *, struct place *, double *);
18 double reduce(double);
19 short getshort(FILE *);
20 char *mapindex(char *);
24 static char *mapdir = "#9/map"; /* default map directory */
30 static struct file dfltfile = {
31 "world", BLACK, SOLID /* default map */
33 static struct file *file = &dfltfile; /* list of map files */
34 static int nfile = 1; /* length of list */
35 static char *currcolor = BLACK; /* current color */
36 static char *gridcolor = BLACK;
37 static char *bordcolor = BLACK;
39 extern struct index index[];
40 int halfwidth = HALFWIDTH;
42 static int (*cut)(struct place *, struct place *, double *);
43 static int (*limb)(double*, double*, double);
44 static void dolimb(void);
47 static double orientation[3] = { 90., 0., 0. }; /* -o option */
48 static int oriented; /* nonzero if -o option occurred */
49 static int upright; /* 1 if orientation[0]==90, -1 if -90, else 0*/
50 static int delta = 1; /* -d setting */
51 static double limits[4] = { /* -l parameters */
52 -90., 90., -180., 180.
54 static double klimits[4] = { /* -k parameters */
55 -90., 90., -180., 180.
58 static double rlimits[4]; /* limits expressed in radians */
59 static double lolat, hilat, lolon, hilon;
60 static double window[4] = { /* option -w */
61 -90., 90., -180., 180.
63 static int windowed; /* nozero if option -w */
64 static struct vert { double x, y; } v[NVERT+2]; /*clipping polygon*/
65 static struct edge { double a, b, c; } e[NVERT]; /* coeffs for linear inequality */
66 static int nvert; /* number of vertices in clipping polygon */
68 static double rwindow[4]; /* window, expressed in radians */
69 static double params[2]; /* projection params */
70 /* bounds on output values before scaling; found by coarse survey */
71 static double xmin = 100.;
72 static double xmax = -100.;
73 static double ymin = 100.;
74 static double ymax = -100.;
75 static double xcent, ycent;
76 static double xoff, yoff;
77 double xrange, yrange;
78 static int left = -HALFWIDTH;
79 static int right = HALFWIDTH;
80 static int bottom = -HALFWIDTH;
81 static int top = HALFWIDTH;
82 static int longlines = SHORTLINES; /* drop longer segments */
83 static int shortlines = SHORTLINES;
84 static int bflag = 1; /* 0 for option -b */
85 static int s1flag = 0; /* 1 for option -s1 */
86 static int s2flag = 0; /* 1 for option -s2 */
87 static int rflag = 0; /* 1 for option -r */
88 static int kflag = 0; /* 1 if option -k occurred */
89 static int xflag = 0; /* 1 for option -x */
90 int vflag = 1; /* -1 if option -v occurred */
91 static double position[3]; /* option -p */
92 static double center[3] = {0., 0., 0.}; /* option -c */
93 static struct coord crot; /* option -c */
94 static double grid[3] = { 10., 10., RESOL }; /* option -g */
95 static double dlat, dlon; /* resolution for tracing grid in lat and lon */
96 static double scaling; /* to compute final integer output */
97 static struct file *track; /* options -t and -u */
98 static int ntrack; /* number of tracks present */
99 static char *symbolfile; /* option -y */
101 void clamp(double *px, double v);
103 double diddle(struct place *, double, double);
104 double diddle(struct place *, double, double);
105 void dobounds(double, double, double, double, int);
106 void dogrid(double, double, double, double);
107 int duple(struct place *, double);
108 #define fmax map_fmax
109 #define fmin map_fmin
110 double fmax(double, double);
111 double fmin(double, double);
112 void getdata(char *);
113 int gridpt(double, double, int);
114 int inpoly(double, double);
115 int inwindow(struct place *);
116 void pathnames(void);
118 void radbds(double *w, double *rw);
119 void revlon(struct place *, double);
120 void satellite(struct file *);
121 int seeable(double, double);
129 if(s[0]=='-' && (s[1]<'0'||s[1]>'9'))
130 return(s[1]!='.'&&s[1]!=0);
136 conv(int k, struct coord *g)
138 g->l = (0.0001/SCALERATIO)*k;
143 main(int argc, char *argv[])
152 abort(); /* getshort() won't work */
153 mapdir = unsharp(mapdir);
157 s = getenv("MAPDIR");
161 error("usage: map projection params options");
162 for(k=0;index[k].name;k++) {
166 if(*s==0) goto found;
171 fprintf(stderr,"projections:\n");
172 for(i=0;index[i].name;i++) {
173 fprintf(stderr,"%s",index[i].name);
174 for(k=0; k<index[i].npar; k++)
175 fprintf(stderr," p%d", k);
176 fprintf(stderr,"\n");
183 limb = index[k].limb;
184 poles = index[k].poles;
185 for(i=0;i<index[k].npar;i++) {
186 if(i>=argc||option(argv[i])) {
187 fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar);
190 params[i] = atof(argv[i]);
194 while(argc>0&&option(argv[0])) {
197 switch(argv[-1][1]) {
199 if(file == &dfltfile) {
203 while(argc && !option(*argv)) {
204 file = realloc(file,(nfile+1)*sizeof(*file));
205 file[nfile].name = *argv;
206 file[nfile].color = currcolor;
207 file[nfile].style = SOLID;
215 for(nvert=0;nvert<NVERT&&argc>=2;nvert++) {
218 v[nvert].x = atof(*argv++);
222 v[nvert].y = atof(*argv++);
226 error("too many clipping vertices");
229 gridcolor = currcolor;
230 for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
231 grid[i] = atof(argv[i]);
234 grid[0] = grid[1] = 0.;
248 while(argc && !option(*argv)) {
249 track = realloc(track,(ntrack+1)*sizeof(*track));
250 track[ntrack].name = *argv;
251 track[ntrack].color = currcolor;
252 track[ntrack].style = style;
262 switch(argv[-1][2]) {
266 case 0: /* compatibility */
272 for(i=0;i<3&&i<argc&&!option(argv[i]);i++)
273 orientation[i] = atof(argv[i]);
279 bordcolor = currcolor;
280 for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
281 limits[i] = atof(argv[i]);
287 for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
288 klimits[i] = atof(argv[i]);
293 if(argc>0&&!option(argv[0])) {
294 delta = atoi(argv[0]);
300 bordcolor = currcolor;
302 for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
303 window[i] = atof(argv[i]);
308 for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
309 center[i] = atof(argv[i]);
314 for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
315 position[i] = atof(argv[i]);
318 if(i!=3||position[2]<=0)
319 error("incomplete positioning");
322 if(argc>0&&!option(argv[0])) {
323 symbolfile = argv[0];
329 if(index[k].limb == 0)
330 error("-v does not apply here");
337 if(argc && !option(*argv)) {
338 currcolor = colorcode(*argv);
346 error("error in arguments");
348 clamp(&limits[0],-90.);
349 clamp(&limits[1],90.);
350 clamp(&klimits[0],-90.);
351 clamp(&klimits[1],90.);
352 clamp(&window[0],-90.);
353 clamp(&window[1],90.);
354 radbds(limits,rlimits);
355 limcase = limits[2]<-180.?0:
359 window[0]>=window[1]||
360 window[2]>=window[3]||
365 error("unreasonable window");
367 radbds(window,rwindow);
368 upright = orientation[0]==90? 1: orientation[0]==-90? -1: 0;
369 if(index[k].spheroid && !upright)
370 error("can't tilt the spheroid");
371 if(limits[2]>limits[3])
374 orientation[2] = (limits[2]+limits[3])/2;
375 orient(orientation[0],orientation[1],orientation[2]);
376 projection = (*index[k].prog)(params[0],params[1]);
378 error("unreasonable projection parameters");
380 grid[0] = fabs(grid[0]);
381 grid[1] = fabs(grid[1]);
384 klimits[i] = limits[i];
385 if(klimits[2]>klimits[3])
391 if(lolon>=hilon||lolat>=hilat||lolat<-90.||hilat>90.)
392 error("unreasonable limits");
393 wlim = kflag? klimits: window;
394 dlat = fmin(hilat-lolat,wlim[1]-wlim[0])/16;
395 dlon = fmin(hilon-lolon,wlim[3]-wlim[2])/32;
396 dd = fmax(dlat,dlon);
397 while(grid[2]>fmin(dlat,dlon)/2)
401 for(lat=klimits[0];lat<klimits[1]+dd-FUZZ;lat+=dd) {
404 for(lon=klimits[2];lon<klimits[3]+dd-FUZZ;lon+=dd) {
405 i = (kflag?posproj:normproj)
406 (lat,lon+(lon<klimits[3]?FUZZ:-FUZZ),
417 for(i=0; i<nvert; i++) {
426 xrange = xmax - xmin;
427 yrange = ymax - ymin;
428 if(xrange<=0||yrange<=0)
429 error("map seems to be empty");
430 scaling = 2; /*plotting area from -1 to 1*/
432 if(posproj(position[0]-.5,position[1],&xcent,&ycent)==0||
433 posproj(position[0]+.5,position[1],&x,&y)==0)
434 error("unreasonable position");
435 scaling /= (position[2]*hypot(x-xcent,y-ycent));
436 if(posproj(position[0],position[1],&xcent,&ycent)==0)
437 error("unreasonable position");
439 scaling /= (xrange>yrange?xrange:yrange);
440 xcent = (xmin+xmax)/2;
441 ycent = (ymin+ymax)/2;
443 xoff = center[0]/scaling;
444 yoff = center[1]/scaling;
445 crot.l = center[2]*RAD;
447 scaling *= HALFWIDTH*0.9;
454 range(left,bottom,right,top);
459 for(lat=ceil(lolat/grid[0])*grid[0];
460 lat<=hilat;lat+=grid[0])
461 dogrid(lat,lat,lolon,hilon);
463 for(lon=ceil(lolon/grid[1])*grid[1];
464 lon<=hilon;lon+=grid[1])
465 dogrid(lolat,hilat,lon,lon);
466 comment("border","");
471 dobounds(lolat,hilat,lolon,hilon,0);
472 dobounds(window[0],window[1],window[2],window[3],1);
474 lolat = floor(limits[0]/10)*10;
475 hilat = ceil(limits[1]/10)*10;
476 lolon = floor(limits[2]/10)*10;
477 hilon = ceil(limits[3]/10)*10;
480 /*do tracks first so as not to lose the standard input*/
481 for(i=0;i<ntrack;i++) {
482 longlines = LONGLINES;
483 satellite(&track[i]);
484 longlines = shortlines;
486 for(i=0;i<nfile;i++) {
487 comment("mapfile",file[i].name);
488 colorx(file[i].color);
490 getdata(file[i].name);
498 /* Out of perverseness (really to recover from a dubious,
499 but documented, convention) the returns from projection
500 functions (-1 unplottable, 0 wrong sheet, 1 good) are
501 recoded into -1 wrong sheet, 0 unplottable, 1 good. */
504 fixproj(struct place *g, double *x, double *y)
506 int i = (*projection)(g,x,y);
507 return i<0? 0: i==0? -1: 1;
511 normproj(double lat, double lon, double *x, double *y)
515 latlon(lat,lon,&geog);
522 i = fixproj(&geog,x,y);
527 fprintf(stderr,"%d %.3f %.3f\n",i,*x,*y);
533 posproj(double lat, double lon, double *x, double *y)
537 latlon(lat,lon,&geog);
539 i = fixproj(&geog,x,y);
546 inwindow(struct place *geog)
548 if(geog->nlat.l<rwindow[0]-FUZZ||
549 geog->nlat.l>rwindow[1]+FUZZ||
550 geog->wlon.l<rwindow[2]-FUZZ||
551 geog->wlon.l>rwindow[3]+FUZZ)
557 inlimits(struct place *g)
559 if(rlimits[0]-FUZZ>g->nlat.l||
560 rlimits[1]+FUZZ<g->nlat.l)
564 if(rlimits[2]+TWOPI-FUZZ>g->wlon.l&&
565 rlimits[3]+FUZZ<g->wlon.l)
569 if(rlimits[2]-FUZZ>g->wlon.l||
570 rlimits[3]+FUZZ<g->wlon.l)
574 if(rlimits[2]>g->wlon.l&&
575 rlimits[3]-TWOPI+FUZZ<g->wlon.l)
586 getdata(char *mapfile)
601 indexfile = mapindex(mapfile);
602 xfile = fopen(indexfile,"r");
604 filerror("can't find map index", indexfile);
606 for(i=0,p=patch[0];i<18*36;i++,p++)
608 while(!feof(xfile) && fscanf(xfile,"%d%d%ld",&i,&j,&b)==3)
609 patch[i+9][j+18] = b;
611 ifile = fopen(mapfile,"r");
613 filerror("can't find map data", mapfile);
614 for(lat=lolat;lat<hilat;lat+=10.)
615 for(lon=lolon;lon<hilon;lon+=10.) {
616 if(!seeable(lat,lon))
620 if((b=patch[i+9][j+18])&1)
623 while((ip=getc(ifile))>=0&&(jp=getc(ifile))>=0){
624 if(ip!=(i&0377)||jp!=(j&0377))
628 if(n > 0) { /* absolute coordinates */
629 kx = ky = 0; /* set */
631 kx = SCALERATIO*getshort(ifile);
632 ky = SCALERATIO*getshort(ifile);
633 if (((k%delta) != 0) && (k != (n-1)))
637 conn = plotpt(&g,conn);
639 } else { /* differential, scaled by SCALERATI0 */
641 kx = SCALERATIO*getshort(ifile);
642 ky = SCALERATIO*getshort(ifile);
645 if(c&0200) c|= ~0177;
648 if(c&0200) c|= ~0177;
650 if(k%delta!=0&&k!=n-1)
654 conn = plotpt(&g,conn);
668 seeable(double lat0, double lon0)
672 for(lat=lat0;lat<=lat0+10;lat+=2*grid[2])
673 for(lon=lon0;lon<=lon0+10;lon+=2*grid[2])
674 if(normproj(lat,lon,&x,&y)*vflag>0)
680 satellite(struct file *t)
692 if(t->name[0]!='-'||t->name[1]!=0) {
694 if((ifile=fopen(t->name,"r"))==NULL)
695 filerror("can't find track", t->name);
697 comment("track",t->name);
702 while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){
703 latlon(lat,lon,&place);
704 if(fscanf(ifile,"%1s",lbl) == 1) {
705 if(strchr("+-.0123456789",*lbl)==0)
709 conn = plotpt(&place,conn);
713 fscanf(ifile,"%[^\n]",lbl+1);
716 if(plotpt(&place,conn))
721 if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1)
723 if(plotpt(&place,conn?conn:-1)) {
724 int r = *lbl=='!'?0:rflag?-1:1;
726 if(putsym(&place,sym,scale,r) == 0)
732 if(plotpt(&place,conn))
745 if(i>=18) return(i-36);
746 if(i<-18) return(i+36);
753 fprintf(stderr,"map: \r\n%s\n",s);
758 filerror(char *s, char *f)
760 fprintf(stderr,"\r\n%s %s\n",s,f);
767 char *t = malloc(strlen(s)+3);
774 static int ox = NOPT;
775 static int oy = NOPT;
778 cpoint(int xi, int yi, int conn)
782 if(!xflag && (xi<left||xi>=right || yi<bottom||yi>=top)) {
786 if(conn == -1) /* isolated plotting symbol */
791 if(dx+dy>longlines) {
792 ox = oy = NOPT; /* don't leap across cuts */
799 return dx+dy<=2? 2: 1; /* 2=very near; see dogrid */
806 plotpt(struct place *g, int conn)
818 switch((*cut)(g,&oldg,&cutlon)) {
821 ret = duple(g,cutlon)|duple(g,cutlon);
827 default: /* prevent diags about bad return value */
830 ret = doproj(g,&kx,&ky);
831 if(ret==0 || !onlimb && ret*vflag<=0)
833 ret = cpoint(kx,ky,conn);
839 doproj(struct place *g, int *kx, int *ky)
843 /*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
844 i = fixproj(g,&x,&y);
849 /*fprintf(stderr,"dopr2 %f %f\n",x,y);*/
855 x = (x1*crot.c - y1*crot.s + xoff)*scaling;
856 y = (x1*crot.s + y1*crot.c + yoff)*scaling;
857 *kx = x + (x>0?.5:-.5);
858 *ky = y + (y>0?.5:-.5);
863 duple(struct place *g, double cutlon)
869 revlon(&oldg,cutlon);
874 if(doproj(g,&kx,&ky)*vflag<=0 ||
875 doproj(&oldg,&okx,&oky)*vflag<=0)
883 revlon(struct place *g, double cutlon)
885 g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon));
890 /* recognize problems of cuts
891 * move a point across cut to side of its predecessor
892 * if its very close to the cut
893 * return(0) if cut interrupts the line
894 * return(1) if line is to be drawn normally
895 * return(2) if line is so close to cut as to
896 * be properly drawn on both sheets
900 picut(struct place *g, struct place *og, double *cutlon)
903 return(ckcut(g,og,PI));
907 nocut(struct place *g, struct place *og, double *cutlon)
921 ckcut(struct place *g1, struct place *g2, double lon)
926 d1 = reduce(g1->wlon.l -lon);
927 d2 = reduce(g2->wlon.l -lon);
928 if((f1=fabs(d1))<FUZZ)
929 d1 = diddle(g1,lon,d2);
930 if((f2=fabs(d2))<FUZZ) {
931 d2 = diddle(g2,lon,d1);
932 if(doproj(g2,&kx,&ky)*vflag>0)
937 if(f1>PI*TWO_THRD||f2>PI*TWO_THRD)
943 diddle(struct place *g, double lon, double d)
949 g->wlon.l = reduce(lon+d1);
965 double tetrapt = 35.26438968; /* atan(1/sqrt(2)) */
968 dogrid(double lat0, double lat1, double lon0, double lon1)
970 double slat,slon,tlat,tlon;
971 register int conn, oconn;
972 slat = tlat = slon = tlon = 0;
974 slat = tlat = fmin(grid[2],dlat);
976 slon = tlon = fmin(grid[2],dlon);;
978 while(lat0<=lat1&&lon0<=lon1) {
979 conn = gridpt(lat0,lon0,conn);
980 if(projection==Xguyou&&slat>0) {
981 if(lat0<-45&&lat0+slat>-45)
982 conn = gridpt(-45.,lon0,conn);
983 else if(lat0<45&&lat0+slat>45)
984 conn = gridpt(45.,lon0,conn);
985 } else if(projection==Xtetra&&slat>0) {
986 if(lat0<-tetrapt&&lat0+slat>-tetrapt) {
987 gridpt(-tetrapt-.001,lon0,conn);
988 conn = gridpt(-tetrapt+.001,lon0,0);
990 else if(lat0<tetrapt&&lat0+slat>tetrapt) {
991 gridpt(tetrapt-.001,lon0,conn);
992 conn = gridpt(tetrapt+.001,lon0,0);
995 if(conn==0 && oconn!=0) {
997 lat0 -= slat; /* steps too big */
998 lon0 -= slon; /* or near bdry */
1001 conn = oconn = gridpt(lat0,lon0,conn);
1015 gridpt(lat1,lon1,conn);
1018 static int gridinv; /* nonzero when doing window bounds */
1021 gridpt(double lat, double lon, int conn)
1024 /*fprintf(stderr,"%f %f\n",lat,lon);*/
1028 return(plotpt(&g,conn));
1031 /* win=0 ordinary grid lines, win=1 window lines */
1034 dobounds(double lolat, double hilat, double lolon, double hilon, int win)
1037 if(lolat>-90 || win && (poles&1)!=0)
1038 dogrid(lolat+FUZZ,lolat+FUZZ,lolon,hilon);
1039 if(hilat<90 || win && (poles&2)!=0)
1040 dogrid(hilat-FUZZ,hilat-FUZZ,lolon,hilon);
1041 if(hilon-lolon<360 || win && cut==picut) {
1042 dogrid(lolat,hilat,lolon+FUZZ,lolon+FUZZ);
1043 dogrid(lolat,hilat,hilon-FUZZ,hilon-FUZZ);
1052 double res = fmin(dlat, dlon)/4;
1057 onlimb = gridinv = 1;
1059 newconn = (*limb)(&lat, &lon, res);
1062 conn = gridpt(lat, lon, conn*newconn);
1064 onlimb = gridinv = 0;
1069 radbds(double *w, double *rw)
1083 double center = orientation[0];
1086 center = 180 - center;
1088 center = -180 - center;
1090 error("unreasonable orientation");
1091 colat = 90 - window[0];
1092 if(center-colat>limits[0])
1093 limits[0] = center - colat;
1094 if(center+colat<limits[1])
1095 limits[1] = center + colat;
1104 r = (c | getc(f)<<8);
1106 r |= ~0xFFFF; /* in case short > 16 bits */
1111 fmin(double x, double y)
1117 fmax(double x, double y)
1123 clamp(double *px, double v)
1125 *px = (v<0?fmax:fmin)(*px,v);
1132 char *t, *indexfile, *name;
1134 for(i=0; i<nfile; i++) {
1135 name = file[i].name;
1138 indexfile = mapindex(name);
1139 /* ansi equiv of unix access() call */
1140 f = fopen(name, "r");
1141 fx = fopen(indexfile, "r");
1147 t = malloc(strlen(name)+strlen(mapdir)+2);
1162 for(i=0; i<nvert; i++) { /*convert latlon to xy*/
1163 if(normproj(v[i].x,v[i].y,&v[i].x,&v[i].y)==0)
1164 error("invisible clipping vertex");
1166 if(nvert==2) { /*rectangle with diag specified*/
1169 v[1].x=v[0].x, v[1].y=v[2].y, v[3].x=v[2].x, v[3].y=v[0].y;
1174 for(i=1; i<=nvert; i++) { /*test for convexity*/
1175 t = (v[i-1].x-v[i].x)*(v[i+1].y-v[i].y) -
1176 (v[i-1].y-v[i].y)*(v[i+1].x-v[i].x);
1177 if(t<-FUZZ && s>=0) s = 1;
1178 if(t>FUZZ && s<=0) s = -1;
1179 if(-FUZZ<=t&&t<=FUZZ || t*s>0) {
1185 error("improper clipping polygon");
1186 for(i=0; i<nvert; i++) { /*edge equation ax+by=c*/
1187 e[i].a = s*(v[i+1].y - v[i].y);
1188 e[i].b = s*(v[i].x - v[i+1].x);
1189 e[i].c = s*(v[i].x*v[i+1].y - v[i].y*v[i+1].x);
1194 inpoly(double x, double y)
1197 for(i=0; i<nvert; i++) {
1198 register struct edge *ei = &e[i];
1199 double val = x*ei->a + y*ei->b - ei->c;
1212 if(cut != picut) /* punt on unusual cuts */
1214 for(lat=window[0]; lat<=window[1]; lat+=grid[2]) {
1227 longlines = shortlines = LONGLINES;
1228 cut = nocut; /* not necessary; small eff. gain */