[6149] | 1 | """Implementation of Redfearn's formula to compute UTM projections from latitude and longitude |
---|
| 2 | |
---|
| 3 | Based in part on spreadsheet |
---|
| 4 | www.icsm.gov.au/gda/gdatm/redfearn.xls |
---|
| 5 | downloaded from INTERGOVERNMENTAL COMMITTEE ON SURVEYING & MAPPING (ICSM) |
---|
| 6 | http://www.icsm.gov.au/icsm/ |
---|
| 7 | |
---|
| 8 | """ |
---|
| 9 | from anuga.coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | def degminsec2decimal_degrees(dd,mm,ss): |
---|
| 13 | assert abs(mm) == mm |
---|
| 14 | assert abs(ss) == ss |
---|
| 15 | |
---|
| 16 | if dd < 0: |
---|
| 17 | sign = -1 |
---|
| 18 | else: |
---|
| 19 | sign = 1 |
---|
| 20 | |
---|
| 21 | return sign * (abs(dd) + mm/60. + ss/3600.) |
---|
| 22 | |
---|
| 23 | def decimal_degrees2degminsec(dec): |
---|
| 24 | |
---|
| 25 | if dec < 0: |
---|
| 26 | sign = -1 |
---|
| 27 | else: |
---|
| 28 | sign = 1 |
---|
| 29 | |
---|
| 30 | dec = abs(dec) |
---|
| 31 | dd = int(dec) |
---|
| 32 | f = dec-dd |
---|
| 33 | |
---|
| 34 | mm = int(f*60) |
---|
| 35 | ss = (f*60-mm)*60 |
---|
| 36 | |
---|
| 37 | return sign*dd, mm, ss |
---|
| 38 | |
---|
[6553] | 39 | def redfearn(lat, lon, false_easting=None, false_northing=None, |
---|
| 40 | zone=None, central_meridian=None, scale_factor=None): |
---|
[6149] | 41 | """Compute UTM projection using Redfearn's formula |
---|
| 42 | |
---|
| 43 | lat, lon is latitude and longitude in decimal degrees |
---|
| 44 | |
---|
| 45 | If false easting and northing are specified they will override |
---|
| 46 | the standard |
---|
| 47 | |
---|
| 48 | If zone is specified reproject lat and long to specified zone instead of |
---|
| 49 | standard zone |
---|
[6553] | 50 | If meridian is specified, reproject lat and lon to that instead of zone. In this case |
---|
| 51 | zone will be set to -1 to indicate non-UTM projection |
---|
| 52 | |
---|
| 53 | Note that zone and meridian cannot both be specifed |
---|
[6149] | 54 | """ |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | from math import pi, sqrt, sin, cos, tan |
---|
| 58 | |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | #GDA Specifications |
---|
| 62 | a = 6378137.0 #Semi major axis |
---|
| 63 | inverse_flattening = 298.257222101 #1/f |
---|
[6553] | 64 | if scale_factor is None: |
---|
| 65 | K0 = 0.9996 #Central scale factor |
---|
| 66 | else: |
---|
| 67 | K0 = scale_factor |
---|
| 68 | #print 'scale', K0 |
---|
[6149] | 69 | zone_width = 6 #Degrees |
---|
| 70 | |
---|
| 71 | longitude_of_central_meridian_zone0 = -183 |
---|
| 72 | longitude_of_western_edge_zone0 = -186 |
---|
| 73 | |
---|
| 74 | if false_easting is None: |
---|
| 75 | false_easting = 500000 |
---|
| 76 | |
---|
| 77 | if false_northing is None: |
---|
| 78 | if lat < 0: |
---|
| 79 | false_northing = 10000000 #Southern hemisphere |
---|
| 80 | else: |
---|
| 81 | false_northing = 0 #Northern hemisphere) |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | #Derived constants |
---|
| 85 | f = 1.0/inverse_flattening |
---|
| 86 | b = a*(1-f) #Semi minor axis |
---|
| 87 | |
---|
| 88 | e2 = 2*f - f*f# = f*(2-f) = (a^2-b^2/a^2 #Eccentricity |
---|
| 89 | e = sqrt(e2) |
---|
| 90 | e2_ = e2/(1-e2) # = (a^2-b^2)/b^2 #Second eccentricity |
---|
| 91 | e_ = sqrt(e2_) |
---|
| 92 | e4 = e2*e2 |
---|
| 93 | e6 = e2*e4 |
---|
| 94 | |
---|
| 95 | #Foot point latitude |
---|
| 96 | n = (a-b)/(a+b) #Same as e2 - why ? |
---|
| 97 | n2 = n*n |
---|
| 98 | n3 = n*n2 |
---|
| 99 | n4 = n2*n2 |
---|
| 100 | |
---|
| 101 | G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180 |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | phi = lat*pi/180 #Convert latitude to radians |
---|
| 105 | |
---|
| 106 | sinphi = sin(phi) |
---|
| 107 | sin2phi = sin(2*phi) |
---|
| 108 | sin4phi = sin(4*phi) |
---|
| 109 | sin6phi = sin(6*phi) |
---|
| 110 | |
---|
| 111 | cosphi = cos(phi) |
---|
| 112 | cosphi2 = cosphi*cosphi |
---|
| 113 | cosphi3 = cosphi*cosphi2 |
---|
| 114 | cosphi4 = cosphi2*cosphi2 |
---|
| 115 | cosphi5 = cosphi*cosphi4 |
---|
| 116 | cosphi6 = cosphi2*cosphi4 |
---|
| 117 | cosphi7 = cosphi*cosphi6 |
---|
| 118 | cosphi8 = cosphi4*cosphi4 |
---|
| 119 | |
---|
| 120 | t = tan(phi) |
---|
| 121 | t2 = t*t |
---|
| 122 | t4 = t2*t2 |
---|
| 123 | t6 = t2*t4 |
---|
| 124 | |
---|
| 125 | #Radius of Curvature |
---|
| 126 | rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5 |
---|
| 127 | nu = a/(1-e2*sinphi*sinphi)**0.5 |
---|
| 128 | psi = nu/rho |
---|
| 129 | psi2 = psi*psi |
---|
| 130 | psi3 = psi*psi2 |
---|
| 131 | psi4 = psi2*psi2 |
---|
| 132 | |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | #Meridian distance |
---|
| 136 | |
---|
| 137 | A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256 |
---|
| 138 | A2 = 3.0/8*(e2+e4/4+15*e6/128) |
---|
| 139 | A4 = 15.0/256*(e4+3*e6/4) |
---|
| 140 | A6 = 35*e6/3072 |
---|
| 141 | |
---|
| 142 | term1 = a*A0*phi |
---|
| 143 | term2 = -a*A2*sin2phi |
---|
| 144 | term3 = a*A4*sin4phi |
---|
| 145 | term4 = -a*A6*sin6phi |
---|
| 146 | |
---|
| 147 | m = term1 + term2 + term3 + term4 #OK |
---|
| 148 | |
---|
[6553] | 149 | if zone is not None and central_meridian is not None: |
---|
| 150 | msg = 'You specified both zone and central_meridian. Provide only one of them' |
---|
| 151 | raise Exception, msg |
---|
| 152 | |
---|
| 153 | # Zone |
---|
[6149] | 154 | if zone is None: |
---|
| 155 | zone = int((lon - longitude_of_western_edge_zone0)/zone_width) |
---|
| 156 | |
---|
[6553] | 157 | # Central meridian |
---|
| 158 | if central_meridian is None: |
---|
| 159 | central_meridian = zone*zone_width+longitude_of_central_meridian_zone0 |
---|
| 160 | else: |
---|
| 161 | zone = -1 |
---|
[6149] | 162 | |
---|
| 163 | omega = (lon-central_meridian)*pi/180 #Relative longitude (radians) |
---|
| 164 | omega2 = omega*omega |
---|
| 165 | omega3 = omega*omega2 |
---|
| 166 | omega4 = omega2*omega2 |
---|
| 167 | omega5 = omega*omega4 |
---|
| 168 | omega6 = omega3*omega3 |
---|
| 169 | omega7 = omega*omega6 |
---|
| 170 | omega8 = omega4*omega4 |
---|
| 171 | |
---|
| 172 | #Northing |
---|
| 173 | term1 = nu*sinphi*cosphi*omega2/2 |
---|
| 174 | term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24 |
---|
| 175 | term3 = nu*sinphi*cosphi5*\ |
---|
| 176 | (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\ |
---|
| 177 | psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720 |
---|
| 178 | term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320 |
---|
| 179 | northing = false_northing + K0*(m + term1 + term2 + term3 + term4) |
---|
| 180 | |
---|
| 181 | #Easting |
---|
| 182 | term1 = nu*omega*cosphi |
---|
| 183 | term2 = nu*cosphi3*(psi-t2)*omega3/6 |
---|
| 184 | term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120 |
---|
| 185 | term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040 |
---|
| 186 | easting = false_easting + K0*(term1 + term2 + term3 + term4) |
---|
| 187 | |
---|
| 188 | return zone, easting, northing |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | def convert_from_latlon_to_utm(points=None, |
---|
| 193 | latitudes=None, |
---|
| 194 | longitudes=None, |
---|
| 195 | false_easting=None, |
---|
| 196 | false_northing=None): |
---|
| 197 | """Convert latitude and longitude data to UTM as a list of coordinates. |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | Input |
---|
| 201 | |
---|
| 202 | points: list of points given in decimal degrees (latitude, longitude) or |
---|
| 203 | latitudes: list of latitudes and |
---|
| 204 | longitudes: list of longitudes |
---|
| 205 | false_easting (optional) |
---|
| 206 | false_northing (optional) |
---|
| 207 | |
---|
| 208 | Output |
---|
| 209 | |
---|
| 210 | points: List of converted points |
---|
| 211 | zone: Common UTM zone for converted points |
---|
| 212 | |
---|
| 213 | |
---|
| 214 | Notes |
---|
| 215 | |
---|
| 216 | Assume the false_easting and false_northing are the same for each list. |
---|
| 217 | If points end up in different UTM zones, an ANUGAerror is thrown. |
---|
| 218 | """ |
---|
| 219 | |
---|
| 220 | old_geo = Geo_reference() |
---|
| 221 | utm_points = [] |
---|
| 222 | if points == None: |
---|
| 223 | assert len(latitudes) == len(longitudes) |
---|
| 224 | points = map(None, latitudes, longitudes) |
---|
| 225 | |
---|
| 226 | for point in points: |
---|
| 227 | |
---|
| 228 | zone, easting, northing = redfearn(float(point[0]), |
---|
| 229 | float(point[1]), |
---|
| 230 | false_easting=false_easting, |
---|
| 231 | false_northing=false_northing) |
---|
| 232 | new_geo = Geo_reference(zone) |
---|
| 233 | old_geo.reconcile_zones(new_geo) |
---|
| 234 | utm_points.append([easting, northing]) |
---|
| 235 | |
---|
| 236 | return utm_points, old_geo.get_zone() |
---|