#!/usr/bin/env python # Lat Long - UTM, UTM - Lat Long conversions from warnings import warn from numpy import array, pi, sin, cos, tan, sqrt #LatLong- UTM conversion..h #definitions for lat/long to UTM and UTM to lat/lng conversions #include _deg2rad = pi / 180.0 _rad2deg = 180.0 / pi _EquatorialRadius = 2 _eccentricitySquared = 3 _ellipsoid = [ # id, Ellipsoid name, Equatorial Radius, square of eccentricity # first once is a placeholder only, To allow array indices to match id numbers [ -1, "Placeholder", 0, 0], [ 1, "Airy", 6377563, 0.00667054], [ 2, "Australian National", 6378160, 0.006694542], [ 3, "Bessel 1841", 6377397, 0.006674372], [ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372], [ 5, "Clarke 1866", 6378206, 0.006768658], [ 6, "Clarke 1880", 6378249, 0.006803511], [ 7, "Everest", 6377276, 0.006637847], [ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422], [ 9, "Fischer 1968", 6378150, 0.006693422], [ 10, "GRS 1967", 6378160, 0.006694605], [ 11, "GRS 1980", 6378137, 0.00669438], [ 12, "Helmert 1906", 6378200, 0.006693422], [ 13, "Hough", 6378270, 0.00672267], [ 14, "International", 6378388, 0.00672267], [ 15, "Krassovsky", 6378245, 0.006693422], [ 16, "Modified Airy", 6377340, 0.00667054], [ 17, "Modified Everest", 6377304, 0.006637847], [ 18, "Modified Fischer 1960", 6378155, 0.006693422], [ 19, "South American 1969", 6378160, 0.006694542], [ 20, "WGS 60", 6378165, 0.006693422], [ 21, "WGS 66", 6378145, 0.006694542], [ 22, "WGS-72", 6378135, 0.006694318], [ 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@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 #def LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long, # double &UTMNorthing, double &UTMEasting, char* UTMZone) def LLtoUTM(ReferenceEllipsoid, Lat, Long, ZoneNumber=None): #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@globalstar.com a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] eccSquaredSquared = eccSquared*eccSquared eccSquaredCubed = eccSquaredSquared*eccSquared k0 = 0.9996 # Make sure the longitude is between -180.00 .. 179.9 LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9 LatRad = Lat*_deg2rad LongRad = LongTemp*_deg2rad if ZoneNumber is None: ZoneNumber = int((LongTemp + 180)/6) + 1 if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0: ZoneNumber = 32 # Special zones for Svalbard if Lat >= 72.0 and Lat < 84.0: if LongTemp >= 0.0 and LongTemp < 9.0: ZoneNumber = 31 elif LongTemp >= 9.0 and LongTemp < 21.0: ZoneNumber = 33 elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35 elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37 LongOriginRad = _deg2rad*((ZoneNumber - 1)*6 - 180 + 3) #+3 puts origin in middle of zone #compute the UTM Zone from the latitude and longitude UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat)) eccPrimeSquared = (eccSquared)/(1-eccSquared) sinLR,cosLR,tanLR = sin(LatRad),cos(LatRad),tan(LatRad) N = a/sqrt(1-eccSquared*sinLR*sinLR) T = tanLR*tanLR C = eccPrimeSquared*cosLR*cosLR A = cosLR*(LongRad-LongOriginRad) A2 = A*A M = a*((1 - eccSquared/4 - 3*eccSquaredSquared/64 - 5*eccSquaredCubed/256)*LatRad - (3*eccSquared/8 + 3*eccSquaredSquared/32 + 45*eccSquaredCubed/1024)*sin(2*LatRad) + (15*eccSquaredSquared/256 + 45*eccSquaredCubed/1024)*sin(4*LatRad) - (35*eccSquaredCubed/3072)*sin(6*LatRad)) UTMEasting = (k0*N*(A+(1-T+C)*A2*A/6 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A2*A2*A/120) + 500000.0) UTMNorthing = (k0*(M+N*tanLR*(A2/2+(5-T+9*C+4*C*C)*A2*A2/24 + (61 -58*T +T*T +600*C -330*eccPrimeSquared)*A2*A2*A2/720))) if Lat < 0: UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere return (UTMZone, UTMEasting, UTMNorthing) def _UTMLetterDesignator(Lat): #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 if 84 >= Lat >= 72: return 'X' elif 72 > Lat >= 64: return 'W' elif 64 > Lat >= 56: return 'V' elif 56 > Lat >= 48: return 'U' elif 48 > Lat >= 40: return 'T' elif 40 > Lat >= 32: return 'S' elif 32 > Lat >= 24: return 'R' elif 24 > Lat >= 16: return 'Q' elif 16 > Lat >= 8: return 'P' elif 8 > Lat >= 0: return 'N' elif 0 > Lat >= -8: return 'M' elif -8> Lat >= -16: return 'L' elif -16 > Lat >= -24: return 'K' elif -24 > Lat >= -32: return 'J' elif -32 > Lat >= -40: return 'H' elif -40 > Lat >= -48: return 'G' elif -48 > Lat >= -56: return 'F' elif -56 > Lat >= -64: return 'E' elif -64 > Lat >= -72: return 'D' elif -72 > Lat >= -80: return 'C' else: return 'Z' # if the Latitude is outside the UTM limits #void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone, # double& Lat, double& Long ) def UTMtoLL(ReferenceEllipsoid, northing, easting, zone): #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 #Converted to Python by Russ Nelson k0 = 0.9996 a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] eccPrimeSquared = (eccSquared)/(1-eccSquared) eccSquaredSquared = eccSquared*eccSquared e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared)) e3 = e1*e1*e1 #NorthernHemisphere; //1 for northern hemispher, 0 for southern x = array([easting]).copy().ravel() - 500000.0 #remove 500,000 meter offset for longitude y = array([northing]).copy().ravel() # assume zone format = '11N' ZoneLetter = zone[-1] ZoneNumber = int(zone[:-1]) if ZoneLetter < 'N': # southern hemisphere < N 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 mu = (y/k0)/(a*(1-eccSquared/4-3*eccSquaredSquared/64-5*eccSquaredSquared*eccSquared/256)) phi1Rad = (mu + (3*e1/2-27*e3/32)*sin(2*mu) + (21*e1*e1/16-55*e3*e1/32)*sin(4*mu) +(151*e3/96)*sin(6*mu)) sinphi = sin(phi1Rad) cosphi = cos(phi1Rad) tanphi = tan(phi1Rad) Ndenom = 1-eccSquared*sinphi*sinphi N1 = a/sqrt(Ndenom) T1 = tanphi*tanphi C1 = eccPrimeSquared*cosphi*cosphi C2 = C1*C1 R1 = a*(1-eccSquared)/pow(Ndenom, 1.5) D = x/(N1*k0) D3 = D*D*D Lat = (phi1Rad - (N1*tanphi/R1)*(D*D/2-(5+3*T1+10*C1-4*C2-9*eccPrimeSquared)*D3*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C2)*D3*D3/720))*_rad2deg Long = LongOrigin + (D-(1+2*T1+C1)*D3/6+(5-2*C1+28*T1-3*C2+8*eccPrimeSquared+24*T1*T1) *D3*D*D/120)/cosphi * _rad2deg return (Lat, Long) if __name__ == '__main__': (z, e, n) = LLtoUTM(23, 40 + (6 + 18.3591/60)/60, -(88 + (13 + 9.52349/60)/60)) print(z, e, n) print(UTMtoLL(23, e, n, z))