o Added UTM conversion functions to RMProjection

@@ -54,4 +54,20 @@
 // forward project latitude/longitude, return meters
 - (RMProjectedPoint)coordinateToProjectedPoint:(CLLocationCoordinate2D)aLatLong;
 
+#pragma mark - UTM conversions
+
++ (void)convertCoordinate:(CLLocationCoordinate2D)coordinate
+          toUTMZoneNumber:(int *)utmZoneNumber
+            utmZoneLetter:(NSString **)utmZoneLetter   // may be NULL
+     isNorthernHemisphere:(BOOL *)isNorthernHemisphere // may be NULL
+                  easting:(double *)easting
+                 northing:(double *)northing;
+
++ (void)convertUTMZoneNumber:(int)utmZoneNumber
+               utmZoneLetter:(NSString *)utmZoneLetter  // #utmZoneLetter will be used for calculations if not nil,
+        isNorthernHemisphere:(BOOL)isNorthernHemisphere // otherwise #isNorthernHemisphere
+                     easting:(double)easting
+                    northing:(double)northing
+                toCoordinate:(CLLocationCoordinate2D *)coordinate;
+
 @end

@@ -45,6 +45,45 @@
 @synthesize planetBounds = _planetBounds;
 @synthesize projectionWrapsHorizontally = _projectionWrapsHorizontally;
 
+#pragma mark - Common projections
+
+static RMProjection *_googleProjection = nil;
+static RMProjection *_latitudeLongitudeProjection = nil;
+
++ (RMProjection *)googleProjection
+{
+    if (_googleProjection)
+    {
+        return _googleProjection;
+    }
+    else
+    {
+        RMProjectedRect theBounds = RMProjectedRectMake(-20037508.34, -20037508.34, 20037508.34 * 2, 20037508.34 * 2);
+
+        _googleProjection = [[RMProjection alloc] initWithString:@"+title= Google Mercator EPSG:900913 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
+                                                        inBounds:theBounds];
+        return _googleProjection;
+    }
+}
+
++ (RMProjection *)EPSGLatLong
+{
+    if (_latitudeLongitudeProjection)
+    {
+        return _latitudeLongitudeProjection;
+    }
+    else
+    {
+        RMProjectedRect theBounds = RMProjectedRectMake(-kMaxLong, -kMaxLat, 360.0, kMaxLong);
+
+        _latitudeLongitudeProjection = [[RMProjection alloc] initWithString:@"+proj=latlong +ellps=WGS84"
+                                                                   inBounds:theBounds];
+        return _latitudeLongitudeProjection;
+    }
+}
+
+#pragma mark -
+
 - (id)initWithString:(NSString *)proj4String inBounds:(RMProjectedRect)projectedBounds
 {
     if (!(self = [super init]))
@@ -149,39 +188,187 @@
     return result_coordinate;
 }
 
-static RMProjection *_googleProjection = nil;
-static RMProjection *_latitudeLongitudeProjection = nil;
+#pragma mark - UTM conversions
 
-+ (RMProjection *)googleProjection
+// This uses code by Chuck Gantz, found at http://www.gpsy.com/gpsinfo/geotoutm/
+// It is limited to WGS84, have a look at the original source code if you need more.
+
+//
+// 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
+//
+
+#define deg2rad (M_PI / 180.0)
+#define rad2deg (180.0 / M_PI)
+
+// This routine determines the correct UTM letter designator for the given latitude.
+// Returns 'Z' if latitude is outside the UTM limits of 84N to 80S
+// Written by Chuck Gantz- chuck.gantz@globalstar.com
++ (NSString *)UTMLetterDesignatorForLatitude:(double)latitude
 {
-    if (_googleProjection)
+    char letterDesignator;
+
+    if ((84 >= latitude) && (latitude >= 72)) letterDesignator = 'X';
+    else if ((72 > latitude) && (latitude >= 64)) letterDesignator = 'W';
+    else if ((64 > latitude) && (latitude >= 56)) letterDesignator = 'V';
+    else if ((56 > latitude) && (latitude >= 48)) letterDesignator = 'U';
+    else if ((48 > latitude) && (latitude >= 40)) letterDesignator = 'T';
+    else if ((40 > latitude) && (latitude >= 32)) letterDesignator = 'S';
+    else if ((32 > latitude) && (latitude >= 24)) letterDesignator = 'R';
+    else if ((24 > latitude) && (latitude >= 16)) letterDesignator = 'Q';
+    else if ((16 > latitude) && (latitude >= 8)) letterDesignator = 'P';
+    else if (( 8 > latitude) && (latitude >= 0)) letterDesignator = 'N';
+    else if (( 0 > latitude) && (latitude >= -8)) letterDesignator = 'M';
+    else if ((-8> latitude) && (latitude >= -16)) letterDesignator = 'L';
+    else if ((-16 > latitude) && (latitude >= -24)) letterDesignator = 'K';
+    else if ((-24 > latitude) && (latitude >= -32)) letterDesignator = 'J';
+    else if ((-32 > latitude) && (latitude >= -40)) letterDesignator = 'H';
+    else if ((-40 > latitude) && (latitude >= -48)) letterDesignator = 'G';
+    else if ((-48 > latitude) && (latitude >= -56)) letterDesignator = 'F';
+    else if ((-56 > latitude) && (latitude >= -64)) letterDesignator = 'E';
+    else if ((-64 > latitude) && (latitude >= -72)) letterDesignator = 'D';
+    else if ((-72 > latitude) && (latitude >= -80)) letterDesignator = 'C';
+    else letterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
+
+    return [NSString stringWithFormat:@"%c", letterDesignator];
+}
+
+// Converts latitude/longitude to UTM coordinates.  Equations from USGS Bulletin 1532.
+// East longitudes are positive, West longitudes are negative.
+// North latitudes are positive, South latitudes are negative.
+// Latitude and longitude are in decimal degrees.
+// Written by Chuck Gantz - chuck.gantz@globalstar.com
++ (void)convertCoordinate:(CLLocationCoordinate2D)coordinate
+          toUTMZoneNumber:(int *)utmZoneNumber
+            utmZoneLetter:(NSString **)utmZoneLetter
+     isNorthernHemisphere:(BOOL *)isNorthernHemisphere
+                  easting:(double *)easting
+                 northing:(double *)northing
+{
+    double a = 6378137.0;
+    double eccSquared = 0.00669438;
+    double k0 = 0.9996;
+
+    double longitudeOrigin, longitudeOriginRad;
+    double eccPrimeSquared;
+    double N, T, C, A, M;
+
+    // Make sure the longitude is between -180.00 .. 179.9
+    double longitudeTemp = (coordinate.longitude + 180.0) - (floor((coordinate.longitude + 180.0) / 360.0) * 360.0) - 180.0;
+    double latitudeRad = coordinate.latitude * deg2rad;
+    double longitudeRad = longitudeTemp * deg2rad;
+
+    *utmZoneNumber = floor((longitudeTemp + 180.0) / 6.0) + 1;
+
+    if (coordinate.latitude >= 56.0 && coordinate.latitude < 64.0 && longitudeTemp >= 3.0 && longitudeTemp < 12.0)
+        *utmZoneNumber = 32;
+
+    // Special zones for Svalbard
+    if (coordinate.latitude >= 72.0 && coordinate.latitude < 84.0)
     {
-        return _googleProjection;
+        if (longitudeTemp >= 0.0 && longitudeTemp < 9.0) *utmZoneNumber = 31;
+        else if (longitudeTemp >= 9.0 && longitudeTemp < 21.0) *utmZoneNumber = 33;
+        else if (longitudeTemp >= 21.0 && longitudeTemp < 33.0) *utmZoneNumber = 35;
+        else if (longitudeTemp >= 33.0 && longitudeTemp < 42.0) *utmZoneNumber = 37;
     }
-    else
-    {
-        RMProjectedRect theBounds = RMProjectedRectMake(-20037508.34, -20037508.34, 20037508.34 * 2, 20037508.34 * 2);
 
-        _googleProjection = [[RMProjection alloc] initWithString:@"+title= Google Mercator EPSG:900913 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
-                                                        inBounds:theBounds];
-        return _googleProjection;
+    longitudeOrigin = (*utmZoneNumber - 1) * 6 - 180 + 3;  //+3 puts origin in middle of zone
+    longitudeOriginRad = longitudeOrigin * deg2rad;
+
+    // Compute the UTM Zone from the latitude and longitude
+    NSString *utmLetterDesignator = [self UTMLetterDesignatorForLatitude:coordinate.latitude];
+
+    if (utmZoneLetter != NULL)
+        *utmZoneLetter = utmLetterDesignator;
+
+    if (isNorthernHemisphere != NULL)
+    {
+        char zoneLetterChar = [utmLetterDesignator UTF8String][0];
+        *isNorthernHemisphere = (zoneLetterChar >= 'N' && zoneLetterChar <= 'X');
     }
+
+    eccPrimeSquared = (eccSquared) / (1.0 - eccSquared);
+
+    N = a / sqrt(1.0 - eccSquared * sin(latitudeRad) * sin(longitudeRad));
+    T = tan(latitudeRad) * tan(latitudeRad);
+    C = eccPrimeSquared * cos(latitudeRad) * cos(latitudeRad);
+    A = cos(latitudeRad) * (longitudeRad - longitudeOriginRad);
+
+    M = a * ((1.0 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256) * latitudeRad
+             - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024) * sin(2*latitudeRad)
+             + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024) * sin(4*latitudeRad)
+             - (35*eccSquared*eccSquared*eccSquared/3072) * sin(6*latitudeRad));
+
+    *easting = (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);
+
+    *northing = (double)(k0 * (M + N*tan(latitudeRad) * (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 (coordinate.latitude < 0)
+        *northing += 10000000.0; //10000000 meter offset for southern hemisphere
 }
 
-+ (RMProjection *)EPSGLatLong
+// Converts UTM coords to latitude/longitude. Equations from USGS Bulletin 1532.
+// East longitudes are positive, West longitudes are negative.
+// North latitudes are positive, South latitudes are negative.
+// Latitude and longitude are in decimal degrees.
+// Written by Chuck Gantz - chuck.gantz@globalstar.com
++ (void)convertUTMZoneNumber:(int)utmZoneNumber
+               utmZoneLetter:(NSString *)utmZoneLetter
+        isNorthernHemisphere:(BOOL)isNorthernHemisphere
+                     easting:(double)easting
+                    northing:(double)northing
+                toCoordinate:(CLLocationCoordinate2D *)coordinate
 {
-    if (_latitudeLongitudeProjection)
+    double k0 = 0.9996;
+    double a = 6378137.0;
+    double eccSquared = 0.00669438;
+
+    double eccPrimeSquared;
+    double e1 = (1 - sqrt(1-eccSquared)) / (1 + sqrt(1-eccSquared));
+    double N1, T1, C1, R1, D, M;
+    double longitudeOrigin;
+    double mu, phi1, phi1Rad;
+    double x, y, latitude, longitude;
+
+    x = easting - 500000.0; // remove 500,000 meter offset for longitude
+    y = northing;
+
+    if (utmZoneLetter != nil)
     {
-        return _latitudeLongitudeProjection;
+        char zoneLetter = [utmZoneLetter UTF8String][0];
+        if ((zoneLetter >= 'c' && zoneLetter <= 'm') || (zoneLetter >= 'C' && zoneLetter <= 'M'))
+            y -= 10000000.0; // remove 10,000,000 meter offset used for southern hemisphere
     }
-    else
+    else if ( ! isNorthernHemisphere)
     {
-        RMProjectedRect theBounds = RMProjectedRectMake(-kMaxLong, -kMaxLat, 360.0, kMaxLong);
-
-        _latitudeLongitudeProjection = [[RMProjection alloc] initWithString:@"+proj=latlong +ellps=WGS84"
-                                                                   inBounds:theBounds];
-        return _latitudeLongitudeProjection;
+            y -= 10000000.0; // remove 10,000,000 meter offset used for southern hemisphere
     }
+
+    longitudeOrigin = (utmZoneNumber - 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);
+
+    latitude = 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);
+    latitude = latitude * rad2deg;
+
+    longitude = (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);
+    longitude = longitudeOrigin + longitude * rad2deg;
+
+    (*coordinate).latitude = latitude;
+    (*coordinate).longitude = longitude;
 }
 
 @end

@@ -140,6 +140,30 @@
 	center.latitude = 47.5635;
 	center.longitude = 10.20981;
 
+//    int zoneNumber;
+//    BOOL isNorthernHemisphere;
+//    NSString *utmZone;
+//    double easting, northing;
+//
+//    [RMProjection convertCoordinate:center
+//                    toUTMZoneNumber:&zoneNumber
+//                      utmZoneLetter:&utmZone
+//               isNorthernHemisphere:&isNorthernHemisphere
+//                            easting:&easting
+//                           northing:&northing];
+//
+//    NSLog(@"{%f,%f} -> %d%@ %.0f %.0f (north: %@)", center.latitude, center.longitude, zoneNumber, utmZone, easting, northing, isNorthernHemisphere ? @"YES" : @"NO");
+//
+//    CLLocationCoordinate2D coordinate;
+//    [RMProjection convertUTMZoneNumber:zoneNumber
+//                         utmZoneLetter:utmZone
+//                  isNorthernHemisphere:isNorthernHemisphere
+//                               easting:easting
+//                              northing:northing
+//                          toCoordinate:&coordinate];
+//
+//    NSLog(@"-> {%f,%f}", coordinate.latitude, coordinate.longitude);
+
 //    [mapView zoomWithLatitudeLongitudeBoundsSouthWest:CLLocationCoordinate2DMake(47.5, 10.15) northEast:CLLocationCoordinate2DMake(47.6, 10.25) animated:NO];
 
 	[mapView setZoom:10.0];