[mapserver-users] Coordinate Conversions

Tue Oct 2 10:28:06 EDT 2001

Kieran -

Here's some code that's been kicking around the net for a while; I've
used it just fine for the TopoZone Web site for the last two years
(NAD27 only, however, so I can't vouch for the other datums, and I don't
need results more precise than 5 meters).  It's C++ code as it stands,
but there's nothing special about it - I ported it to JavaScript in
about 10 minutes once.  If you are only interested in one specific
datum/ellipsoid, you can just grab the equatorial radius and
eccentricity squared values from the table and toss out all the
Ellipsoid class stuff, simplifying things even more.

	- Ed

Ed McNierney
Chief Mapmaker
TopoZone.com
ed at topozone.com
(978) 251-4242

const double PI = 3.14159265;
const double FOURTHPI = PI / 4;
const double deg2rad = PI / 180;
const double rad2deg = 180.0 / PI;

#define CLARKE_1866		5
#define NAD_83		   11

class Ellipsoid
{
public:
	Ellipsoid(){};
	Ellipsoid(int Id, char* name, double radius, double ecc)
	{
		id = Id; ellipsoidName = name; 
		EquatorialRadius = radius; eccentricitySquared = ecc;
	}

	int id;
	char* ellipsoidName;
	double EquatorialRadius; 
	double eccentricitySquared;  

};

static Ellipsoid ellipsoid[] = 
{//  id, Ellipsoid name, Equatorial Radius, square of eccentricity	
	Ellipsoid( -1, "Placeholder", 0, 0),//placeholder only, To allow
array indices to match id numbers
	Ellipsoid( 1, "Airy", 6377563, 0.00667054),
	Ellipsoid( 2, "Australian National", 6378160, 0.006694542),
	Ellipsoid( 3, "Bessel 1841", 6377397, 0.006674372),
	Ellipsoid( 4, "Bessel 1841 (Nambia) ", 6377484, 0.006674372),
	Ellipsoid( 5, "Clarke 1866", 6378206, 0.006768658),
	Ellipsoid( 6, "Clarke 1880", 6378249, 0.006803511),
	Ellipsoid( 7, "Everest", 6377276, 0.006637847),
	Ellipsoid( 8, "Fischer 1960 (Mercury) ", 6378166, 0.006693422),
	Ellipsoid( 9, "Fischer 1968", 6378150, 0.006693422),
	Ellipsoid( 10, "GRS 1967", 6378160, 0.006694605),
	Ellipsoid( 11, "GRS 1980", 6378137, 0.00669438),
	Ellipsoid( 12, "Helmert 1906", 6378200, 0.006693422),
	Ellipsoid( 13, "Hough", 6378270, 0.00672267),
	Ellipsoid( 14, "International", 6378388, 0.00672267),
	Ellipsoid( 15, "Krassovsky", 6378245, 0.006693422),
	Ellipsoid( 16, "Modified Airy", 6377340, 0.00667054),
	Ellipsoid( 17, "Modified Everest", 6377304, 0.006637847),
	Ellipsoid( 18, "Modified Fischer 1960", 6378155, 0.006693422),
	Ellipsoid( 19, "South American 1969", 6378160, 0.006694542),
	Ellipsoid( 20, "WGS 60", 6378165, 0.006693422),
	Ellipsoid( 21, "WGS 66", 6378145, 0.006694542),
	Ellipsoid( 22, "WGS-72", 6378135, 0.006694318),
	Ellipsoid( 23, "WGS-84", 6378137, 0.00669438)
};
/*Reference ellipsoids derived from Peter H. Dana's website- 
http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
Department of Geography, University of Texas at Austin
Internet: pdana at mail.utexas.edu
3/22/95

Source
Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to
Department of Defense World Geodetic System
1984 Technical Report. Part I and II. Washington, DC: Defense Mapping
Agency
*/

void LLtoUTM(int ReferenceEllipsoid, const double Lat, const double
Long, 
			 double &UTMNorthing, double &UTMEasting, int
&UTMZone)
{
//converts lat/long to UTM coords.  Equations from USGS Bulletin 1532 
//East Longitudes are positive, West longitudes are negative. 
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees
	//Written by Chuck Gantz- chuck.gantz at globalstar.com

	double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
	double eccSquared =
ellipsoid[ReferenceEllipsoid].eccentricitySquared;
	double k0 = 0.9996;

	double LongOrigin;
	double eccPrimeSquared;
	double N, T, C, A, M;

//Make sure the longitude is between -180.00 .. 179.9
	double LongTemp = (Long+180)-int((Long+180)/360)*360-180; //
-180.00 .. 179.9;

	double LatRad = Lat*deg2rad;
	double LongRad = LongTemp*deg2rad;
	double LongOriginRad;
	int    ZoneNumber;

	ZoneNumber = int((LongTemp + 180)/6) + 1;

	if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp <
12.0 )
		ZoneNumber = 32;

  // Special zones for Svalbard
	if( Lat >= 72.0 && Lat < 84.0 ) 
	{
	  if(      LongTemp >= 0.0  && LongTemp <  9.0 ) ZoneNumber =
31;
	  else if( LongTemp >= 9.0  && LongTemp < 21.0 ) ZoneNumber =
33;
	  else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber =
35;
	  else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber =
37;
	 }
	LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in
middle of zone
	LongOriginRad = LongOrigin * deg2rad;

	UTMZone = ZoneNumber;

	eccPrimeSquared = (eccSquared)/(1-eccSquared);

	N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
	T = tan(LatRad)*tan(LatRad);
	C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
	A = cos(LatRad)*(LongRad-LongOriginRad);

	M = a*((1	- eccSquared/4		-
3*eccSquared*eccSquared/64	-
5*eccSquared*eccSquared*eccSquared/256)*LatRad 
				- (3*eccSquared/8	+
3*eccSquared*eccSquared/32	+
45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)

+ (15*eccSquared*eccSquared/256 +
45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad) 

- (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad));

	UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
					+
(5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
					+ 500000.0);

	UTMNorthing =
(double)(k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
				 +
(61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
	if(Lat < 0)
		UTMNorthing += 10000000.0; //10000000 meter offset for
southern hemisphere
}

-----Original Message-----
From: Kieran J. Ames [mailto:kames at keyspanenergy.com]
Sent: Tuesday, October 02, 2001 9:29 AM
To: mapserver-users
Subject: [mapserver-users] Coordinate Conversions

List,
I've been struggling in an attempt to find/decipher/build a tool that I
can invoke from the command line that will convert decimal degrees to
UTM, returning Easting, Northing and Zone. I need to invoke it in a Perl
cgi process so that I can use to redirect a browser client.
eg: for argument's sake, a pgm called ll2utm would be called like this:
system "ll2utm 40.123456 -71.654321"
and I would be returned (at least) the 3 elements above.

An executable or Perl sub routine would be my dream come true!
Can anyone point me to the right place or forward such a program?
Thanks so much for any assistance.
Kieran