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

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

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