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

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

Experimented with scale factors in redfearn
File size: 16.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    def Xtest_nonstandard_meridian(self):
277        """test_nonstandard_meridian
278
279        This test will verify that redfearn can be used to project
280        points using an arbitrary central meridian.
281        """
282
283        # The file projection_test_points.csv contains 10 points
284        # which straddle the boundary between UTM zones 53 and 54.
285        # They have been projected using a central meridian of 137.5
286        # degrees (the boundary is 138 so it is pretty much right
287        # in the middle of zones 53 and 54).
288
289        path = get_pathname_from_package('anuga.coordinate_transforms')
290        datafile = join(path, 'projection_test_points.csv')
291        fid = open(datafile)
292
293        for line in fid.readlines()[1:]:
294            fields = line.strip().split(',')
295
296            lon = float(fields[1])
297            lat = float(fields[2])
298            x = float(fields[3])
299            y = float(fields[4])           
300
301            zone, easting, northing = redfearn(lat, lon,
302                                               central_meridian=137.5,
303                                               scale_factor=0.9998154)
304
305            print
306            print 'Lat', lat
307            print 'Lon', lon
308            print 'Zone', zone
309            print 'Ref x', x, 'Computed x', easting, 'Close enough:', num.allclose(x, easting)
310            print 'Ref y', y, 'Computed y', northing, 'Close enough:', num.allclose(y, northing)
311
312            # Check calculation
313            assert zone == -1 # Indicates non UTM projection
314            print 
315            #assert num.allclose(x, easting)
316            #assert num.allclose(y, northing)
317
318        # Test that zone and meridian can't both be specified
319        try:
320            zone, easting, northing = redfearn(lat, lon,
321                                               zone=50, 
322                                               central_meridian=137.5,
323                                               scale_factor=0.968)
324        except:
325            pass
326        else:
327            msg = 'Should have raised exception'
328            raise Exception, msg
329
330    def test_convert_lats_longs(self):
331
332        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
333        #Zone:   56   
334        #Easting:  222908.705  Northing: 6233785.284
335        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
336        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
337
338        lat_gong = degminsec2decimal_degrees(-34,30,0.)
339        lon_gong = degminsec2decimal_degrees(150,55,0.)
340       
341        lat_2 = degminsec2decimal_degrees(-34,00,0.)
342        lon_2 = degminsec2decimal_degrees(150,00,0.)
343       
344        lats = [lat_gong, lat_2]
345        longs = [lon_gong, lon_2]
346        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
347
348        assert num.allclose(points[0][0], 308728.009)
349        assert num.allclose(points[0][1], 6180432.601)
350        assert num.allclose(points[1][0],  222908.705)
351        assert num.allclose(points[1][1], 6233785.284)
352        self.failUnless(zone == 56,
353                        'Bad zone error!')
354       
355    def test_convert_lats_longs2(self):
356
357        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
358        #Zone:   56   
359        #Easting:  222908.705  Northing: 6233785.284
360        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
361        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
362
363        lat_gong = degminsec2decimal_degrees(-34,30,0.)
364        lon_gong = degminsec2decimal_degrees(150,55,0.)
365       
366        lat_2 = degminsec2decimal_degrees(34,00,0.)
367        lon_2 = degminsec2decimal_degrees(100,00,0.)
368       
369        lats = [lat_gong, lat_2]
370        longs = [lon_gong, lon_2]
371       
372        try:
373            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
374        except ANUGAError:
375            pass
376        else:
377            self.failUnless(False,
378                            'Error not thrown error!')
379           
380    def test_convert_lats_longs3(self):
381
382        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
383        #Zone:   56   
384        #Easting:  222908.705  Northing: 6233785.284
385        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
386        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
387
388        lat_gong = "-34.5"
389        lon_gong = "150.916666667"
390        lat_2 = degminsec2decimal_degrees(34,00,0.)
391        lon_2 = degminsec2decimal_degrees(100,00,0.)
392       
393        lats = [lat_gong, lat_2]
394        longs = [lon_gong, lon_2]
395        try:
396            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
397        except ANUGAError:
398            pass
399        else:
400            self.failUnless(False,
401                            'Error not thrown error!')
402
403    def test_convert_latlon_to_UTM1(self):
404
405        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
406        #Zone:   56   
407        #Easting:  222908.705  Northing: 6233785.284
408        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
409        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
410
411        lat_gong = degminsec2decimal_degrees(-34,30,0.)
412        lon_gong = degminsec2decimal_degrees(150,55,0.)
413       
414        lat_2 = degminsec2decimal_degrees(-34,00,0.)
415        lon_2 = degminsec2decimal_degrees(150,00,0.)
416       
417        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
418        points, zone = convert_from_latlon_to_utm(points=points)
419        #print "points",points
420        assert num.allclose(points[0][0], 308728.009)
421        assert num.allclose(points[0][1], 6180432.601)
422        assert num.allclose(points[1][0],  222908.705)
423        assert num.allclose(points[1][1], 6233785.284)
424        self.failUnless(zone == 56,
425                        'Bad zone error!')
426
427    def test_convert_latlon_to_UTM2(self):       
428
429        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
430        #Zone:   56   
431        #Easting:  222908.705  Northing: 6233785.284
432        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
433        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
434
435        lat_gong = degminsec2decimal_degrees(-34,30,0.)
436        lon_gong = degminsec2decimal_degrees(150,55,0.)
437       
438        lat_2 = degminsec2decimal_degrees(34,00,0.)
439        lon_2 = degminsec2decimal_degrees(100,00,0.)
440
441        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
442
443        try:
444            points, zone = convert_from_latlon_to_utm(points=points)           
445        except ANUGAError:
446            pass
447        else:
448            self.fail('Error not thrown error!')
449
450    def test_convert_latlon_to_UTM3(self):           
451
452        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
453        #Zone:   56   
454        #Easting:  222908.705  Northing: 6233785.284
455        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
456        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
457
458        lat_gong = "-34.5"
459        lon_gong = "150.916666667"
460        lat_2 = degminsec2decimal_degrees(34,00,0.)
461        lon_2 = degminsec2decimal_degrees(100,00,0.)
462
463        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
464
465        try:
466            points, zone = convert_from_latlon_to_utm(points=points)           
467        except ANUGAError:
468            pass
469        else:
470            self.fail('Error not thrown error!')
471
472           
473#-------------------------------------------------------------
474if __name__ == "__main__":
475
476    #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
477    #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')
478    mysuite = unittest.makeSuite(TestCase,'test')
479    runner = unittest.TextTestRunner()
480    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.