ric@rioja.ifs.umich.edu (Richard Campbell) (07/17/90)
I've got a copy of Brian Reid's netmaps. I'd like to have an example of the algorithm used to translate between a given latitude and longitude to the position on the map. I've tried some empirical forumals, but nothing is quite working out. thanks, ric
reid@wrl.dec.com (Brian Reid) (08/15/90)
In article <1990Jul16.193115.26400@terminator.cc.umich.edu> ric@rioja.ifs.umich.edu (Richard Campbell) writes: >I've got a copy of Brian Reid's netmaps. I'd like >to have an example of the algorithm used to translate >between a given latitude and longitude to the >position on the map. I've tried some empirical forumals, >but nothing is quite working out. /* These are the functions that implement various cartographic projections. At the moment there are implementations of Mercator projection and Lambert Conformal projection, though it is a simple matter of adding 4 functions to implement any new kind of projection. Brian K. Reid, DEC Western Research, 1986 Copyright (c) 1986, 1990 Digital Equipment Corporation All Rights Reserved */ #include "netmap.h" /* ------------------------------------------------------------------------ */ /* This is the user-level function that projects a latitude and longitude onto the cartesian plane, storing into x and y. */ project(lat,longit,x,y) double lat,longit,*x,*y; { double PolarR, PolarTheta, T; switch(Projection) { case Mercator: mercat(lat,longit,x,y); break; case Lambert: lambert(lat,longit,x,y); break; case Gall: gallpr(lat,longit,x,y); break; } *x -= CenterX; *y -= CenterY; switch (RotDegrees) { case 0: break; case 90: T = *y; *y = - *x; *x = T; break; case -90: T = *y; *y = *x; *x = -T; break; default: PolarR = sqrt (*x * *x + *y * *y); PolarTheta = rotate + atan2(*x,*y); *x = PolarR*sin(PolarTheta); *y = PolarR*cos(PolarTheta); } #ifdef DEBUG { double ddlat,ddlong; fprintf(stderr,"P(%3.1f,%3.1f) --> [%3.1f,%3.1f]\n",lat,longit,*x,*y); unproject(*x,*y,&ddlat,&ddlong); } #endif DEBUG } /* ------------------------------------------------------------------------ */ /* This is a reverse projection: given a cartesian point, return the latitude and longitude that project to that point. */ unproject(xin,yin,lat,longit) double *lat,*longit,xin,yin; { double Rpolar,Theta,x=xin,y=yin; if (rotate != 0) { Rpolar = sqrt(x*x + y*y); Theta = atan2(x, y) - rotate; x = Rpolar*sin(Theta); y = Rpolar*cos(Theta); } x += CenterX; y += CenterY; switch(Projection) { case Mercator: invmerc(x,y,lat,longit); break; case Lambert: invlamb(x,y,lat,longit); break; case Gall: invgall(x,y,lat,longit); break; } *longit = fixlong(*longit); #ifdef DEBUG fprintf(stderr,"U(%3.1f,%3.1f) <-- [%3.1f,%3.1f]\n",*lat,*longit,xin,yin); #endif DEBUG } /* ------------------------------------------------------------------------ */ /* This function is called from init to print an identification string for the map caption. */ identify() { switch(Projection) { case Mercator: identMerc(); break; case Lambert: identLamb(); break; case Gall: identGall(); break; } } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Lambert Projection: the paper is a cone whose axis is vertical and which intersects the sphere at two places, called Lam1 and Lam2. If Lam1=Lam2 it is an ordinary conic projection. */ /* These are precomputed intermediate values */ static double L, /* half-ht of lambert zone in radians */ alpha, /* latitude of L w.r.t. equator */ sinalpha, cosalpha, ralpha, /* length of perp from center to L */ H; /* height of cone */ #define R 1.0 /* radius of the sphere */ initLamb() { if (Lamb1==0 || Lamb2==0) { fputs("LC error: -l option missing, or bad values.\n",stderr); exit(2); } if (abs(Lamb2) > abs(Lamb1)) {double T; T=Lamb1; Lamb1=Lamb2; Lamb2=T; } /* Precompute as much as we can for the lambert conversion routine */ L = rads((Lamb1 - Lamb2)/2); alpha = (rads(Lamb2) + L); ralpha = R * cos(L); sinalpha = sin(alpha); cosalpha = cos(alpha); H = ralpha / sinalpha; } identLamb() { printf("Lambert Conformal Projection [%s,",prll(Lamb1,1,1)); printf("%s]",prll(Lamb2,1,1)); } /* ---------------------------------------------------------------------- This function converts "inlat" and "inlong" into Lambert Conformal coordinates that are stored into "xout" and "yout". The Lambert cone intersects the earth at latitude Lam1 and Lam2; we are assuming always that Lam1 is farther north than Lam2. The map is centered on LongCenter. */ lambert(inlat,inlong,xout,yout) double inlat,inlong, *xout, *yout; { double beta, /* relative lat of inlat w.r.t. alpha */ rcone, /* radius to cone at [lat,long] intersec */ PolarR, /* distance down cone from vertex */ PolarTheta; /* relative longitude */ if (sign(inlat)==sign(Lamb1)) { beta = rads(inlat) - alpha; PolarR = H * cosalpha - (ralpha * tan(beta)); rcone = PolarR*sinalpha; PolarTheta = (rads(fixlong(inlong-LongCenter)))*(rcone/PolarR); *xout = MapScale*PolarR*sin(PolarTheta); *yout = - MapScale*(PolarR * cos(PolarTheta)); } else { beta = rads(-inlat) - alpha; PolarR = H * cosalpha - (ralpha * tan(beta)); rcone = PolarR*sinalpha; PolarTheta = (rads(fixlong(inlong-LongCenter)))*(rcone/PolarR); *xout = MapScale*PolarR*sin(PolarTheta); *yout = MapScale*((PolarR * cos(PolarTheta))-2.0*H); } } /* ---------------------------------------------------------------------- Do an inverse Lambert transformation. We need this so that we can get a rough approximation to the clipping region in lat/long space, so that we don't have to transform every point to see if we need to clip it. */ invlamb(inx,iny,latout,longout) double inx,iny,*latout,*longout; { double Rpolar,Theta,Phi,beta,xx,yy,rcone; xx = inx/MapScale; yy = -iny/MapScale; Rpolar = sqrt(xx*xx + yy*yy); if (H<0) Rpolar = -Rpolar; /* get southern hemisphere right */ rcone = Rpolar*sinalpha; Theta = atan2(xx, yy); if (Theta > halfpi) Theta -= pi; if (Theta < -halfpi) Theta += pi; beta = atan((H*cosalpha-Rpolar)/ralpha); Phi = alpha + beta; *latout = degs(Phi); *longout = LongCenter + degs(Theta*abs(Rpolar)/rcone); } /* ----------------- End of Lambert projection code ----------------------- */ /* ------------------- Mercator Projection -------------------------------- */ /* Mercator projection: the paper is a cylinder whose axis is vertical and that is tangent to the sphere at the equator. */ initMerc() { } identMerc() { printf("Mercator projection"); } mercat(inlat,inlong,xout,yout) double inlat,inlong, *xout, *yout; { *xout = MapScale*R*(rads(fixlong(inlong - LongCenter))); *yout = MapScale*R*tan(rads(inlat)); } invmerc(inx,iny,latout,longout) double inx,iny,*latout,*longout; { *longout = LongCenter + degs(inx/(R*MapScale)); *latout = degs(atan2(iny,R*MapScale)); } /* ---------------------- Gall Projection --------------------------------- */ /* Gall Stereographic projection: the paper is a cylinder whose axis is vertical and that intersects the sphere at a latitude of 45 degrees. */ static double RX,Ry; initGall() { RX = R/sqrt(2.0); Ry = R + RX; } identGall() { printf("Gall Stereographic Projection"); } gallpr(inlat,inlong,xout,yout) double inlat,inlong, *xout, *yout; { *xout = MapScale*RX*(rads(fixlong(inlong - LongCenter))); *yout = MapScale*(Ry*sin(rads(inlat))/(1.0 + cos(rads(inlat)))); } invgall(inx,iny,latout,longout) double inx,iny,*latout,*longout; { double alpha; *longout = LongCenter + degs(inx/(RX*MapScale)); iny = iny/MapScale; alpha = atan2(iny, Ry); if (abs(iny) > Ry) *latout = 90*sign(iny); else *latout = degs(asin(2*cos(alpha)*sin(alpha))); }