paul@hpldola.HP.COM (Paul Bame) (03/22/89)
To satisfy both those who wanted the USGS map projections book reference and
the Albers projection code, here is that code. The first comment contains
the reference to the USGS book which may be obtained (it's pretty cheap
even) from either of the two main USGS places (I got it from Denver). This
code could stand to be re-arranged a bit but should suffice for most uses
as-is. No warranty implied or expressed by me or Hewlett-Packard - but note
that it passed the test cases in the book - for me anyway. Enjoy,
-Paul Bame
hplabs!hpldola!paul
paul@hpldola.hp.com
#!/bin/sh
# This is a shell archive. Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by paul at hpldola on Tue Mar 21 09:33:29 1989
#
# This archive contains:
# albers.c mapconstants.h
#
echo x - albers.c
cat >albers.c <<'@EOF'
#include <stdio.h>
#include <math.h>
#define CLARKE1866
#include "mapconstants.h"
/*
* Albers Conic Equal-Area Projection
*
* Formulae taken from "Map Projections Used by the U.S. Geological Survey"
* Bulletin #1532
*
* Many variable names taken from the reference.
*
* To use, call albers_setup() once and then albers_project() for each
* coordinate pair of interest.
*
* AUTHOR: Paul Bame hplabs!hpldola!paul paul@hpldola.hp.com 719-590-5557
*
*/
static double eradius, middlelon, bigC, cone_const, r0;
static
calc_q_msq(lat, qp, msqp)
double lat;
double *qp;
double *msqp;
/*
* Given latitude, calculate 'q' [eq. 3-12] and if msqp is != NULL, m^2
* [eq. 12-15].
*/
{
double s, c, es;
s = sin(lat);
es = s * CONST_Ec;
*qp = (1.0 - CONST_Ecsq) * ((s / (1 - es * es)) -
(1 / (2 * CONST_Ec)) * log((1 - es) / (1 + es)));
if (msqp != NULL)
{
c = cos(lat);
*msqp = c * c / (1 - es * es);
}
}
albers_setup(x_eradius, southlat, northlat, middlelat, x_middlelon)
double x_eradius;
double southlat, northlat;
double middlelat;
double x_middlelon;
/*
* x_eradius is the desired radius of the Earth in any units
*
* *lat are the standard parallels used for the Albers projection in radians
*
* x_middlelon is the longitude of the center of the projection
*/
{
double q1, q2, q0;
double m1sq, m2sq;
middlelon = x_middlelon;
eradius = x_eradius;
calc_q_msq(southlat, &q1, &m1sq);
calc_q_msq(northlat, &q2, &m2sq);
calc_q_msq(middlelat, &q0, NULL);
cone_const = (m1sq - m2sq) / (q2 - q1);
bigC = m1sq + cone_const * q1;
r0 = eradius * sqrt(bigC - cone_const * q0) / cone_const;
}
albers_project(lon, lat, xp, yp)
double lon, lat;
double *xp, *yp;
/*
* Project lon/lat (in radians) according to albers_setup and return the
* results via xp, yp.
*/
{
double q, r, theta;
calc_q_msq(lat, &q, NULL);
theta = cone_const * (lon - middlelon);
r = eradius * sqrt(bigC - cone_const * q) / cone_const;
*xp = r * sin(theta);
*yp = r0 - r * cos(theta);
}
@EOF
chmod 644 albers.c
echo x - mapconstants.h
cat >mapconstants.h <<'@EOF'
/*
* This somewhat obtuse header file gives constants required for the
* Albers projection (Ec and Ecsq) and several Earth radii from which
* to choose.
*
* CLARKE1866 means to use the Clarke 1866 Earth spheroid which is often
* used by the USGS - there are other spheroids however - see the book.
*
* AUTHOR: Paul Bame hplabs!hpldola!paul paul@hpldola.hp.com 719-590-5557
*/
#ifdef XXX
#define CONST_Ecsq (0.006695662) /* */
#define CONST_Ec (0.081827028) /* */
#endif
#ifdef CLARKE1866
#define CONST_EradiusKM (6378.2064) /* Kilometers */
#define CONST_Eradius (CONST_EradiusKM/1.609)
#define CONST_Ecsq (0.006768658) /* */
#define CONST_Ec (0.082271854) /* */
#endif
#ifndef CONST_Ecsq
#define CONST_Ecsq (CONST_Ec * CONST_Ec)
#endif
#define CONST_Eradiussq (CONST_Eradius * CONST_Eradius)
@EOF
chmod 666 mapconstants.h
exit 0