source: branches/numpy/anuga/coordinate_transforms/test_redfearn.py @ 7195

Last change on this file since 7195 was 6553, checked in by rwilson, 16 years ago

Merged trunk into numpy, all tests and validations work.
File size: 17.0 KB
RevLine 
[6149]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
[6553]13from anuga.utilities.system_tools import get_pathname_from_package
14from os.path import join
[6304]15import numpy as num
[6149]16
17
18#-------------------------------------------------------------
19
20class TestCase(unittest.TestCase):
21
22    def test_decimal_degrees_conversion(Self):
23        lat = degminsec2decimal_degrees(-37,39,10.15610)
24        lon = degminsec2decimal_degrees(143,55,35.38390) 
25        assert num.allclose(lat, -37.65282114)
26        assert num.allclose(lon, 143.9264955)
27
28        dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
29        assert dd==-37
30        assert mm==39
31        assert num.allclose(ss, 10.15610)
32
33        dd,mm,ss = decimal_degrees2degminsec(143.9264955)
34        assert dd==143
35        assert mm==55
36        assert num.allclose(ss, 35.38390) 
37
38
39    def test_UTM_1(self):
40        #latitude:  -37 39' 10.15610"
41        #Longitude: 143 55' 35.38390"
42        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
43        #Zone:   54   
44        #Easting:  758173.797  Northing: 5828674.340
45        #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 ''
46        #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107
47
48        lat = degminsec2decimal_degrees(-37,39,10.15610)
49        lon = degminsec2decimal_degrees(143,55,35.38390) 
50        assert num.allclose(lat, -37.65282114)
51        assert num.allclose(lon, 143.9264955)
52
53
54        zone, easting, northing = redfearn(lat,lon)
55
56        assert zone == 54
57        assert num.allclose(easting, 758173.797)
58        assert num.allclose(northing, 5828674.340)
59
60
61    def test_UTM_2(self):
62        #TEST 2
63
64        #Latitude:  -37 57 03.7203
65        #Longitude: 144 25 29.5244
66        #Zone:   55   
67        #Easting:  273741.297  Northing: 5796489.777
68        #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 ''
69        #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056
70
71        lat = degminsec2decimal_degrees(-37,57,03.7203)
72        lon = degminsec2decimal_degrees(144,25,29.5244) 
73        #print lat, lon
74
75        zone, easting, northing = redfearn(lat,lon)
76
77        assert zone == 55
78        assert num.allclose(easting, 273741.297)
79        assert num.allclose(northing, 5796489.777)
80
81       
82    def test_UTM_3(self):
83        #Test 3
84        lat = degminsec2decimal_degrees(-60,0,0)
85        lon = degminsec2decimal_degrees(130,0,0) 
86
87        zone, easting, northing = redfearn(lat,lon)
88        #print zone, easting, northing
89
90        assert zone == 52
91        assert num.allclose(easting, 555776.267)
92        assert num.allclose(northing, 3348167.264)
93
94
95    def test_UTM_4(self):
96        #Test 4 (Kobenhavn, Northern hemisphere)
97        lat = 55.70248
98        dd,mm,ss = decimal_degrees2degminsec(lat)
99
100        lon = 12.58364
101        dd,mm,ss = decimal_degrees2degminsec(lon)
102
103        zone, easting, northing = redfearn(lat,lon)
104
105
106        assert zone == 33
107        assert num.allclose(easting, 348157.631)
108        assert num.allclose(northing, 6175612.993) 
109
110
111    def test_UTM_5(self):
112        #Test 5 (Wollongong)
113
114        lat = degminsec2decimal_degrees(-34,30,0.)
115        lon = degminsec2decimal_degrees(150,55,0.) 
116       
117        zone, easting, northing = redfearn(lat,lon)
118
119        #print zone, easting, northing
120
121        assert zone == 56
122        assert num.allclose(easting, 308728.009)
123        assert num.allclose(northing, 6180432.601)
124
125    def test_UTM_6_nonstandard_projection(self):
[6553]126        """test_UTM_6_nonstandard_projection
[6149]127
[6553]128        Test that projections can be forced to
129        use other than native zone.
130
131        Data is from Geraldton, WA
132        """
133       
134
135        # First test native projection (zone 50)
[6149]136        zone, easting, northing = redfearn(-29.233299999,114.05)
137
138        assert zone == 50
139        assert num.allclose(easting, 213251.040253)
140        assert num.allclose(northing, 6762559.15978)
141
142        #Testing using the native zone
143        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
144
145        assert zone == 50
146        assert num.allclose(easting, 213251.040253)
147        assert num.allclose(northing, 6762559.15978)
148
149        #Then project to zone 49
150        zone, easting, northing = redfearn(-29.233299999,114.05,zone=49)
151
152        assert zone == 49
153        assert num.allclose(easting, 796474.020057)
154        assert num.allclose(northing, 6762310.25162)
155
156       
157
158       
159
160        #First test native projection (zone 49)
161        zone, easting, northing = redfearn(-29.1333,113.9667)
162
163        assert zone == 49
164        assert num.allclose(easting, 788653.192779)
165        assert num.allclose(northing, 6773605.46384)
166
167        #Then project to zone 50
168        zone, easting, northing = redfearn(-29.1333,113.9667,zone=50)
169
170        assert zone == 50
171        assert num.allclose(easting, 204863.606467)
172        assert num.allclose(northing, 6773440.04726)
173
174        #Testing point on zone boundary
175        #First test native projection (zone 50)
176        zone, easting, northing = redfearn(-29.1667,114)
177
178        assert zone == 50
179        assert num.allclose(easting, 208199.768268)
180        assert num.allclose(northing, 6769820.01453)
181
182        #Then project to zone 49
183        zone, easting, northing = redfearn(-29.1667,114,zone=49)
184
185        assert zone == 49
186        assert num.allclose(easting, 791800.231817)
187        assert num.allclose(northing, 6769820.01453)
188
189        #Testing furthest point in Geraldton scenario)
190        #First test native projection (zone 49)
191        zone, easting, northing = redfearn(-28.2167,113.4167)
192
193        assert zone == 49
194        assert num.allclose(easting, 737178.16131)
195        assert num.allclose(northing, 6876426.38578)
196
197        #Then project to zone 50
198        zone, easting, northing = redfearn(-28.2167,113.4167,zone=50)
199
200        assert zone == 50
201        assert num.allclose(easting, 148260.567427)
202        assert num.allclose(northing, 6873587.50926)
203
204        #Testing outside GDA zone (New Zeland)
205        #First test native projection (zone 60)
206        zone, easting, northing = redfearn(-44,178)
207
208        assert zone == 60
209        assert num.allclose(easting, 580174.259843)
210        assert num.allclose(northing, 5127641.114461)
211
212        #Then project to zone 59
213        zone, easting, northing = redfearn(-44,178,zone=59)
214
215        assert zone == 59
216        assert num.allclose(easting, 1061266.922118)
217        assert num.allclose(northing, 5104249.395469)
218
219        #Then skip three zones 57 (native 60)
220        zone, easting, northing = redfearn(-44,178,zone=57)
221
222        assert zone == 57
223        assert num.allclose(easting, 2023865.527497)
224        assert num.allclose(northing, 4949253.934967)
225
226        # Note Projecting into the Northern Hemisphere Does not coincide
227        # redfearn or ArcMap conversions. The difference lies in
228        # False Northing which is 0 for the Northern hemisphere
229        # in the redfearn implementation but 10 million in Arc.
230        #
231        # But the redfearn implementation does coincide with
232        # Google Earth (he he)
233       
234        #Testing outside GDA zone (Northern Hemisphere)
235        #First test native projection (zone 57)
236        zone, easting, northing = redfearn(44,156)
237
238        # Google Earth interpretation
239        assert zone == 57
240        assert num.allclose(easting, 259473.69)
241        assert num.allclose(northing, 4876249.13)
242
243        # ArcMap's interpretation
244        #assert zone == 57       
245        #assert num.allclose(easting, 259473.678944)
246        #assert num.allclose(northing, 14876249.1268)
247       
248        #Then project to zone 56
249        zone, easting, northing = redfearn(44,156,zone=56)
250
251        assert zone == 56
252        assert num.allclose(easting, 740526.321055)
253        assert num.allclose(northing, 4876249.13)
254
255       
256
257
258    #def test_UTM_6(self):
259    #    """Test 6 (Don's Wollongong file's ref point)
260    #    """
261    #
262    #    lat = -34.490286785873
263    #    lon = 150.79712139578
264    #
265    #    dd,mm,ss = decimal_degrees2degminsec(lat)
266    #    print dd,mm,ss
267    #    dd,mm,ss = decimal_degrees2degminsec(lon)       
268    #    print dd,mm,ss
269    #     
270    #    zone, easting, northing = redfearn(lat,lon)
271    #
272    #    print zone, easting, northing
273    #
274    #    assert zone == 56
275    #    #assert allclose(easting, 297717.36468927) #out by 10m
276    #    #assert allclose(northing, 6181725.1724276)
277
[6553]278   
279    def test_nonstandard_meridian_coinciding_with_native(self):
280        """test_nonstandard_meridian_coinciding_with_native
281
282        This test will verify that redfearn can be used to project
283        points using an arbitrary central meridian that happens to
284        coincide with the standard meridian at the center of a UTM zone.
285        This is a preliminary test before testing this functionality
286        with a truly arbitrary non-standard meridian.
287        """
288
289        # The file projection_test_points_z53.csv contains 10 points
290        # which straddle the boundary between UTM zones 53 and 54.
291        # They have been projected to zone 53 irrespective of where they
292        # belong.
293
294        path = get_pathname_from_package('anuga.coordinate_transforms')
295       
296        for forced_zone in [53, 54]:
297       
298            datafile = join(path, 'projection_test_points_z%d.csv' % forced_zone)
299            fid = open(datafile)
300
301            for line in fid.readlines()[1:]:
302                fields = line.strip().split(',')
303               
304                lon = float(fields[1])
305                lat = float(fields[2])
306                x = float(fields[3])
307                y = float(fields[4])           
308
309                zone, easting, northing = redfearn(lat, lon,
310                                                   zone=forced_zone)
311               
312                # Check calculation
313                assert zone == forced_zone
314                assert num.allclose(x, easting)
315                assert num.allclose(y, northing)
316
317   
318   
319   
320    def test_nonstandard_meridian(self):
321        """test_nonstandard_meridian
322
323        This test will verify that redfearn can be used to project
324        points using an arbitrary central meridian.
325        """
326
327        # The file projection_test_points.csv contains 10 points
328        # which straddle the boundary between UTM zones 53 and 54.
329        # They have been projected using a central meridian of 137.5
330        # degrees (the boundary is 138 so it is pretty much right
331        # in the middle of zones 53 and 54).
332
333        path = get_pathname_from_package('anuga.coordinate_transforms')
334        datafile = join(path, 'projection_test_points.csv')
335        fid = open(datafile)
336
337        for line in fid.readlines()[1:]:
338            fields = line.strip().split(',')
339
340            lon = float(fields[1])
341            lat = float(fields[2])
342            x = float(fields[3])
343            y = float(fields[4])           
344
345            zone, easting, northing = redfearn(lat, lon,
346                                               central_meridian=137.5,
347                                               scale_factor=0.9996)
348
349            assert zone == -1 # Indicates non UTM projection
350            assert num.allclose(x, easting)
351            assert num.allclose(y, northing)
352
353        # Test that zone and meridian can't both be specified
354        try:
355            zone, easting, northing = redfearn(lat, lon,
356                                               zone=50, 
357                                               central_meridian=137.5)
358        except:
359            pass
360        else:
361            msg = 'Should have raised exception'
362            raise Exception, msg
363
364           
[6149]365    def test_convert_lats_longs(self):
366
367        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
368        #Zone:   56   
369        #Easting:  222908.705  Northing: 6233785.284
370        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
371        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
372
373        lat_gong = degminsec2decimal_degrees(-34,30,0.)
374        lon_gong = degminsec2decimal_degrees(150,55,0.)
375       
376        lat_2 = degminsec2decimal_degrees(-34,00,0.)
377        lon_2 = degminsec2decimal_degrees(150,00,0.)
378       
379        lats = [lat_gong, lat_2]
380        longs = [lon_gong, lon_2]
381        points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
382
383        assert num.allclose(points[0][0], 308728.009)
384        assert num.allclose(points[0][1], 6180432.601)
385        assert num.allclose(points[1][0],  222908.705)
386        assert num.allclose(points[1][1], 6233785.284)
387        self.failUnless(zone == 56,
388                        'Bad zone error!')
389       
390    def test_convert_lats_longs2(self):
391
392        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
393        #Zone:   56   
394        #Easting:  222908.705  Northing: 6233785.284
395        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
396        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
397
398        lat_gong = degminsec2decimal_degrees(-34,30,0.)
399        lon_gong = degminsec2decimal_degrees(150,55,0.)
400       
401        lat_2 = degminsec2decimal_degrees(34,00,0.)
402        lon_2 = degminsec2decimal_degrees(100,00,0.)
403       
404        lats = [lat_gong, lat_2]
405        longs = [lon_gong, lon_2]
406       
407        try:
408            points, zone = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
409        except ANUGAError:
410            pass
411        else:
412            self.failUnless(False,
413                            'Error not thrown error!')
414           
415    def test_convert_lats_longs3(self):
416
417        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
418        #Zone:   56   
419        #Easting:  222908.705  Northing: 6233785.284
420        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
421        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
422
423        lat_gong = "-34.5"
424        lon_gong = "150.916666667"
425        lat_2 = degminsec2decimal_degrees(34,00,0.)
426        lon_2 = degminsec2decimal_degrees(100,00,0.)
427       
428        lats = [lat_gong, lat_2]
429        longs = [lon_gong, lon_2]
430        try:
431            points, zone  = convert_from_latlon_to_utm(latitudes=lats, longitudes=longs)
432        except ANUGAError:
433            pass
434        else:
435            self.failUnless(False,
436                            'Error not thrown error!')
437
438    def test_convert_latlon_to_UTM1(self):
439
440        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
441        #Zone:   56   
442        #Easting:  222908.705  Northing: 6233785.284
443        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
444        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
445
446        lat_gong = degminsec2decimal_degrees(-34,30,0.)
447        lon_gong = degminsec2decimal_degrees(150,55,0.)
448       
449        lat_2 = degminsec2decimal_degrees(-34,00,0.)
450        lon_2 = degminsec2decimal_degrees(150,00,0.)
451       
452        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
453        points, zone = convert_from_latlon_to_utm(points=points)
454        #print "points",points
455        assert num.allclose(points[0][0], 308728.009)
456        assert num.allclose(points[0][1], 6180432.601)
457        assert num.allclose(points[1][0],  222908.705)
458        assert num.allclose(points[1][1], 6233785.284)
459        self.failUnless(zone == 56,
460                        'Bad zone error!')
461
462    def test_convert_latlon_to_UTM2(self):       
463
464        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
465        #Zone:   56   
466        #Easting:  222908.705  Northing: 6233785.284
467        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
468        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
469
470        lat_gong = degminsec2decimal_degrees(-34,30,0.)
471        lon_gong = degminsec2decimal_degrees(150,55,0.)
472       
473        lat_2 = degminsec2decimal_degrees(34,00,0.)
474        lon_2 = degminsec2decimal_degrees(100,00,0.)
475
476        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
477
478        try:
479            points, zone = convert_from_latlon_to_utm(points=points)           
480        except ANUGAError:
481            pass
482        else:
483            self.fail('Error not thrown error!')
484
485    def test_convert_latlon_to_UTM3(self):           
486
487        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
488        #Zone:   56   
489        #Easting:  222908.705  Northing: 6233785.284
490        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
491        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
492
493        lat_gong = "-34.5"
494        lon_gong = "150.916666667"
495        lat_2 = degminsec2decimal_degrees(34,00,0.)
496        lon_2 = degminsec2decimal_degrees(100,00,0.)
497
498        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
499
500        try:
501            points, zone = convert_from_latlon_to_utm(points=points)           
502        except ANUGAError:
503            pass
504        else:
505            self.fail('Error not thrown error!')
506
507           
508#-------------------------------------------------------------
[6410]509
[6149]510if __name__ == "__main__":
511    mysuite = unittest.makeSuite(TestCase,'test')
512    runner = unittest.TextTestRunner()
513    runner.run(mysuite)
Note: See TracBrowser for help on using the repository browser.