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

Last change on this file since 6440 was 6440, checked in by ole, 14 years ago

Added more test points for nonstandard projections and a new test.
Both non standard project tests are disabled and don't pass fully.

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