🚀 go-pugleaf

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)));
}