source: anuga_core/source/anuga/coordinate_transforms/test_redfearn.py @ 6595

Last change on this file since 6595 was 6595, checked in by ole, 16 years ago

Added the nonstandard central meridian option, tested and implemented it in the Adeladie Phase II scripts.
File size: 17.0 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.utilities.anuga_exceptions import ANUGAError
12from anuga.utilities.system_tools import get_pathname_from_package
13from os.path import join
14import Numeric as num
15
16
17#-------------------------------------------------------------
18
19class TestCase(unittest.TestCase):
20
21    def test_decimal_degrees_conversion(Self):
22        lat = degminsec2decimal_degrees(-37,39,10.15610)
23        lon = degminsec2decimal_degrees(143,55,35.38390) 
24        assert num.allclose(lat, -37.65282114)
25        assert num.allclose(lon, 143.9264955)
26
27        dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
28        assert dd==-37
29        assert mm==39
30        assert num.allclose(ss, 10.15610)
31
32        dd,mm,ss = decimal_degrees2degminsec(143.9264955)
33        assert dd==143
34        assert mm==55
35        assert num.allclose(ss, 35.38390) 
36
37
38    def test_UTM_1(self):
39        #latitude:  -37 39' 10.15610"
40        #Longitude: 143 55' 35.38390"
41        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
42        #Zone:   54   
43        #Easting:  758173.797  Northing: 5828674.340
44        #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 ''
45        #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107
46
47        lat = degminsec2decimal_degrees(-37,39,10.15610)
48        lon = degminsec2decimal_degrees(143,55,35.38390) 
49        assert num.allclose(lat, -37.65282114)
50        assert num.allclose(lon, 143.9264955)
51
52
53        zone, easting, northing = redfearn(lat,lon)
54
55        assert zone == 54
56        assert num.allclose(easting, 758173.797)
57        assert num.allclose(northing, 5828674.340)
58
59
60    def test_UTM_2(self):
61        #TEST 2
62
63        #Latitude:  -37 57 03.7203
64        #Longitude: 144 25 29.5244
65        #Zone:   55   
66        #Easting:  273741.297  Northing: 5796489.777
67        #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 ''
68        #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056
69
70        lat = degminsec2decimal_degrees(-37,57,03.7203)
71        lon = degminsec2decimal_degrees(144,25,29.5244) 
72        #print lat, lon
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        #print zone, easting, northing
88
89        assert zone == 52
90        assert num.allclose(easting, 555776.267)
91        assert num.allclose(northing, 3348167.264)
92
93
94    def test_UTM_4(self):
95        #Test 4 (Kobenhavn, Northern hemisphere)
96        lat = 55.70248
97        dd,mm,ss = decimal_degrees2degminsec(lat)
98
99        lon = 12.58364
100        dd,mm,ss = decimal_degrees2degminsec(lon)
101
102        zone, easting, northing = redfearn(lat,lon)
103
104
105        assert zone == 33
106        assert num.allclose(easting, 348157.631)
107        assert num.allclose(northing, 6175612.993) 
108
109
110    def test_UTM_5(self):
111        #Test 5 (Wollongong)
112
113        lat = degminsec2decimal_degrees(-34,30,0.)
114        lon = degminsec2decimal_degrees(150,55,0.) 
115       
116        zone, easting, northing = redfearn(lat,lon)
117
118        #print zone, easting, northing
119
120        assert zone == 56
121        assert num.allclose(easting, 308728.009)
122        assert num.allclose(northing, 6180432.601)
123
124    def test_UTM_6_nonstandard_projection(self):
125        """test_UTM_6_nonstandard_projection
126
127        Test that projections can be forced to
128        use other than native zone.
129
130        Data is from Geraldton, WA
131        """
132       
133
134        # First test native projection (zone 50)
135        zone, easting, northing = redfearn(-29.233299999,114.05)
136
137        assert zone == 50
138        assert num.allclose(easting, 213251.040253)
139        assert num.allclose(northing, 6762559.15978)
140
141        # Testing using the native zone
142        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
143
144        assert zone == 50
145        assert num.allclose(easting, 213251.040253)
146        assert num.allclose(northing, 6762559.15978)
147
148        # Then project to zone 49
149        zone, easting, northing = redfearn(-29.233299999,114.05, zone=49)
150
151        assert zone == 49
152        assert num.allclose(easting, 796474.020057)
153        assert num.allclose(northing, 6762310.25162)
154
155       
156
157        # First test native projection (zone 49)
158        zone, easting, northing = redfearn(-29.1333,113.9667)
159
160        assert zone == 49
161        assert num.allclose(easting, 788653.192779)
162        assert num.allclose(northing, 6773605.46384)
163
164        # Then project to zone 50
165        zone, easting, northing = redfearn(-29.1333,113.9667, zone=50)
166
167        assert zone == 50
168        assert num.allclose(easting, 204863.606467)
169        assert num.allclose(northing, 6773440.04726)
170
171        # Testing point on zone boundary
172        # First test native projection (zone 50)
173        zone, easting, northing = redfearn(-29.1667,114)
174
175        assert zone == 50
176        assert num.allclose(easting, 208199.768268)
177        assert num.allclose(northing, 6769820.01453)
178
179        # Then project to zone 49
180        zone, easting, northing = redfearn(-29.1667,114, zone=49)
181
182        assert zone == 49
183        assert num.allclose(easting, 791800.231817)
184        assert num.allclose(northing, 6769820.01453)
185
186        # Testing furthest point in Geraldton scenario)
187        # First test native projection (zone 49)
188        zone, easting, northing = redfearn(-28.2167,113.4167)
189
190        assert zone == 49
191        assert num.allclose(easting, 737178.16131)
192        assert num.allclose(northing, 6876426.38578)
193
194        # Then project to zone 50
195        zone, easting, northing = redfearn(-28.2167,113.4167,zone=50)
196
197        assert zone == 50
198        assert num.allclose(easting, 148260.567427)
199        assert num.allclose(northing, 6873587.50926)
200
201        # Testing outside GDA zone (New Zeland)
202        # First test native projection (zone 60)
203        zone, easting, northing = redfearn(-44,178)
204
205        assert zone == 60
206        assert num.allclose(easting, 580174.259843)
207        assert num.allclose(northing, 5127641.114461)
208
209        # Then project to zone 59
210        zone, easting, northing = redfearn(-44,178,zone=59)
211
212        assert zone == 59
213        assert num.allclose(easting, 1061266.922118)
214        assert num.allclose(northing, 5104249.395469)
215
216        # Then skip three zones 57 (native 60)
217        zone, easting, northing = redfearn(-44,178,zone=57)
218
219        assert zone == 57
220        assert num.allclose(easting, 2023865.527497)
221        assert num.allclose(northing, 4949253.934967)
222
223        # Note Projecting into the Northern Hemisphere Does not coincide
224        # redfearn or ArcMap conversions. The difference lies in
225        # False Northing which is 0 for the Northern hemisphere
226        # in the redfearn implementation but 10 million in Arc.
227        #
228        # But the redfearn implementation does coincide with
229        # Google Earth (he he)
230       
231        # Testing outside GDA zone (Northern Hemisphere)
232        # First test native projection (zone 57)
233        zone, easting, northing = redfearn(44,156)
234
235        # Google Earth interpretation
236        assert zone == 57
237        assert num.allclose(easting, 259473.69)
238        assert num.allclose(northing, 4876249.13)
239
240        # ArcMap's interpretation
241        #assert zone == 57       
242        #assert num.allclose(easting, 259473.678944)
243        #assert num.allclose(northing, 14876249.1268)
244       
245        # Then project to zone 56
246        zone, easting, northing = redfearn(44,156,zone=56)
247
248        assert zone == 56
249        assert num.allclose(easting, 740526.321055)
250        assert num.allclose(northing, 4876249.13)
251
252       
253
254
255    #def test_UTM_6(self):
256    #    """Test 6 (Don's Wollongong file's ref point)
257    #    """
258    #
259    #    lat = -34.490286785873
260    #    lon = 150.79712139578
261    #
262    #    dd,mm,ss = decimal_degrees2degminsec(lat)
263    #    print dd,mm,ss
264    #    dd,mm,ss = decimal_degrees2degminsec(lon)       
265    #    print dd,mm,ss
266    #     
267    #    zone, easting, northing = redfearn(lat,lon)
268    #
269    #    print zone, easting, northing
270    #
271    #    assert zone == 56
272    #    #assert allclose(easting, 297717.36468927) #out by 10m
273    #    #assert allclose(northing, 6181725.1724276)
274
275
276   
277    def test_nonstandard_meridian_coinciding_with_native(self):
278        """test_nonstandard_meridian_coinciding_with_native
279
280        This test will verify that redfearn can be used to project
281        points using an arbitrary central meridian that happens to
282        coincide with the standard meridian at the center of a UTM zone.
283        This is a preliminary test before testing this functionality
284        with a truly arbitrary non-standard meridian.
285        """
286
287        # The file projection_test_points_z53.csv contains 10 points
288        # which straddle the boundary between UTM zones 53 and 54.
289        # They have been projected to zone 53 irrespective of where they
290        # belong.
291
292        path = get_pathname_from_package('anuga.coordinate_transforms')
293       
294        for forced_zone in [53, 54]:
295       
296            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
297            fid = open(datafile)
298
299            for line in fid.readlines()[1:]:
300                fields = line.strip().split(',')
301               
302                lon = float(fields[1])
303                lat = float(fields[2])
304                x = float(fields[3])
305                y = float(fields[4])           
306
307                zone, easting, northing = redfearn(lat, lon,
308                                                   zone=forced_zone)
309               
310                # Check calculation
311                assert zone == forced_zone
312                assert num.allclose(x, easting)
313                assert num.allclose(y, northing)
314
315   
316   
317   
318    def test_nonstandard_meridian(self):
319        """test_nonstandard_meridian
320
321        This test will verify that redfearn can be used to project
322        points using an arbitrary central meridian.
323        """
324
325        # The file projection_test_points.csv contains 10 points
326        # which straddle the boundary between UTM zones 53 and 54.
327        # They have been projected using a central meridian of 137.5
328        # degrees (the boundary is 138 so it is pretty much right
329        # in the middle of zones 53 and 54).
330
331        path = get_pathname_from_package('anuga.coordinate_transforms')
332        datafile = join(path, 'projection_test_points.csv')
333        fid = open(datafile)
334
335        for line in fid.readlines()[1:]:
336            fields = line.strip().split(',')
337
338            lon = float(fields[1])
339            lat = float(fields[2])
340            x = float(fields[3])
341            y = float(fields[4])           
342
343            zone, easting, northing = redfearn(lat, lon,
344                                               central_meridian=137.5,
345                                               scale_factor=0.9996)
346
347            assert zone == -1 # Indicates non UTM projection
348            assert num.allclose(x, easting)
349            assert num.allclose(y, northing)
350
351        # Test that zone and meridian can't both be specified
352        try:
353            zone, easting, northing = redfearn(lat, lon,
354                                               zone=50, 
355                                               central_meridian=137.5)
356        except:
357            pass
358        else:
359            msg = 'Should have raised exception'
360            raise Exception, msg
361
362           
363    def test_convert_lats_longs(self):
364
365        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
366        #Zone:   56   
367        #Easting:  222908.705  Northing: 6233785.284
368        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
369        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
370
371        lat_gong = degminsec2decimal_degrees(-34,30,0.)
372        lon_gong = degminsec2decimal_degrees(150,55,0.)
373       
374        lat_2 = degminsec2decimal_degrees(-34,00,0.)
375        lon_2 = degminsec2decimal_degrees(150,00,0.)
376       
377        lats = [lat_gong, lat_2]
378        longs = [lon_gong, lon_2]
379        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
380
381        assert num.allclose(points[0][0], 308728.009)
382        assert num.allclose(points[0][1], 6180432.601)
383        assert num.allclose(points[1][0],  222908.705)
384        assert num.allclose(points[1][1], 6233785.284)
385        self.failUnless(zone == 56,
386                        'Bad zone error!')
387       
388    def test_convert_lats_longs2(self):
389
390        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
391        #Zone:   56   
392        #Easting:  222908.705  Northing: 6233785.284
393        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
394        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
395
396        lat_gong = degminsec2decimal_degrees(-34,30,0.)
397        lon_gong = degminsec2decimal_degrees(150,55,0.)
398       
399        lat_2 = degminsec2decimal_degrees(34,00,0.)
400        lon_2 = degminsec2decimal_degrees(100,00,0.)
401       
402        lats = [lat_gong, lat_2]
403        longs = [lon_gong, lon_2]
404       
405        try:
406            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
407        except ANUGAError:
408            pass
409        else:
410            self.failUnless(False,
411                            'Error not thrown error!')
412           
413    def test_convert_lats_longs3(self):
414
415        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
416        #Zone:   56   
417        #Easting:  222908.705  Northing: 6233785.284
418        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
419        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
420
421        lat_gong = "-34.5"
422        lon_gong = "150.916666667"
423        lat_2 = degminsec2decimal_degrees(34,00,0.)
424        lon_2 = degminsec2decimal_degrees(100,00,0.)
425       
426        lats = [lat_gong, lat_2]
427        longs = [lon_gong, lon_2]
428        try:
429            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
430        except ANUGAError:
431            pass
432        else:
433            self.failUnless(False,
434                            'Error not thrown error!')
435
436    def test_convert_latlon_to_UTM1(self):
437
438        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
439        #Zone:   56   
440        #Easting:  222908.705  Northing: 6233785.284
441        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
442        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
443
444        lat_gong = degminsec2decimal_degrees(-34,30,0.)
445        lon_gong = degminsec2decimal_degrees(150,55,0.)
446       
447        lat_2 = degminsec2decimal_degrees(-34,00,0.)
448        lon_2 = degminsec2decimal_degrees(150,00,0.)
449       
450        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
451        points, zone = convert_from_latlon_to_utm(points=points)
452        #print "points",points
453        assert num.allclose(points[0][0], 308728.009)
454        assert num.allclose(points[0][1], 6180432.601)
455        assert num.allclose(points[1][0],  222908.705)
456        assert num.allclose(points[1][1], 6233785.284)
457        self.failUnless(zone == 56,
458                        'Bad zone error!')
459
460    def test_convert_latlon_to_UTM2(self):       
461
462        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
463        #Zone:   56   
464        #Easting:  222908.705  Northing: 6233785.284
465        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
466        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
467
468        lat_gong = degminsec2decimal_degrees(-34,30,0.)
469        lon_gong = degminsec2decimal_degrees(150,55,0.)
470       
471        lat_2 = degminsec2decimal_degrees(34,00,0.)
472        lon_2 = degminsec2decimal_degrees(100,00,0.)
473
474        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
475
476        try:
477            points, zone = convert_from_latlon_to_utm(points=points)           
478        except ANUGAError:
479            pass
480        else:
481            self.fail('Error not thrown error!')
482
483    def test_convert_latlon_to_UTM3(self):           
484
485        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
486        #Zone:   56   
487        #Easting:  222908.705  Northing: 6233785.284
488        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
489        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
490
491        lat_gong = "-34.5"
492        lon_gong = "150.916666667"
493        lat_2 = degminsec2decimal_degrees(34,00,0.)
494        lon_2 = degminsec2decimal_degrees(100,00,0.)
495
496        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
497
498        try:
499            points, zone = convert_from_latlon_to_utm(points=points)           
500        except ANUGAError:
501            pass
502        else:
503            self.fail('Error not thrown error!')
504
505           
506#-------------------------------------------------------------
507if __name__ == "__main__":
508
509    mysuite = unittest.makeSuite(TestCase,'test')
510    runner = unittest.TextTestRunner()
511    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.