source: branches/numpy/anuga/coordinate_transforms/test_redfearn.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 16 years ago

Initial commit of numpy changes. Still a long way to go.
File size: 13.7 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
12
13import numpy as num
14
15
16#-------------------------------------------------------------
17
18class 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):
124        #Test 6 (Geraldton, WA)
125
126        #First test native projection (zone 50)
127        zone, easting, northing = redfearn(-29.233299999,114.05)
128
129        assert zone == 50
130        assert num.allclose(easting, 213251.040253)
131        assert num.allclose(northing, 6762559.15978)
132
133        #Testing using the native zone
134        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
135
136        assert zone == 50
137        assert num.allclose(easting, 213251.040253)
138        assert num.allclose(northing, 6762559.15978)
139
140        #Then project to zone 49
141        zone, easting, northing = redfearn(-29.233299999,114.05,zone=49)
142
143        assert zone == 49
144        assert num.allclose(easting, 796474.020057)
145        assert num.allclose(northing, 6762310.25162)
146
147       
148
149       
150
151        #First test native projection (zone 49)
152        zone, easting, northing = redfearn(-29.1333,113.9667)
153
154        assert zone == 49
155        assert num.allclose(easting, 788653.192779)
156        assert num.allclose(northing, 6773605.46384)
157
158        #Then project to zone 50
159        zone, easting, northing = redfearn(-29.1333,113.9667,zone=50)
160
161        assert zone == 50
162        assert num.allclose(easting, 204863.606467)
163        assert num.allclose(northing, 6773440.04726)
164
165        #Testing point on zone boundary
166        #First test native projection (zone 50)
167        zone, easting, northing = redfearn(-29.1667,114)
168
169        assert zone == 50
170        assert num.allclose(easting, 208199.768268)
171        assert num.allclose(northing, 6769820.01453)
172
173        #Then project to zone 49
174        zone, easting, northing = redfearn(-29.1667,114,zone=49)
175
176        assert zone == 49
177        assert num.allclose(easting, 791800.231817)
178        assert num.allclose(northing, 6769820.01453)
179
180        #Testing furthest point in Geraldton scenario)
181        #First test native projection (zone 49)
182        zone, easting, northing = redfearn(-28.2167,113.4167)
183
184        assert zone == 49
185        assert num.allclose(easting, 737178.16131)
186        assert num.allclose(northing, 6876426.38578)
187
188        #Then project to zone 50
189        zone, easting, northing = redfearn(-28.2167,113.4167,zone=50)
190
191        assert zone == 50
192        assert num.allclose(easting, 148260.567427)
193        assert num.allclose(northing, 6873587.50926)
194
195        #Testing outside GDA zone (New Zeland)
196        #First test native projection (zone 60)
197        zone, easting, northing = redfearn(-44,178)
198
199        assert zone == 60
200        assert num.allclose(easting, 580174.259843)
201        assert num.allclose(northing, 5127641.114461)
202
203        #Then project to zone 59
204        zone, easting, northing = redfearn(-44,178,zone=59)
205
206        assert zone == 59
207        assert num.allclose(easting, 1061266.922118)
208        assert num.allclose(northing, 5104249.395469)
209
210        #Then skip three zones 57 (native 60)
211        zone, easting, northing = redfearn(-44,178,zone=57)
212
213        assert zone == 57
214        assert num.allclose(easting, 2023865.527497)
215        assert num.allclose(northing, 4949253.934967)
216
217        # Note Projecting into the Northern Hemisphere Does not coincide
218        # redfearn or ArcMap conversions. The difference lies in
219        # False Northing which is 0 for the Northern hemisphere
220        # in the redfearn implementation but 10 million in Arc.
221        #
222        # But the redfearn implementation does coincide with
223        # Google Earth (he he)
224       
225        #Testing outside GDA zone (Northern Hemisphere)
226        #First test native projection (zone 57)
227        zone, easting, northing = redfearn(44,156)
228
229        # Google Earth interpretation
230        assert zone == 57
231        assert num.allclose(easting, 259473.69)
232        assert num.allclose(northing, 4876249.13)
233
234        # ArcMap's interpretation
235        #assert zone == 57       
236        #assert num.allclose(easting, 259473.678944)
237        #assert num.allclose(northing, 14876249.1268)
238       
239        #Then project to zone 56
240        zone, easting, northing = redfearn(44,156,zone=56)
241
242        assert zone == 56
243        assert num.allclose(easting, 740526.321055)
244        assert num.allclose(northing, 4876249.13)
245
246       
247
248
249    #def test_UTM_6(self):
250    #    """Test 6 (Don's Wollongong file's ref point)
251    #    """
252    #
253    #    lat = -34.490286785873
254    #    lon = 150.79712139578
255    #
256    #    dd,mm,ss = decimal_degrees2degminsec(lat)
257    #    print dd,mm,ss
258    #    dd,mm,ss = decimal_degrees2degminsec(lon)       
259    #    print dd,mm,ss
260    #     
261    #    zone, easting, northing = redfearn(lat,lon)
262    #
263    #    print zone, easting, northing
264    #
265    #    assert zone == 56
266    #    #assert allclose(easting, 297717.36468927) #out by 10m
267    #    #assert allclose(northing, 6181725.1724276)
268
269    def test_convert_lats_longs(self):
270
271        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
272        #Zone:   56   
273        #Easting:  222908.705  Northing: 6233785.284
274        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
275        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
276
277        lat_gong = degminsec2decimal_degrees(-34,30,0.)
278        lon_gong = degminsec2decimal_degrees(150,55,0.)
279       
280        lat_2 = degminsec2decimal_degrees(-34,00,0.)
281        lon_2 = degminsec2decimal_degrees(150,00,0.)
282       
283        lats = [lat_gong, lat_2]
284        longs = [lon_gong, lon_2]
285        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
286
287        assert num.allclose(points[0][0], 308728.009)
288        assert num.allclose(points[0][1], 6180432.601)
289        assert num.allclose(points[1][0],  222908.705)
290        assert num.allclose(points[1][1], 6233785.284)
291        self.failUnless(zone == 56,
292                        'Bad zone error!')
293       
294    def test_convert_lats_longs2(self):
295
296        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
297        #Zone:   56   
298        #Easting:  222908.705  Northing: 6233785.284
299        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
300        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
301
302        lat_gong = degminsec2decimal_degrees(-34,30,0.)
303        lon_gong = degminsec2decimal_degrees(150,55,0.)
304       
305        lat_2 = degminsec2decimal_degrees(34,00,0.)
306        lon_2 = degminsec2decimal_degrees(100,00,0.)
307       
308        lats = [lat_gong, lat_2]
309        longs = [lon_gong, lon_2]
310       
311        try:
312            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
313        except ANUGAError:
314            pass
315        else:
316            self.failUnless(False,
317                            'Error not thrown error!')
318           
319    def test_convert_lats_longs3(self):
320
321        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
322        #Zone:   56   
323        #Easting:  222908.705  Northing: 6233785.284
324        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
325        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
326
327        lat_gong = "-34.5"
328        lon_gong = "150.916666667"
329        lat_2 = degminsec2decimal_degrees(34,00,0.)
330        lon_2 = degminsec2decimal_degrees(100,00,0.)
331       
332        lats = [lat_gong, lat_2]
333        longs = [lon_gong, lon_2]
334        try:
335            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
336        except ANUGAError:
337            pass
338        else:
339            self.failUnless(False,
340                            'Error not thrown error!')
341
342    def test_convert_latlon_to_UTM1(self):
343
344        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
345        #Zone:   56   
346        #Easting:  222908.705  Northing: 6233785.284
347        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
348        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
349
350        lat_gong = degminsec2decimal_degrees(-34,30,0.)
351        lon_gong = degminsec2decimal_degrees(150,55,0.)
352       
353        lat_2 = degminsec2decimal_degrees(-34,00,0.)
354        lon_2 = degminsec2decimal_degrees(150,00,0.)
355       
356        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
357        points, zone = convert_from_latlon_to_utm(points=points)
358        #print "points",points
359        assert num.allclose(points[0][0], 308728.009)
360        assert num.allclose(points[0][1], 6180432.601)
361        assert num.allclose(points[1][0],  222908.705)
362        assert num.allclose(points[1][1], 6233785.284)
363        self.failUnless(zone == 56,
364                        'Bad zone error!')
365
366    def test_convert_latlon_to_UTM2(self):       
367
368        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
369        #Zone:   56   
370        #Easting:  222908.705  Northing: 6233785.284
371        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
372        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
373
374        lat_gong = degminsec2decimal_degrees(-34,30,0.)
375        lon_gong = degminsec2decimal_degrees(150,55,0.)
376       
377        lat_2 = degminsec2decimal_degrees(34,00,0.)
378        lon_2 = degminsec2decimal_degrees(100,00,0.)
379
380        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
381
382        try:
383            points, zone = convert_from_latlon_to_utm(points=points)           
384        except ANUGAError:
385            pass
386        else:
387            self.fail('Error not thrown error!')
388
389    def test_convert_latlon_to_UTM3(self):           
390
391        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
392        #Zone:   56   
393        #Easting:  222908.705  Northing: 6233785.284
394        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
395        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
396
397        lat_gong = "-34.5"
398        lon_gong = "150.916666667"
399        lat_2 = degminsec2decimal_degrees(34,00,0.)
400        lon_2 = degminsec2decimal_degrees(100,00,0.)
401
402        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
403
404        try:
405            points, zone = convert_from_latlon_to_utm(points=points)           
406        except ANUGAError:
407            pass
408        else:
409            self.fail('Error not thrown error!')
410
411           
412#-------------------------------------------------------------
413if __name__ == "__main__":
414
415    #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
416    #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')
417    mysuite = unittest.makeSuite(TestCase,'test')
418    runner = unittest.TextTestRunner()
419    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.