source: trunk/anuga_core/source/anuga/coordinate_transforms/test_redfearn.py @ 7967

Last change on this file since 7967 was 7779, checked in by hudson, 15 years ago

Fixed wrong exception module import paths.
File size: 16.8 KB
Line 
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
8import unittest
9
10from redfearn import *
11from anuga.anuga_exceptions import ANUGAError
12
13from anuga.utilities.system_tools import get_pathname_from_package
14from os.path import join
15import numpy as num
16
17
18#-------------------------------------------------------------
19
20class TestCase(unittest.TestCase):
21
22    def test_decimal_degrees_conversion(Self):
23        lat = degminsec2decimal_degrees(-37,39,10.15610)
24        lon = degminsec2decimal_degrees(143,55,35.38390) 
25        assert num.allclose(lat, -37.65282114)
26        assert num.allclose(lon, 143.9264955)
27
28        dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
29        assert dd==-37
30        assert mm==39
31        assert num.allclose(ss, 10.15610)
32
33        dd,mm,ss = decimal_degrees2degminsec(143.9264955)
34        assert dd==143
35        assert mm==55
36        assert num.allclose(ss, 35.38390) 
37
38
39    def test_UTM_1(self):
40        #latitude:  -37 39' 10.15610"
41        #Longitude: 143 55' 35.38390"
42        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
43        #Zone:   54   
44        #Easting:  758173.797  Northing: 5828674.340
45        #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 ''
46        #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107
47
48        lat = degminsec2decimal_degrees(-37,39,10.15610)
49        lon = degminsec2decimal_degrees(143,55,35.38390) 
50        assert num.allclose(lat, -37.65282114)
51        assert num.allclose(lon, 143.9264955)
52
53
54        zone, easting, northing = redfearn(lat,lon)
55
56        assert zone == 54
57        assert num.allclose(easting, 758173.797)
58        assert num.allclose(northing, 5828674.340)
59
60
61    def test_UTM_2(self):
62        #TEST 2
63
64        #Latitude:  -37 57 03.7203
65        #Longitude: 144 25 29.5244
66        #Zone:   55   
67        #Easting:  273741.297  Northing: 5796489.777
68        #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 ''
69        #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056
70
71        lat = degminsec2decimal_degrees(-37,57,03.7203)
72        lon = degminsec2decimal_degrees(144,25,29.5244) 
73
74        zone, easting, northing = redfearn(lat,lon)
75
76        assert zone == 55
77        assert num.allclose(easting, 273741.297)
78        assert num.allclose(northing, 5796489.777)
79
80       
81    def test_UTM_3(self):
82        #Test 3
83        lat = degminsec2decimal_degrees(-60,0,0)
84        lon = degminsec2decimal_degrees(130,0,0) 
85
86        zone, easting, northing = redfearn(lat,lon)
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        assert zone == 56
118        assert num.allclose(easting, 308728.009)
119        assert num.allclose(northing, 6180432.601)
120
121    def test_UTM_6_nonstandard_projection(self):
122        """test_UTM_6_nonstandard_projection
123
124        Test that projections can be forced to
125        use other than native zone.
126
127        Data is from Geraldton, WA
128        """
129       
130
131        # First test native projection (zone 50)
132        zone, easting, northing = redfearn(-29.233299999,114.05)
133
134        assert zone == 50
135        assert num.allclose(easting, 213251.040253)
136        assert num.allclose(northing, 6762559.15978)
137
138        #Testing using the native zone
139        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
140
141        assert zone == 50
142        assert num.allclose(easting, 213251.040253)
143        assert num.allclose(northing, 6762559.15978)
144
145        #Then project to zone 49
146        zone, easting, northing = redfearn(-29.233299999,114.05,zone=49)
147
148        assert zone == 49
149        assert num.allclose(easting, 796474.020057)
150        assert num.allclose(northing, 6762310.25162)
151
152       
153
154       
155
156        #First test native projection (zone 49)
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
163        #Then project to zone 50
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
170        #Testing point on zone boundary
171        #First test native projection (zone 50)
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
178        #Then project to zone 49
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
185        #Testing furthest point in Geraldton scenario)
186        #First test native projection (zone 49)
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
193        #Then project to zone 50
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
200        #Testing outside GDA zone (New Zeland)
201        #First test native projection (zone 60)
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
208        #Then project to zone 59
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
215        #Then skip three zones 57 (native 60)
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       
230        #Testing outside GDA zone (Northern Hemisphere)
231        #First test native projection (zone 57)
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       
244        #Then project to zone 56
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
274   
275    def test_nonstandard_meridian_coinciding_with_native(self):
276        """test_nonstandard_meridian_coinciding_with_native
277
278        This test will verify that redfearn can be used to project
279        points using an arbitrary central meridian that happens to
280        coincide with the standard meridian at the center of a UTM zone.
281        This is a preliminary test before testing this functionality
282        with a truly arbitrary non-standard meridian.
283        """
284
285        # The file projection_test_points_z53.csv contains 10 points
286        # which straddle the boundary between UTM zones 53 and 54.
287        # They have been projected to zone 53 irrespective of where they
288        # belong.
289
290        path = get_pathname_from_package('anuga.coordinate_transforms')
291       
292        for forced_zone in [53, 54]:
293       
294            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
295            fid = open(datafile)
296
297            for line in fid.readlines()[1:]:
298                fields = line.strip().split(',')
299               
300                lon = float(fields[1])
301                lat = float(fields[2])
302                x = float(fields[3])
303                y = float(fields[4])           
304
305                zone, easting, northing = redfearn(lat, lon,
306                                                   zone=forced_zone)
307               
308                # Check calculation
309                assert zone == forced_zone
310                assert num.allclose(x, easting)
311                assert num.allclose(y, northing)
312
313   
314   
315   
316    def test_nonstandard_meridian(self):
317        """test_nonstandard_meridian
318
319        This test will verify that redfearn can be used to project
320        points using an arbitrary central meridian.
321        """
322
323        # The file projection_test_points.csv contains 10 points
324        # which straddle the boundary between UTM zones 53 and 54.
325        # They have been projected using a central meridian of 137.5
326        # degrees (the boundary is 138 so it is pretty much right
327        # in the middle of zones 53 and 54).
328
329        path = get_pathname_from_package('anuga.coordinate_transforms')
330        datafile = join(path, 'projection_test_points.csv')
331        fid = open(datafile)
332
333        for line in fid.readlines()[1:]:
334            fields = line.strip().split(',')
335
336            lon = float(fields[1])
337            lat = float(fields[2])
338            x = float(fields[3])
339            y = float(fields[4])           
340
341            zone, easting, northing = redfearn(lat, lon,
342                                               central_meridian=137.5,
343                                               scale_factor=0.9996)
344
345            assert zone == -1 # Indicates non UTM projection
346            assert num.allclose(x, easting)
347            assert num.allclose(y, northing)
348
349        # Test that zone and meridian can't both be specified
350        try:
351            zone, easting, northing = redfearn(lat, lon,
352                                               zone=50, 
353                                               central_meridian=137.5)
354        except:
355            pass
356        else:
357            msg = 'Should have raised exception'
358            raise Exception, msg
359
360           
361    def test_convert_lats_longs(self):
362
363        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
364        #Zone:   56   
365        #Easting:  222908.705  Northing: 6233785.284
366        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
367        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
368
369        lat_gong = degminsec2decimal_degrees(-34,30,0.)
370        lon_gong = degminsec2decimal_degrees(150,55,0.)
371       
372        lat_2 = degminsec2decimal_degrees(-34,00,0.)
373        lon_2 = degminsec2decimal_degrees(150,00,0.)
374       
375        lats = [lat_gong, lat_2]
376        longs = [lon_gong, lon_2]
377        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
378
379        assert num.allclose(points[0][0], 308728.009)
380        assert num.allclose(points[0][1], 6180432.601)
381        assert num.allclose(points[1][0],  222908.705)
382        assert num.allclose(points[1][1], 6233785.284)
383        self.failUnless(zone == 56,
384                        'Bad zone error!')
385       
386    def test_convert_lats_longs2(self):
387
388        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
389        #Zone:   56   
390        #Easting:  222908.705  Northing: 6233785.284
391        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
392        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
393
394        lat_gong = degminsec2decimal_degrees(-34,30,0.)
395        lon_gong = degminsec2decimal_degrees(150,55,0.)
396       
397        lat_2 = degminsec2decimal_degrees(34,00,0.)
398        lon_2 = degminsec2decimal_degrees(100,00,0.)
399       
400        lats = [lat_gong, lat_2]
401        longs = [lon_gong, lon_2]
402       
403        try:
404            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
405        except ANUGAError:
406            pass
407        else:
408            self.failUnless(False,
409                            'Error not thrown error!')
410           
411    def test_convert_lats_longs3(self):
412
413        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
414        #Zone:   56   
415        #Easting:  222908.705  Northing: 6233785.284
416        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
417        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
418
419        lat_gong = "-34.5"
420        lon_gong = "150.916666667"
421        lat_2 = degminsec2decimal_degrees(34,00,0.)
422        lon_2 = degminsec2decimal_degrees(100,00,0.)
423       
424        lats = [lat_gong, lat_2]
425        longs = [lon_gong, lon_2]
426        try:
427            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
428        except ANUGAError:
429            pass
430        else:
431            self.failUnless(False,
432                            'Error not thrown error!')
433
434    def test_convert_latlon_to_UTM1(self):
435
436        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
437        #Zone:   56   
438        #Easting:  222908.705  Northing: 6233785.284
439        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
440        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
441
442        lat_gong = degminsec2decimal_degrees(-34,30,0.)
443        lon_gong = degminsec2decimal_degrees(150,55,0.)
444       
445        lat_2 = degminsec2decimal_degrees(-34,00,0.)
446        lon_2 = degminsec2decimal_degrees(150,00,0.)
447       
448        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
449        points, zone = convert_from_latlon_to_utm(points=points)
450        assert num.allclose(points[0][0], 308728.009)
451        assert num.allclose(points[0][1], 6180432.601)
452        assert num.allclose(points[1][0],  222908.705)
453        assert num.allclose(points[1][1], 6233785.284)
454        self.failUnless(zone == 56,
455                        'Bad zone error!')
456
457    def test_convert_latlon_to_UTM2(self):       
458
459        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
460        #Zone:   56   
461        #Easting:  222908.705  Northing: 6233785.284
462        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
463        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
464
465        lat_gong = degminsec2decimal_degrees(-34,30,0.)
466        lon_gong = degminsec2decimal_degrees(150,55,0.)
467       
468        lat_2 = degminsec2decimal_degrees(34,00,0.)
469        lon_2 = degminsec2decimal_degrees(100,00,0.)
470
471        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
472
473        try:
474            points, zone = convert_from_latlon_to_utm(points=points)           
475        except ANUGAError:
476            pass
477        else:
478            self.fail('Error not thrown error!')
479
480    def test_convert_latlon_to_UTM3(self):           
481
482        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
483        #Zone:   56   
484        #Easting:  222908.705  Northing: 6233785.284
485        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
486        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
487
488        lat_gong = "-34.5"
489        lon_gong = "150.916666667"
490        lat_2 = degminsec2decimal_degrees(34,00,0.)
491        lon_2 = degminsec2decimal_degrees(100,00,0.)
492
493        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
494
495        try:
496            points, zone = convert_from_latlon_to_utm(points=points)           
497        except ANUGAError:
498            pass
499        else:
500            self.fail('Error not thrown error!')
501
502           
503#-------------------------------------------------------------
504
505if __name__ == "__main__":
506    mysuite = unittest.makeSuite(TestCase,'test')
507    runner = unittest.TextTestRunner()
508    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.