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

Last change on this file since 5662 was 5662, checked in by kristy, 16 years ago

Update so that we can pass a different zone into the redfearn equation (test_UTM_6_nonstandard_projection)
File size: 13.5 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
217##        #Testing outside GDA zone (Northern Hemisphere)
218##        #First test native projection (zone 57)
219##        zone, easting, northing = redfearn(44,156)
220##
221##        assert zone == 57
222##        assert allclose(easting, 259473.678944)
223##        assert allclose(northing, 14876249.1268)
224##
225##        #Then project to zone 59
226##        zone, easting, northing = redfearn(44,156,zone=56)
227##
228##        assert zone == 56
229##        assert allclose(easting, 740526.321055)
230##        assert allclose(northing, 14876249.1268)
231
232
233       
234
235
236    #def test_UTM_6(self):
237    #    """Test 6 (Don's Wollongong file's ref point)
238    #    """
239    #
240    #    lat = -34.490286785873
241    #    lon = 150.79712139578
242    #
243    #    dd,mm,ss = decimal_degrees2degminsec(lat)
244    #    print dd,mm,ss
245    #    dd,mm,ss = decimal_degrees2degminsec(lon)       
246    #    print dd,mm,ss
247    #     
248    #    zone, easting, northing = redfearn(lat,lon)
249    #
250    #    print zone, easting, northing
251    #
252    #    assert zone == 56
253    #    #assert allclose(easting, 297717.36468927) #out by 10m
254    #    #assert allclose(northing, 6181725.1724276)
255
256    def test_convert_lats_longs(self):
257
258        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
259        #Zone:   56   
260        #Easting:  222908.705  Northing: 6233785.284
261        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
262        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
263
264        lat_gong = degminsec2decimal_degrees(-34,30,0.)
265        lon_gong = degminsec2decimal_degrees(150,55,0.)
266       
267        lat_2 = degminsec2decimal_degrees(-34,00,0.)
268        lon_2 = degminsec2decimal_degrees(150,00,0.)
269       
270        lats = [lat_gong, lat_2]
271        longs = [lon_gong, lon_2]
272        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
273
274        assert allclose(points[0][0], 308728.009)
275        assert allclose(points[0][1], 6180432.601)
276        assert allclose(points[1][0],  222908.705)
277        assert allclose(points[1][1], 6233785.284)
278        self.failUnless(zone == 56,
279                        'Bad zone error!')
280       
281    def test_convert_lats_longs2(self):
282
283        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
284        #Zone:   56   
285        #Easting:  222908.705  Northing: 6233785.284
286        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
287        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
288
289        lat_gong = degminsec2decimal_degrees(-34,30,0.)
290        lon_gong = degminsec2decimal_degrees(150,55,0.)
291       
292        lat_2 = degminsec2decimal_degrees(34,00,0.)
293        lon_2 = degminsec2decimal_degrees(100,00,0.)
294       
295        lats = [lat_gong, lat_2]
296        longs = [lon_gong, lon_2]
297       
298        try:
299            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
300        except ANUGAError:
301            pass
302        else:
303            self.failUnless(False,
304                            'Error not thrown error!')
305           
306    def test_convert_lats_longs3(self):
307
308        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
309        #Zone:   56   
310        #Easting:  222908.705  Northing: 6233785.284
311        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
312        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
313
314        lat_gong = "-34.5"
315        lon_gong = "150.916666667"
316        lat_2 = degminsec2decimal_degrees(34,00,0.)
317        lon_2 = degminsec2decimal_degrees(100,00,0.)
318       
319        lats = [lat_gong, lat_2]
320        longs = [lon_gong, lon_2]
321        try:
322            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
323        except ANUGAError:
324            pass
325        else:
326            self.failUnless(False,
327                            'Error not thrown error!')
328
329    def test_convert_latlon_to_UTM1(self):
330
331        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
332        #Zone:   56   
333        #Easting:  222908.705  Northing: 6233785.284
334        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
335        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
336
337        lat_gong = degminsec2decimal_degrees(-34,30,0.)
338        lon_gong = degminsec2decimal_degrees(150,55,0.)
339       
340        lat_2 = degminsec2decimal_degrees(-34,00,0.)
341        lon_2 = degminsec2decimal_degrees(150,00,0.)
342       
343        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
344        points, zone = convert_from_latlon_to_utm(points=points)
345        #print "points",points
346        assert allclose(points[0][0], 308728.009)
347        assert allclose(points[0][1], 6180432.601)
348        assert allclose(points[1][0],  222908.705)
349        assert allclose(points[1][1], 6233785.284)
350        self.failUnless(zone == 56,
351                        'Bad zone error!')
352
353    def test_convert_latlon_to_UTM2(self):       
354
355        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
356        #Zone:   56   
357        #Easting:  222908.705  Northing: 6233785.284
358        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
359        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
360
361        lat_gong = degminsec2decimal_degrees(-34,30,0.)
362        lon_gong = degminsec2decimal_degrees(150,55,0.)
363       
364        lat_2 = degminsec2decimal_degrees(34,00,0.)
365        lon_2 = degminsec2decimal_degrees(100,00,0.)
366
367        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
368
369        try:
370            points, zone = convert_from_latlon_to_utm(points=points)           
371        except ANUGAError:
372            pass
373        else:
374            self.fail('Error not thrown error!')
375
376    def test_convert_latlon_to_UTM3(self):           
377
378        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
379        #Zone:   56   
380        #Easting:  222908.705  Northing: 6233785.284
381        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
382        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
383
384        lat_gong = "-34.5"
385        lon_gong = "150.916666667"
386        lat_2 = degminsec2decimal_degrees(34,00,0.)
387        lon_2 = degminsec2decimal_degrees(100,00,0.)
388
389        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
390
391        try:
392            points, zone = convert_from_latlon_to_utm(points=points)           
393        except ANUGAError:
394            pass
395        else:
396            self.fail('Error not thrown error!')
397
398           
399#-------------------------------------------------------------
400if __name__ == "__main__":
401
402    #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
403    #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')
404    mysuite = unittest.makeSuite(TestCase,'test')
405    runner = unittest.TextTestRunner()
406    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.