[6149] | 1 | |
---|
| 2 | #Test of redfearns formula. Tests can be verified at |
---|
| 3 | # |
---|
| 4 | #http://www.cellspark.com/UTM.html |
---|
| 5 | #http://www.ga.gov.au/nmd/geodesy/datums/redfearn_geo_to_grid.jsp |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | import unittest |
---|
| 9 | |
---|
| 10 | from redfearn import * |
---|
| 11 | from anuga.utilities.anuga_exceptions import ANUGAError |
---|
| 12 | |
---|
| 13 | import Numeric as num |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | #------------------------------------------------------------- |
---|
| 17 | |
---|
| 18 | class TestCase(unittest.TestCase): |
---|
| 19 | |
---|
| 20 | def test_decimal_degrees_conversion(Self): |
---|
| 21 | lat = degminsec2decimal_degrees(-37,39,10.15610) |
---|
| 22 | lon = degminsec2decimal_degrees(143,55,35.38390) |
---|
| 23 | assert num.allclose(lat, -37.65282114) |
---|
| 24 | assert num.allclose(lon, 143.9264955) |
---|
| 25 | |
---|
| 26 | dd,mm,ss = decimal_degrees2degminsec(-37.65282114) |
---|
| 27 | assert dd==-37 |
---|
| 28 | assert mm==39 |
---|
| 29 | assert num.allclose(ss, 10.15610) |
---|
| 30 | |
---|
| 31 | dd,mm,ss = decimal_degrees2degminsec(143.9264955) |
---|
| 32 | assert dd==143 |
---|
| 33 | assert mm==55 |
---|
| 34 | assert num.allclose(ss, 35.38390) |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | def test_UTM_1(self): |
---|
| 38 | #latitude: -37 39' 10.15610" |
---|
| 39 | #Longitude: 143 55' 35.38390" |
---|
| 40 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 41 | #Zone: 54 |
---|
| 42 | #Easting: 758173.797 Northing: 5828674.340 |
---|
| 43 | #Latitude: -37 39 ' 10.15610 '' Longitude: 143 55 ' 35.38390 '' |
---|
| 44 | #Grid Convergence: 1 47 ' 19.36 '' Point Scale: 1.00042107 |
---|
| 45 | |
---|
| 46 | lat = degminsec2decimal_degrees(-37,39,10.15610) |
---|
| 47 | lon = degminsec2decimal_degrees(143,55,35.38390) |
---|
| 48 | assert num.allclose(lat, -37.65282114) |
---|
| 49 | assert num.allclose(lon, 143.9264955) |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | zone, easting, northing = redfearn(lat,lon) |
---|
| 53 | |
---|
| 54 | assert zone == 54 |
---|
| 55 | assert num.allclose(easting, 758173.797) |
---|
| 56 | assert num.allclose(northing, 5828674.340) |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | def test_UTM_2(self): |
---|
| 60 | #TEST 2 |
---|
| 61 | |
---|
| 62 | #Latitude: -37 57 03.7203 |
---|
| 63 | #Longitude: 144 25 29.5244 |
---|
| 64 | #Zone: 55 |
---|
| 65 | #Easting: 273741.297 Northing: 5796489.777 |
---|
| 66 | #Latitude: -37 57 ' 3.72030 '' Longitude: 144 25 ' 29.52440 '' |
---|
| 67 | #Grid Convergence: -1 35 ' 3.65 '' Point Scale: 1.00023056 |
---|
| 68 | |
---|
| 69 | lat = degminsec2decimal_degrees(-37,57,03.7203) |
---|
| 70 | lon = degminsec2decimal_degrees(144,25,29.5244) |
---|
| 71 | #print lat, lon |
---|
| 72 | |
---|
| 73 | zone, easting, northing = redfearn(lat,lon) |
---|
| 74 | |
---|
| 75 | assert zone == 55 |
---|
| 76 | assert num.allclose(easting, 273741.297) |
---|
| 77 | assert num.allclose(northing, 5796489.777) |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | def test_UTM_3(self): |
---|
| 81 | #Test 3 |
---|
| 82 | lat = degminsec2decimal_degrees(-60,0,0) |
---|
| 83 | lon = degminsec2decimal_degrees(130,0,0) |
---|
| 84 | |
---|
| 85 | zone, easting, northing = redfearn(lat,lon) |
---|
| 86 | #print zone, easting, northing |
---|
| 87 | |
---|
| 88 | assert zone == 52 |
---|
| 89 | assert num.allclose(easting, 555776.267) |
---|
| 90 | assert num.allclose(northing, 3348167.264) |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | def test_UTM_4(self): |
---|
| 94 | #Test 4 (Kobenhavn, Northern hemisphere) |
---|
| 95 | lat = 55.70248 |
---|
| 96 | dd,mm,ss = decimal_degrees2degminsec(lat) |
---|
| 97 | |
---|
| 98 | lon = 12.58364 |
---|
| 99 | dd,mm,ss = decimal_degrees2degminsec(lon) |
---|
| 100 | |
---|
| 101 | zone, easting, northing = redfearn(lat,lon) |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | assert zone == 33 |
---|
| 105 | assert num.allclose(easting, 348157.631) |
---|
| 106 | assert num.allclose(northing, 6175612.993) |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | def test_UTM_5(self): |
---|
| 110 | #Test 5 (Wollongong) |
---|
| 111 | |
---|
| 112 | lat = degminsec2decimal_degrees(-34,30,0.) |
---|
| 113 | lon = degminsec2decimal_degrees(150,55,0.) |
---|
| 114 | |
---|
| 115 | zone, easting, northing = redfearn(lat,lon) |
---|
| 116 | |
---|
| 117 | #print zone, easting, northing |
---|
| 118 | |
---|
| 119 | assert zone == 56 |
---|
| 120 | assert num.allclose(easting, 308728.009) |
---|
| 121 | assert num.allclose(northing, 6180432.601) |
---|
| 122 | |
---|
| 123 | def test_UTM_6_nonstandard_projection(self): |
---|
[6400] | 124 | """test_UTM_6_nonstandard_projection |
---|
[6149] | 125 | |
---|
[6400] | 126 | Test that projections can be forced to |
---|
| 127 | use other than native zone. |
---|
| 128 | |
---|
| 129 | Data is from Geraldton, WA |
---|
| 130 | """ |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | # First test native projection (zone 50) |
---|
[6149] | 134 | zone, easting, northing = redfearn(-29.233299999,114.05) |
---|
| 135 | |
---|
| 136 | assert zone == 50 |
---|
| 137 | assert num.allclose(easting, 213251.040253) |
---|
| 138 | assert num.allclose(northing, 6762559.15978) |
---|
| 139 | |
---|
[6400] | 140 | # Testing using the native zone |
---|
[6149] | 141 | zone, easting, northing = redfearn(-29.233299999,114.05, zone=50) |
---|
| 142 | |
---|
| 143 | assert zone == 50 |
---|
| 144 | assert num.allclose(easting, 213251.040253) |
---|
| 145 | assert num.allclose(northing, 6762559.15978) |
---|
| 146 | |
---|
[6400] | 147 | # Then project to zone 49 |
---|
[6149] | 148 | zone, easting, northing = redfearn(-29.233299999,114.05,zone=49) |
---|
| 149 | |
---|
| 150 | assert zone == 49 |
---|
| 151 | assert num.allclose(easting, 796474.020057) |
---|
| 152 | assert num.allclose(northing, 6762310.25162) |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | |
---|
[6400] | 156 | # First test native projection (zone 49) |
---|
[6149] | 157 | zone, easting, northing = redfearn(-29.1333,113.9667) |
---|
| 158 | |
---|
| 159 | assert zone == 49 |
---|
| 160 | assert num.allclose(easting, 788653.192779) |
---|
| 161 | assert num.allclose(northing, 6773605.46384) |
---|
| 162 | |
---|
[6400] | 163 | # Then project to zone 50 |
---|
[6149] | 164 | zone, easting, northing = redfearn(-29.1333,113.9667,zone=50) |
---|
| 165 | |
---|
| 166 | assert zone == 50 |
---|
| 167 | assert num.allclose(easting, 204863.606467) |
---|
| 168 | assert num.allclose(northing, 6773440.04726) |
---|
| 169 | |
---|
[6400] | 170 | # Testing point on zone boundary |
---|
| 171 | # First test native projection (zone 50) |
---|
[6149] | 172 | zone, easting, northing = redfearn(-29.1667,114) |
---|
| 173 | |
---|
| 174 | assert zone == 50 |
---|
| 175 | assert num.allclose(easting, 208199.768268) |
---|
| 176 | assert num.allclose(northing, 6769820.01453) |
---|
| 177 | |
---|
[6400] | 178 | # Then project to zone 49 |
---|
[6149] | 179 | zone, easting, northing = redfearn(-29.1667,114,zone=49) |
---|
| 180 | |
---|
| 181 | assert zone == 49 |
---|
| 182 | assert num.allclose(easting, 791800.231817) |
---|
| 183 | assert num.allclose(northing, 6769820.01453) |
---|
| 184 | |
---|
[6400] | 185 | # Testing furthest point in Geraldton scenario) |
---|
| 186 | # First test native projection (zone 49) |
---|
[6149] | 187 | zone, easting, northing = redfearn(-28.2167,113.4167) |
---|
| 188 | |
---|
| 189 | assert zone == 49 |
---|
| 190 | assert num.allclose(easting, 737178.16131) |
---|
| 191 | assert num.allclose(northing, 6876426.38578) |
---|
| 192 | |
---|
[6400] | 193 | # Then project to zone 50 |
---|
[6149] | 194 | zone, easting, northing = redfearn(-28.2167,113.4167,zone=50) |
---|
| 195 | |
---|
| 196 | assert zone == 50 |
---|
| 197 | assert num.allclose(easting, 148260.567427) |
---|
| 198 | assert num.allclose(northing, 6873587.50926) |
---|
| 199 | |
---|
[6400] | 200 | # Testing outside GDA zone (New Zeland) |
---|
| 201 | # First test native projection (zone 60) |
---|
[6149] | 202 | zone, easting, northing = redfearn(-44,178) |
---|
| 203 | |
---|
| 204 | assert zone == 60 |
---|
| 205 | assert num.allclose(easting, 580174.259843) |
---|
| 206 | assert num.allclose(northing, 5127641.114461) |
---|
| 207 | |
---|
[6400] | 208 | # Then project to zone 59 |
---|
[6149] | 209 | zone, easting, northing = redfearn(-44,178,zone=59) |
---|
| 210 | |
---|
| 211 | assert zone == 59 |
---|
| 212 | assert num.allclose(easting, 1061266.922118) |
---|
| 213 | assert num.allclose(northing, 5104249.395469) |
---|
| 214 | |
---|
[6400] | 215 | # Then skip three zones 57 (native 60) |
---|
[6149] | 216 | zone, easting, northing = redfearn(-44,178,zone=57) |
---|
| 217 | |
---|
| 218 | assert zone == 57 |
---|
| 219 | assert num.allclose(easting, 2023865.527497) |
---|
| 220 | assert num.allclose(northing, 4949253.934967) |
---|
| 221 | |
---|
| 222 | # Note Projecting into the Northern Hemisphere Does not coincide |
---|
| 223 | # redfearn or ArcMap conversions. The difference lies in |
---|
| 224 | # False Northing which is 0 for the Northern hemisphere |
---|
| 225 | # in the redfearn implementation but 10 million in Arc. |
---|
| 226 | # |
---|
| 227 | # But the redfearn implementation does coincide with |
---|
| 228 | # Google Earth (he he) |
---|
| 229 | |
---|
[6400] | 230 | # Testing outside GDA zone (Northern Hemisphere) |
---|
| 231 | # First test native projection (zone 57) |
---|
[6149] | 232 | zone, easting, northing = redfearn(44,156) |
---|
| 233 | |
---|
| 234 | # Google Earth interpretation |
---|
| 235 | assert zone == 57 |
---|
| 236 | assert num.allclose(easting, 259473.69) |
---|
| 237 | assert num.allclose(northing, 4876249.13) |
---|
| 238 | |
---|
| 239 | # ArcMap's interpretation |
---|
| 240 | #assert zone == 57 |
---|
| 241 | #assert num.allclose(easting, 259473.678944) |
---|
| 242 | #assert num.allclose(northing, 14876249.1268) |
---|
| 243 | |
---|
[6400] | 244 | # Then project to zone 56 |
---|
[6149] | 245 | zone, easting, northing = redfearn(44,156,zone=56) |
---|
| 246 | |
---|
| 247 | assert zone == 56 |
---|
| 248 | assert num.allclose(easting, 740526.321055) |
---|
| 249 | assert num.allclose(northing, 4876249.13) |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | #def test_UTM_6(self): |
---|
| 255 | # """Test 6 (Don's Wollongong file's ref point) |
---|
| 256 | # """ |
---|
| 257 | # |
---|
| 258 | # lat = -34.490286785873 |
---|
| 259 | # lon = 150.79712139578 |
---|
| 260 | # |
---|
| 261 | # dd,mm,ss = decimal_degrees2degminsec(lat) |
---|
| 262 | # print dd,mm,ss |
---|
| 263 | # dd,mm,ss = decimal_degrees2degminsec(lon) |
---|
| 264 | # print dd,mm,ss |
---|
| 265 | # |
---|
| 266 | # zone, easting, northing = redfearn(lat,lon) |
---|
| 267 | # |
---|
| 268 | # print zone, easting, northing |
---|
| 269 | # |
---|
| 270 | # assert zone == 56 |
---|
| 271 | # #assert allclose(easting, 297717.36468927) #out by 10m |
---|
| 272 | # #assert allclose(northing, 6181725.1724276) |
---|
| 273 | |
---|
[6400] | 274 | |
---|
| 275 | def test_nonstandard_meridian(self): |
---|
| 276 | """test_nonstandard_meridian |
---|
| 277 | |
---|
| 278 | This test will verify that redfearn can be used to project |
---|
| 279 | points using an arbitrary central meridian. |
---|
| 280 | """ |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | # FIXME: To do using csv file |
---|
| 284 | pass |
---|
| 285 | |
---|
| 286 | |
---|
| 287 | # The file projection_test_points.csv contains 10 points |
---|
| 288 | # which straddle the boundary between UTM zones 53 and 54. |
---|
| 289 | # They have been projected using a central meridian of 137.5 |
---|
| 290 | # degrees (the boundary is 138 so it is pretty much right |
---|
| 291 | # in the middle of zones 53 and 54). |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | |
---|
| 296 | |
---|
[6149] | 297 | def test_convert_lats_longs(self): |
---|
| 298 | |
---|
| 299 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 300 | #Zone: 56 |
---|
| 301 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 302 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 303 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 304 | |
---|
| 305 | lat_gong = degminsec2decimal_degrees(-34,30,0.) |
---|
| 306 | lon_gong = degminsec2decimal_degrees(150,55,0.) |
---|
| 307 | |
---|
| 308 | lat_2 = degminsec2decimal_degrees(-34,00,0.) |
---|
| 309 | lon_2 = degminsec2decimal_degrees(150,00,0.) |
---|
| 310 | |
---|
| 311 | lats = [lat_gong, lat_2] |
---|
| 312 | longs = [lon_gong, lon_2] |
---|
| 313 | points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs) |
---|
| 314 | |
---|
| 315 | assert num.allclose(points[0][0], 308728.009) |
---|
| 316 | assert num.allclose(points[0][1], 6180432.601) |
---|
| 317 | assert num.allclose(points[1][0], 222908.705) |
---|
| 318 | assert num.allclose(points[1][1], 6233785.284) |
---|
| 319 | self.failUnless(zone == 56, |
---|
| 320 | 'Bad zone error!') |
---|
| 321 | |
---|
| 322 | def test_convert_lats_longs2(self): |
---|
| 323 | |
---|
| 324 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 325 | #Zone: 56 |
---|
| 326 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 327 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 328 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 329 | |
---|
| 330 | lat_gong = degminsec2decimal_degrees(-34,30,0.) |
---|
| 331 | lon_gong = degminsec2decimal_degrees(150,55,0.) |
---|
| 332 | |
---|
| 333 | lat_2 = degminsec2decimal_degrees(34,00,0.) |
---|
| 334 | lon_2 = degminsec2decimal_degrees(100,00,0.) |
---|
| 335 | |
---|
| 336 | lats = [lat_gong, lat_2] |
---|
| 337 | longs = [lon_gong, lon_2] |
---|
| 338 | |
---|
| 339 | try: |
---|
| 340 | points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs) |
---|
| 341 | except ANUGAError: |
---|
| 342 | pass |
---|
| 343 | else: |
---|
| 344 | self.failUnless(False, |
---|
| 345 | 'Error not thrown error!') |
---|
| 346 | |
---|
| 347 | def test_convert_lats_longs3(self): |
---|
| 348 | |
---|
| 349 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 350 | #Zone: 56 |
---|
| 351 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 352 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 353 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 354 | |
---|
| 355 | lat_gong = "-34.5" |
---|
| 356 | lon_gong = "150.916666667" |
---|
| 357 | lat_2 = degminsec2decimal_degrees(34,00,0.) |
---|
| 358 | lon_2 = degminsec2decimal_degrees(100,00,0.) |
---|
| 359 | |
---|
| 360 | lats = [lat_gong, lat_2] |
---|
| 361 | longs = [lon_gong, lon_2] |
---|
| 362 | try: |
---|
| 363 | points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs) |
---|
| 364 | except ANUGAError: |
---|
| 365 | pass |
---|
| 366 | else: |
---|
| 367 | self.failUnless(False, |
---|
| 368 | 'Error not thrown error!') |
---|
| 369 | |
---|
| 370 | def test_convert_latlon_to_UTM1(self): |
---|
| 371 | |
---|
| 372 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 373 | #Zone: 56 |
---|
| 374 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 375 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 376 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 377 | |
---|
| 378 | lat_gong = degminsec2decimal_degrees(-34,30,0.) |
---|
| 379 | lon_gong = degminsec2decimal_degrees(150,55,0.) |
---|
| 380 | |
---|
| 381 | lat_2 = degminsec2decimal_degrees(-34,00,0.) |
---|
| 382 | lon_2 = degminsec2decimal_degrees(150,00,0.) |
---|
| 383 | |
---|
| 384 | points = [[lat_gong, lon_gong], [lat_2, lon_2]] |
---|
| 385 | points, zone = convert_from_latlon_to_utm(points=points) |
---|
| 386 | #print "points",points |
---|
| 387 | assert num.allclose(points[0][0], 308728.009) |
---|
| 388 | assert num.allclose(points[0][1], 6180432.601) |
---|
| 389 | assert num.allclose(points[1][0], 222908.705) |
---|
| 390 | assert num.allclose(points[1][1], 6233785.284) |
---|
| 391 | self.failUnless(zone == 56, |
---|
| 392 | 'Bad zone error!') |
---|
| 393 | |
---|
| 394 | def test_convert_latlon_to_UTM2(self): |
---|
| 395 | |
---|
| 396 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 397 | #Zone: 56 |
---|
| 398 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 399 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 400 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 401 | |
---|
| 402 | lat_gong = degminsec2decimal_degrees(-34,30,0.) |
---|
| 403 | lon_gong = degminsec2decimal_degrees(150,55,0.) |
---|
| 404 | |
---|
| 405 | lat_2 = degminsec2decimal_degrees(34,00,0.) |
---|
| 406 | lon_2 = degminsec2decimal_degrees(100,00,0.) |
---|
| 407 | |
---|
| 408 | points = [[lat_gong, lon_gong], [lat_2, lon_2]] |
---|
| 409 | |
---|
| 410 | try: |
---|
| 411 | points, zone = convert_from_latlon_to_utm(points=points) |
---|
| 412 | except ANUGAError: |
---|
| 413 | pass |
---|
| 414 | else: |
---|
| 415 | self.fail('Error not thrown error!') |
---|
| 416 | |
---|
| 417 | def test_convert_latlon_to_UTM3(self): |
---|
| 418 | |
---|
| 419 | #Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) |
---|
| 420 | #Zone: 56 |
---|
| 421 | #Easting: 222908.705 Northing: 6233785.284 |
---|
| 422 | #Latitude: -34 0 ' 0.00000 '' Longitude: 150 0 ' 0.00000 '' |
---|
| 423 | #Grid Convergence: -1 40 ' 43.13 '' Point Scale: 1.00054660 |
---|
| 424 | |
---|
| 425 | lat_gong = "-34.5" |
---|
| 426 | lon_gong = "150.916666667" |
---|
| 427 | lat_2 = degminsec2decimal_degrees(34,00,0.) |
---|
| 428 | lon_2 = degminsec2decimal_degrees(100,00,0.) |
---|
| 429 | |
---|
| 430 | points = [[lat_gong, lon_gong], [lat_2, lon_2]] |
---|
| 431 | |
---|
| 432 | try: |
---|
| 433 | points, zone = convert_from_latlon_to_utm(points=points) |
---|
| 434 | except ANUGAError: |
---|
| 435 | pass |
---|
| 436 | else: |
---|
| 437 | self.fail('Error not thrown error!') |
---|
| 438 | |
---|
| 439 | |
---|
| 440 | #------------------------------------------------------------- |
---|
| 441 | if __name__ == "__main__": |
---|
| 442 | |
---|
| 443 | #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1') |
---|
| 444 | #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection') |
---|
| 445 | mysuite = unittest.makeSuite(TestCase,'test') |
---|
| 446 | runner = unittest.TextTestRunner() |
---|
| 447 | runner.run(mysuite) |
---|