"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude """ from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE from Numeric import array def degminsec2decimal_degrees(dd,mm,ss): assert abs(mm) == mm assert abs(ss) == ss if dd < 0: sign = -1 else: sign = 1 return sign * (abs(dd) + mm/60. + ss/3600.) def decimal_degrees2degminsec(dec): if dec < 0: sign = -1 else: sign = 1 dec = abs(dec) dd = int(dec) f = dec-dd mm = int(f*60) ss = (f*60-mm)*60 return sign*dd, mm, ss def redfearn(lat, lon, false_easting=None, false_northing=None): """Compute UTM projection using Redfearn's formula lat, lon is latitude and longitude in decimal degrees If false easting and northing are specified they will override the standard """ from math import pi, sqrt, sin, cos, tan #GDA Specifications a = 6378137.0 #Semi major axis inverse_flattening = 298.257222101 #1/f K0 = 0.9996 #Central scale factor zone_width = 6 #Degrees longitude_of_central_meridian_zone0 = -183 longitude_of_western_edge_zone0 = -186 if false_easting is None: false_easting = 500000 if false_northing is None: if lat < 0: false_northing = 10000000 #Southern hemisphere else: false_northing = 0 #Northern hemisphere) #Derived constants f = 1.0/inverse_flattening b = a*(1-f) #Semi minor axis e2 = 2*f - f*f# = f*(2-f) = (a^2-b^2/a^2 #Eccentricity e = sqrt(e2) e2_ = e2/(1-e2) # = (a^2-b^2)/b^2 #Second eccentricity e_ = sqrt(e2_) e4 = e2*e2 e6 = e2*e4 #Foot point latitude n = (a-b)/(a+b) #Same as e2 - why ? n2 = n*n n3 = n*n2 n4 = n2*n2 G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180 phi = lat*pi/180 #Convert latitude to radians sinphi = sin(phi) sin2phi = sin(2*phi) sin4phi = sin(4*phi) sin6phi = sin(6*phi) cosphi = cos(phi) cosphi2 = cosphi*cosphi cosphi3 = cosphi*cosphi2 cosphi4 = cosphi2*cosphi2 cosphi5 = cosphi*cosphi4 cosphi6 = cosphi2*cosphi4 cosphi7 = cosphi*cosphi6 cosphi8 = cosphi4*cosphi4 t = tan(phi) t2 = t*t t4 = t2*t2 t6 = t2*t4 #Radius of Curvature rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5 nu = a/(1-e2*sinphi*sinphi)**0.5 psi = nu/rho psi2 = psi*psi psi3 = psi*psi2 psi4 = psi2*psi2 #Meridian distance A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256 A2 = 3.0/8*(e2+e4/4+15*e6/128) A4 = 15.0/256*(e4+3*e6/4) A6 = 35*e6/3072 term1 = a*A0*phi term2 = -a*A2*sin2phi term3 = a*A4*sin4phi term4 = -a*A6*sin6phi m = term1 + term2 + term3 + term4 #OK #Zone zone = int((lon - longitude_of_western_edge_zone0)/zone_width) central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 omega = (lon-central_meridian)*pi/180 #Relative longitude (radians) omega2 = omega*omega omega3 = omega*omega2 omega4 = omega2*omega2 omega5 = omega*omega4 omega6 = omega3*omega3 omega7 = omega*omega6 omega8 = omega4*omega4 #Northing term1 = nu*sinphi*cosphi*omega2/2 term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24 term3 = nu*sinphi*cosphi5*\ (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\ psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720 term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320 northing = false_northing + K0*(m + term1 + term2 + term3 + term4) #Easting term1 = nu*omega*cosphi term2 = nu*cosphi3*(psi-t2)*omega3/6 term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120 term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040 easting = false_easting + K0*(term1 + term2 + term3 + term4) return zone, easting, northing # FIXME (Ole) def convert_points_from_latlon_to_utm(latitudes, longitudes, false_easting=None, false_northing=None): """Wrapper """ return convert_lats_longs(latitudes, longitudes, false_easting, false_northing) def convert_lats_longs(latitudes, longitudes, false_easting=None, false_northing=None): """ Redfearn a list of latitudes and longitudes. Assume the false_easting and false_northing are the same for each list. Note, if a zone turns out to be different, an ANUGAerror is thrown. """ # Using geo_ref so there is only one point in the code where zones # are checked. old_geo = Geo_reference() points = [] for lat, long in map(None, latitudes, longitudes,): zone, easting, northing = redfearn(float(lat), float(long), false_easting=false_easting, false_northing=false_northing) new_geo = Geo_reference(zone) old_geo.reconcile_zones(new_geo) #print "old_geo.get_zone()", old_geo.get_zone() points.append([easting, northing]) return old_geo.get_zone(), points