| 1 | |
|---|
| 2 | |
|---|
| 3 | def degminsecs2decimal_degrees(dd,mm,ss): |
|---|
| 4 | assert abs(mm) == mm |
|---|
| 5 | assert abs(ss) == ss |
|---|
| 6 | |
|---|
| 7 | if dd < 0: |
|---|
| 8 | sign = -1 |
|---|
| 9 | else: |
|---|
| 10 | sign = 1 |
|---|
| 11 | |
|---|
| 12 | return sign * (abs(dd) + mm/60. + ss/3600.) |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | def redfearn(lat, lon): |
|---|
| 16 | """Compute UTM projection using Redfearn's formula |
|---|
| 17 | |
|---|
| 18 | lat, lon is latitude and longitude in decimal degrees |
|---|
| 19 | """ |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | from math import pi, sqrt, sin, cos, tan |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | #Convert to radians |
|---|
| 26 | phi = lat*pi/180 |
|---|
| 27 | lam = lon*pi/180 |
|---|
| 28 | |
|---|
| 29 | #GDA Specifications |
|---|
| 30 | a = 6378137.0 #Semi major axis |
|---|
| 31 | inverse_flattening = 298.257222101 #1/f |
|---|
| 32 | zone_width = 6 #Degrees |
|---|
| 33 | false_easting = 500000 |
|---|
| 34 | false_northing = 10000000 #In the southern hemisphere (0 in the northern) |
|---|
| 35 | longitude_of_central_meridian_zone1 = -177 #Needed??? |
|---|
| 36 | longitude_of_central_meridian_zone0 = -183 |
|---|
| 37 | longitude_of_western_edge_zone0 = -186 |
|---|
| 38 | K0 = 0.9996 #Central scale factor |
|---|
| 39 | |
|---|
| 40 | #Derived constants |
|---|
| 41 | f = 1.0/inverse_flattening |
|---|
| 42 | b = a*(1-f) #Semi minor axis |
|---|
| 43 | |
|---|
| 44 | e2 = 2*f - f*f# = f*(2-f) = (a^2-b^2/a^2 #Eccentricity |
|---|
| 45 | e = sqrt(e2) |
|---|
| 46 | e2_ = e2/(1-e2) # = (a^2-b^2)/b^2 #Second eccentricity |
|---|
| 47 | e_ = sqrt(e2_) |
|---|
| 48 | e4 = e2*e2 |
|---|
| 49 | e6 = e2*e4 |
|---|
| 50 | |
|---|
| 51 | #Foot point latitude |
|---|
| 52 | n = (a-b)/(a+b) #Same as e2 - why ? |
|---|
| 53 | n2 = n*n |
|---|
| 54 | n3 = n*n2 |
|---|
| 55 | n4 = n2*n2 |
|---|
| 56 | |
|---|
| 57 | G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180 |
|---|
| 58 | |
|---|
| 59 | sinphi = sin(phi) #E9 |
|---|
| 60 | sin2phi = sin(2*phi) |
|---|
| 61 | sin4phi = sin(4*phi) |
|---|
| 62 | sin6phi = sin(6*phi) |
|---|
| 63 | |
|---|
| 64 | cosphi = cos(phi) #E22 |
|---|
| 65 | cosphi2 = cosphi*cosphi |
|---|
| 66 | cosphi3 = cosphi*cosphi2 |
|---|
| 67 | cosphi4 = cosphi2*cosphi2 |
|---|
| 68 | cosphi5 = cosphi*cosphi4 |
|---|
| 69 | cosphi6 = cosphi2*cosphi4 |
|---|
| 70 | cosphi7 = cosphi*cosphi6 |
|---|
| 71 | cosphi8 = cosphi4*cosphi4 |
|---|
| 72 | |
|---|
| 73 | t = tan(phi) |
|---|
| 74 | t2 = t*t |
|---|
| 75 | t4 = t2*t2 |
|---|
| 76 | t6 = t2*t4 |
|---|
| 77 | |
|---|
| 78 | #Radius of Curvature |
|---|
| 79 | rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5 |
|---|
| 80 | nu = a/(1-e2*sinphi*sinphi)**0.5 |
|---|
| 81 | psi = nu/rho |
|---|
| 82 | psi2 = psi*psi |
|---|
| 83 | psi3 = psi*psi2 |
|---|
| 84 | psi4 = psi2*psi2 |
|---|
| 85 | |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | #Meridian distance |
|---|
| 89 | |
|---|
| 90 | A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256 |
|---|
| 91 | A2 = 3.0/8*(e2+e4/4+15*e6/128) |
|---|
| 92 | A4 = 15.0/256*(e4+3*e6/4) |
|---|
| 93 | A6 = 35*e6/3072 |
|---|
| 94 | |
|---|
| 95 | term1 = a*A0*phi |
|---|
| 96 | term2 = -a*A2*sin2phi |
|---|
| 97 | term3 = a*A4*sin4phi |
|---|
| 98 | term4 = -a*A6*sin6phi |
|---|
| 99 | |
|---|
| 100 | m = term1 + term2 + term3 + term4 #OK |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | zone_r = (lon - longitude_of_western_edge_zone0)/zone_width |
|---|
| 104 | zone = int(zone_r) |
|---|
| 105 | central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 |
|---|
| 106 | |
|---|
| 107 | omega = (lon-central_meridian)*pi/180 #Relative longitude (radians) |
|---|
| 108 | omega2 = omega*omega |
|---|
| 109 | omega3 = omega*omega2 |
|---|
| 110 | omega4 = omega2*omega2 |
|---|
| 111 | omega5 = omega*omega4 |
|---|
| 112 | omega6 = omega3*omega3 |
|---|
| 113 | omega7 = omega*omega6 |
|---|
| 114 | omega8 = omega4*omega4 |
|---|
| 115 | |
|---|
| 116 | #Northing |
|---|
| 117 | term1 = nu*sinphi*cosphi*omega2/2 |
|---|
| 118 | term2 = nu*sinphi*cosphi3*omega4/24*(4*psi2+psi-t2) |
|---|
| 119 | term3 = nu*sinphi*cosphi5*omega6/720*\ |
|---|
| 120 | (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\ |
|---|
| 121 | psi2*(1-32*t2)-psi*2*t2+t4-t2) |
|---|
| 122 | term4 = nu*sinphi*cosphi7*omega8/40320*(1385-3111*t2+543*t4-t6) |
|---|
| 123 | sum = m + term1 + term2 + term3 + term4 |
|---|
| 124 | northing = false_northing + K0*sum |
|---|
| 125 | |
|---|
| 126 | print northing, zone |
|---|
| 127 | |
|---|
| 128 | import sys; sys.exit() |
|---|
| 129 | |
|---|
| 130 | return zone, easting, northing |
|---|
