Add:Core:Added support for utm projection

git-svn-id: http://svn.code.sf.net/p/navit/code/trunk/navit@2291 ffa7fe5e-494d-0410-b361-a75ebd5db220

1 files changed, 69 insertions, 0 deletions

diff --git a/navit/transform.c b/navit/transform.c
index a2711f6b..c336e01f 100644
--- a/navit/transform.c
+++ b/navit/transform.c
@@ -228,6 +228,7 @@ static const navit_float geo2gar_units = 1/(360.0/(1<<24));
 void
 transform_to_geo(enum projection pro, struct coord *c, struct coord_geo *g)
 {
+	int x,y,northern,zone;
 	switch (pro) {
 	case projection_mg:
 		g->lng=c->x/6371000.0/M_PI*180;
@@ -237,6 +238,17 @@ transform_to_geo(enum projection pro, struct coord *c, struct coord_geo *g)
 		g->lng=c->x*gar2geo_units;
 		g->lat=c->y*gar2geo_units;
 		break;
+	case projection_utm:
+		x=c->x;
+		y=c->y;
+		northern=y >= 0;
+		if (!northern) {
+			y+=10000000;
+		}
+		zone=(x/1000000);
+		x=x%1000000;
+		transform_utm_to_geo(x, y, zone, northern, g);
+		break;	
 	default:
 		break;
 	}
@@ -298,6 +310,63 @@ transform_cart_to_geo(struct coord_geo_cart *cart, navit_float a, navit_float b,
 
 
 void
+transform_utm_to_geo(const double UTMEasting, const double UTMNorthing, int ZoneNumber, int NorthernHemisphere, struct coord_geo *geo)
+{
+//converts UTM coords to lat/long.  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@globalstar.com
+
+	double Lat, Long;
+	double k0 = 0.99960000000000004;
+	double a = 6378137;
+	double eccSquared = 0.0066943799999999998;
+	double eccPrimeSquared;
+	double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
+	double N1, T1, C1, R1, D, M;
+	double LongOrigin;
+	double mu, phi1, phi1Rad;
+	double x, y;
+	double rad2deg = 180/M_PI;
+
+	x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
+	y = UTMNorthing;
+
+	if (!NorthernHemisphere) {
+		y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
+	}
+
+	LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
+
+	eccPrimeSquared = (eccSquared)/(1-eccSquared);
+
+	M = y / k0;
+	mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
+	phi1Rad = mu	+ (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) 
+				+ (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
+				+(151*e1*e1*e1/96)*sin(6*mu);
+	phi1 = phi1Rad*rad2deg;
+
+	N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
+	T1 = tan(phi1Rad)*tan(phi1Rad);
+	C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
+	R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
+	D = x/(N1*k0);
+
+	Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
+					+(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
+	Lat = Lat * rad2deg;
+
+	Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
+					*D*D*D*D*D/120)/cos(phi1Rad);
+	Long = LongOrigin + Long * rad2deg;
+
+	geo->lat=Lat;
+	geo->lng=Long;
+}
+
+void
 transform_datum(struct coord_geo *from, enum map_datum from_datum, struct coord_geo *to, enum map_datum to_datum)
 {
 }