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

Last change on this file since 6400 was 6400, checked in by ole, 15 years ago

Added data from Leharne for use with arbitrary meridian.

File size: 14.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
9
10from redfearn import *
11from anuga.utilities.anuga_exceptions import ANUGAError
12
13import Numeric 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_UTM_6_nonstandard_projection
125
126        Test that projections can be forced to
127        use other than native zone.
128
129        Data is from Geraldton, WA
130        """
131       
132
133        # First test native projection (zone 50)
134        zone, easting, northing = redfearn(-29.233299999,114.05)
135
136        assert zone == 50
137        assert num.allclose(easting, 213251.040253)
138        assert num.allclose(northing, 6762559.15978)
139
140        # Testing using the native zone
141        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
142
143        assert zone == 50
144        assert num.allclose(easting, 213251.040253)
145        assert num.allclose(northing, 6762559.15978)
146
147        # Then project to zone 49
148        zone, easting, northing = redfearn(-29.233299999,114.05,zone=49)
149
150        assert zone == 49
151        assert num.allclose(easting, 796474.020057)
152        assert num.allclose(northing, 6762310.25162)
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(self):
276        """test_nonstandard_meridian
277
278        This test will verify that redfearn can be used to project
279        points using an arbitrary central meridian.
280        """
281
282
283        # FIXME: To do using csv file
284        pass
285
286
287        # The file projection_test_points.csv contains 10 points
288        # which straddle the boundary between UTM zones 53 and 54.
289        # They have been projected using a central meridian of 137.5
290        # degrees (the boundary is 138 so it is pretty much right
291        # in the middle of zones 53 and 54).
292
293       
294   
295   
296
297    def test_convert_lats_longs(self):
298
299        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
300        #Zone:   56   
301        #Easting:  222908.705  Northing: 6233785.284
302        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
303        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
304
305        lat_gong = degminsec2decimal_degrees(-34,30,0.)
306        lon_gong = degminsec2decimal_degrees(150,55,0.)
307       
308        lat_2 = degminsec2decimal_degrees(-34,00,0.)
309        lon_2 = degminsec2decimal_degrees(150,00,0.)
310       
311        lats = [lat_gong, lat_2]
312        longs = [lon_gong, lon_2]
313        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
314
315        assert num.allclose(points[0][0], 308728.009)
316        assert num.allclose(points[0][1], 6180432.601)
317        assert num.allclose(points[1][0],  222908.705)
318        assert num.allclose(points[1][1], 6233785.284)
319        self.failUnless(zone == 56,
320                        'Bad zone error!')
321       
322    def test_convert_lats_longs2(self):
323
324        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
325        #Zone:   56   
326        #Easting:  222908.705  Northing: 6233785.284
327        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
328        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
329
330        lat_gong = degminsec2decimal_degrees(-34,30,0.)
331        lon_gong = degminsec2decimal_degrees(150,55,0.)
332       
333        lat_2 = degminsec2decimal_degrees(34,00,0.)
334        lon_2 = degminsec2decimal_degrees(100,00,0.)
335       
336        lats = [lat_gong, lat_2]
337        longs = [lon_gong, lon_2]
338       
339        try:
340            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
341        except ANUGAError:
342            pass
343        else:
344            self.failUnless(False,
345                            'Error not thrown error!')
346           
347    def test_convert_lats_longs3(self):
348
349        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
350        #Zone:   56   
351        #Easting:  222908.705  Northing: 6233785.284
352        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
353        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
354
355        lat_gong = "-34.5"
356        lon_gong = "150.916666667"
357        lat_2 = degminsec2decimal_degrees(34,00,0.)
358        lon_2 = degminsec2decimal_degrees(100,00,0.)
359       
360        lats = [lat_gong, lat_2]
361        longs = [lon_gong, lon_2]
362        try:
363            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
364        except ANUGAError:
365            pass
366        else:
367            self.failUnless(False,
368                            'Error not thrown error!')
369
370    def test_convert_latlon_to_UTM1(self):
371
372        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
373        #Zone:   56   
374        #Easting:  222908.705  Northing: 6233785.284
375        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
376        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
377
378        lat_gong = degminsec2decimal_degrees(-34,30,0.)
379        lon_gong = degminsec2decimal_degrees(150,55,0.)
380       
381        lat_2 = degminsec2decimal_degrees(-34,00,0.)
382        lon_2 = degminsec2decimal_degrees(150,00,0.)
383       
384        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
385        points, zone = convert_from_latlon_to_utm(points=points)
386        #print "points",points
387        assert num.allclose(points[0][0], 308728.009)
388        assert num.allclose(points[0][1], 6180432.601)
389        assert num.allclose(points[1][0],  222908.705)
390        assert num.allclose(points[1][1], 6233785.284)
391        self.failUnless(zone == 56,
392                        'Bad zone error!')
393
394    def test_convert_latlon_to_UTM2(self):       
395
396        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
397        #Zone:   56   
398        #Easting:  222908.705  Northing: 6233785.284
399        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
400        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
401
402        lat_gong = degminsec2decimal_degrees(-34,30,0.)
403        lon_gong = degminsec2decimal_degrees(150,55,0.)
404       
405        lat_2 = degminsec2decimal_degrees(34,00,0.)
406        lon_2 = degminsec2decimal_degrees(100,00,0.)
407
408        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
409
410        try:
411            points, zone = convert_from_latlon_to_utm(points=points)           
412        except ANUGAError:
413            pass
414        else:
415            self.fail('Error not thrown error!')
416
417    def test_convert_latlon_to_UTM3(self):           
418
419        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
420        #Zone:   56   
421        #Easting:  222908.705  Northing: 6233785.284
422        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
423        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
424
425        lat_gong = "-34.5"
426        lon_gong = "150.916666667"
427        lat_2 = degminsec2decimal_degrees(34,00,0.)
428        lon_2 = degminsec2decimal_degrees(100,00,0.)
429
430        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
431
432        try:
433            points, zone = convert_from_latlon_to_utm(points=points)           
434        except ANUGAError:
435            pass
436        else:
437            self.fail('Error not thrown error!')
438
439           
440#-------------------------------------------------------------
441if __name__ == "__main__":
442
443    #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
444    #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')
445    mysuite = unittest.makeSuite(TestCase,'test')
446    runner = unittest.TextTestRunner()
447    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.