paul@hpldola.HP.COM (Paul Bame) (07/27/88)
I'm casting about for a map of the continental US which includes state
boundaries and whose coordinate system is longitude/latitude. Since I
work for HP, and HP is licensed for the Harvard versions of the World
Databank II stuff that's OK if you have it. In fact, I have the projected
version of the map I want from the Harvard stuff (called state.something) and
nobody seems to have the states.something version which is unprojected.
In TRADE, I'm willing to offer my routines to do the "Albers Equal-Area
Conic" projection which is the one typically used to project maps of
the continental US. You could look this up in PB-270 304 from NTIS/CIA
but I did this and found error(s?) so there is actually some value
added here. This projection is purportedly the one used to project
the map I have - so I could back-project it in theory - but either my code
or their projection (or both) has a problem.
This is all so I can plot cities on a map of the US with state boundaries
and connecting arcs - for a packet-radio map. I'll append the code since
it's not long - nor is it particularly beautiful - but it works and can
be beautified if desired.
If you have an appropriate map, send me mail. I have internet access too
if that helps. If I get a map which is *not* derived from the Harvard
WDBII I'd be glad to share it - maybe via comp.sources.misc if it's not
too big (and it shouldn't be *too* big). THANKS!!!
--Paul Bame
hplabs!hpldola!paul
paul@hpldola.hp.com
(719) 590-5557 [has answering machine]
P.S. The maps on simtel20.arpa/PD1:<MSDOS.WORLDMAP> are missing state lines
#!/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 Jul 26 19:36:02 1988
#
# This archive contains:
# albers.c
# mapconstants.h
#
echo x - albers.c
cat >albers.c <<'@EOF'
#include <math.h>
#include "mapconstants.h"
/*
* Albers Conic Equal-Area Projection
*/
/* Some constants. Yes, I know the log() isn't constant but I was tired */
#define csqtemp1 (0.5 / (1.0 - CONST_Ecsq))
#define csqtemp2 (0.25 / CONST_Ec * log((1.0 + CONST_Ec)/(1.0 - CONST_Ec)))
#define CSQ (CONST_Eradiussq * (1.0 - CONST_Ecsq) * (csqtemp1 + csqtemp2))
#define EDENOM (1.0 + 2.0 / 3.0 * CONST_Ecsq + 3.0 / 5.0 * CONST_Ec4 + 4.0 / 7.0 * CONST_Ec6)
static double albers_l, ksq, r0;
static double sinB1, sinB2;
static double centrallon;
static
double
calc_lx(angle)
double angle;
{
double cosasq, sinasq;
sinasq = sin(angle);
sinasq *= sinasq;
cosasq = cos(angle);
cosasq *= cosasq;
return ((CONST_Eradiussq * cosasq) / (1 - CONST_Ecsq * sinasq));
}
static
double
calc_B1B2(parallel)
double parallel;
{
double tempsin, tempsinsq, tempsin4, tempsin6, numerator;
tempsin = sin(parallel);
tempsinsq = tempsin * tempsin;
tempsin4 = tempsinsq * tempsinsq;
tempsin6 = tempsin4 * tempsinsq;
numerator = 1.0 + 2 / 3 * CONST_Ecsq * tempsinsq +
3 / 5 * CONST_Ec4 * tempsin4 +
4 / 7 * CONST_Ec6 * tempsin6;
return (sin(parallel) * (numerator / EDENOM));
}
albers_setup(sp1, sp2, x_centrallon)
double sp1, sp2;
double x_centrallon;
/*
* sp* are the standard parallels used for the Albers projection in radians
* sp1 is the southerly and sp2 the northerly
*
* centrallon is the longitude line in the center of the projection
*/
{
double temp;
double sinB0;
centrallon = x_centrallon;
sinB1 = calc_B1B2(sp1);
sinB2 = calc_B1B2(sp2);
albers_l = (calc_lx(sp1) - calc_lx(sp2)) /
(2 * CSQ * (sinB2 - sinB1));
temp = sin(sp1);
ksq = (CONST_Eradius * cos(sp1)) /
(albers_l * sqrt(1 - CONST_Ecsq * temp * temp));
ksq *= ksq;
sinB0 = calc_B1B2((sp1 + sp2) / 2.0);
r0 = sqrt(ksq + 2 * CSQ / albers_l * (sinB1 - sinB0));
}
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 rm, sinB;
sinB = calc_B1B2(lat);
rm = sqrt(ksq + (2 * CSQ / albers_l) * (sinB1 - sinB));
*xp = rm * sin(albers_l * (lon - centrallon));
*yp = r0 - rm * cos(albers_l * (lon - centrallon));
}
@EOF
chmod 644 albers.c
echo x - mapconstants.h
cat >mapconstants.h <<'@EOF'
/* Eccentricity of the Earth - from the Dover Textbook on AstroDynamics */
#define CONST_Ec (0.08182)
#define CONST_Ecsq (CONST_Ec * CONST_Ec)
#define CONST_Ec4 (CONST_Ecsq * CONST_Ecsq)
#define CONST_Ec6 (CONST_Ec4 * CONST_Ecsq)
#define CONST_Eradius (3963.195563) /* Statute Miles */
/* #define CONST_Eradius (6378.145) /* Kilometers */
/* #define CONST_Eradius (3443.922786) /* Nautical Miles */
#define CONST_Eradiussq (CONST_Eradius * CONST_Eradius)
@EOF
chmod 666 mapconstants.h
exit 0