source: anuga_core/source/anuga/coordinate_transforms/test_geo_reference.py @ 7063

Last change on this file since 7063 was 7063, checked in by rwilson, 15 years ago

Fixed a 'numpy-ism' in Numeric code.

File size: 18.4 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import tempfile
6import os
7
8from geo_reference import *
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10
11import Numeric as num
12
13
14class geo_referenceTestCase(unittest.TestCase):
15    def setUp(self):
16        pass
17       
18    def tearDown(self):
19        pass
20
21   
22   
23    def test_get_origin(self):
24        g = Geo_reference(56,1.9,1.9)
25        (z,x,y) = g.get_origin()
26
27        self.failUnless(z == g.get_zone(), ' failed')
28        self.failUnless(x == g.get_xllcorner(), ' failed')
29        self.failUnless(y == g.get_yllcorner(), ' failed') 
30       
31    def test_read_write_NetCDF(self):
32        from Scientific.IO.NetCDF import NetCDFFile
33        g = Geo_reference(56,1.9,1.9)
34        file_name = tempfile.mktemp(".geo_referenceTest")
35       
36        out_file = NetCDFFile(file_name, netcdf_mode_w)
37        g.write_NetCDF(out_file)
38        out_file.close()
39       
40        in_file = NetCDFFile(file_name, netcdf_mode_r)
41        new_g = Geo_reference(NetCDFObject=in_file)
42        in_file.close()
43        os.remove(file_name)
44
45        self.failUnless(g == new_g, 'test_read_write_NetCDF failed') 
46       
47    def test_read_NetCDFI(self):
48        # test if read_NetCDF
49        from Scientific.IO.NetCDF import NetCDFFile
50        g = Geo_reference(56,1.9,1.9)
51        file_name = tempfile.mktemp(".geo_referenceTest")
52       
53        outfile = NetCDFFile(file_name, netcdf_mode_w)
54        outfile.xllcorner = g.get_xllcorner() 
55        outfile.yllcorner =  g.get_yllcorner() 
56        outfile.zone = g.get_zone()
57        outfile.close()
58       
59        in_file = NetCDFFile(file_name, netcdf_mode_r)
60        new_g = Geo_reference(NetCDFObject=in_file)
61        in_file.close()
62        os.remove(file_name)
63
64        self.failUnless(g == new_g, ' failed')
65       
66    def test_read_write_ASCII(self):
67        from Scientific.IO.NetCDF import NetCDFFile
68        g = Geo_reference(56,1.9,1.9)
69        file_name = tempfile.mktemp(".geo_referenceTest")
70        fd = open(file_name,'w')
71        g.write_ASCII(fd)
72        fd.close()
73       
74        fd = open(file_name,'r')
75        new_g = Geo_reference(ASCIIFile=fd)
76        fd.close()
77        os.remove(file_name)
78
79        self.failUnless(g == new_g, 'test_read_write_ASCII failed') 
80   
81    def test_read_write_ASCII2(self):
82        from Scientific.IO.NetCDF import NetCDFFile
83        g = Geo_reference(56,1.9,1.9)
84        file_name = tempfile.mktemp(".geo_referenceTest")
85        fd = open(file_name,'w')
86        g.write_ASCII(fd)
87        fd.close()       
88        fd = open(file_name,'r')
89        line = fd.readline()
90        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
91        fd.close()
92        os.remove(file_name)
93
94        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
95       
96    def test_read_write_ASCII3(self):
97        from Scientific.IO.NetCDF import NetCDFFile
98        g = Geo_reference(56,1.9,1.9)
99        file_name = tempfile.mktemp(".geo_referenceTest")
100        fd = open(file_name,'w')
101        g.write_ASCII(fd)
102        fd.close()       
103        fd = open(file_name,'r')
104        line = fd.readline()
105        line = "fail !!"
106        try:
107            new_g = Geo_reference(ASCIIFile=fd, read_title=line)
108            fd.close()
109            os.remove(file_name)
110        except TitleError:
111            fd.close()
112            os.remove(file_name)
113        else:
114            self.failUnless(0 ==1,
115                        'bad text file did not raise error!')
116           
117    def test_change_points_geo_ref(self):
118        x = 433.0
119        y = 3.0
120        g = Geo_reference(56,x,y)
121        lofl = [[3.0,311.0], [677.0,6.0]]
122        new_lofl = g.change_points_geo_ref(lofl)
123        #print "lofl",lofl
124        #print "new_lofl",new_lofl
125
126        self.failUnless(type(new_lofl) == types.ListType, ' failed')
127        self.failUnless(type(new_lofl) == type(lofl), ' failed')
128        for point,new_point in map(None,lofl,new_lofl):
129            self.failUnless(point[0]-x==new_point[0], ' failed')
130            self.failUnless(point[1]-y==new_point[1], ' failed')
131         
132       
133    def test_change_points_geo_ref2(self):
134        x = 3.0
135        y = 543.0
136        g = Geo_reference(56,x,y)
137        lofl = [[3.0,388.0]]
138        new_lofl = g.change_points_geo_ref(lofl)
139        #print "lofl",lofl
140        #print "new_lofl",new_lofl
141
142        self.failUnless(type(new_lofl) == types.ListType, ' failed')
143        self.failUnless(type(new_lofl) == type(lofl), ' failed')
144        for point,new_point in map(None,lofl,new_lofl):
145            self.failUnless(point[0]-x==new_point[0], ' failed')
146            self.failUnless(point[1]-y==new_point[1], ' failed')
147       
148    def test_change_points_geo_ref3(self):
149        x = 3.0
150        y = 443.0
151        g = Geo_reference(56,x,y)
152        lofl = [3.0,345.0]
153        new_lofl = g.change_points_geo_ref(lofl)
154        #print "lofl",lofl
155        #print "new_lofl",new_lofl
156
157        self.failUnless(type(new_lofl) == types.ListType, ' failed')
158        self.failUnless(type(new_lofl) == type(lofl), ' failed')
159        for point,new_point in map(None,[lofl],new_lofl):
160            self.failUnless(point[0]-x==new_point[0], ' failed')
161            self.failUnless(point[1]-y==new_point[1], ' failed')
162       
163   
164    def test_change_points_geo_ref4(self):
165        x = 3.0
166        y = 443.0
167        g = Geo_reference(56,x,y)
168        lofl = num.array([[3.0,323.0], [6.0,645.0]])
169        new_lofl = g.change_points_geo_ref(lofl)
170        #print "4 lofl",lofl
171        #print "4 new_lofl",new_lofl
172
173        self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
174        self.failUnless(type(new_lofl) == type(lofl), ' failed')
175        lofl[:,0] -= x
176        lofl[:,1] -= y
177        assert num.allclose(lofl,new_lofl)
178       
179    def test_change_points_geo_ref5(self):
180        x = 103.0
181        y = 3.0
182        g = Geo_reference(56,x,y)
183        lofl = num.array([[3.0,323.0]])
184
185       
186        #print "5 lofl before",lofl         
187        new_lofl = g.change_points_geo_ref(lofl.copy())
188        #print "5 lofl",lofl
189        #print "5 new_lofl",new_lofl
190
191        self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
192        self.failUnless(type(new_lofl) == type(lofl), ' failed')
193
194
195        for point,new_point in map(None,lofl,new_lofl):
196            self.failUnless(point[0]-x==new_point[0], ' failed')
197            self.failUnless(point[1]-y==new_point[1], ' failed')
198       
199    def test_change_points_geo_ref6(self):
200        x = 53.0
201        y = 3.0
202        g = Geo_reference(56,x,y)
203        lofl = num.array([355.0,3.0])
204        new_lofl = g.change_points_geo_ref(lofl.copy())       
205        #print "lofl",lofl
206        #print "new_lofl",new_lofl
207
208        self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
209        self.failUnless(type(new_lofl) == type(lofl), ' failed')
210        for point,new_point in map(None,[lofl],new_lofl):
211            self.failUnless(point[0]-x==new_point[0], ' failed')
212            self.failUnless(point[1]-y==new_point[1], ' failed')
213     
214    def test_change_points_geo_ref7(self):
215        x = 23.0
216        y = 3.0
217        point_x = 9.0
218        point_y = -60.0
219        g = Geo_reference(56,x,y)
220        points_geo_ref = Geo_reference(56,point_x,point_y)
221        lofl = [[3.0,30.0], [67.0,6.0]]
222        new_lofl = g.change_points_geo_ref(lofl,points_geo_ref=points_geo_ref)
223        #print "lofl",lofl
224        #print "new_lofl",new_lofl
225
226        self.failUnless(type(new_lofl) == types.ListType, ' failed')
227        self.failUnless(type(new_lofl) == type(lofl), ' failed')
228        for point,new_point in map(None,lofl,new_lofl):
229            self.failUnless(point[0]+point_x-x==new_point[0], ' failed')
230            self.failUnless(point[1]+point_y-y==new_point[1], ' failed')
231     
232
233    def test_get_absolute(self):
234        x = 7.0
235        y = 3.0
236       
237        g = Geo_reference(56,x,y)
238        lofl = [[3.0,34.0], [64.0,6.0]]
239        new_lofl = g.get_absolute(lofl)
240        #print "lofl",lofl
241        #print "new_lofl",new_lofl
242
243        self.failUnless(type(new_lofl) == types.ListType, ' failed')
244        self.failUnless(type(new_lofl) == type(lofl), ' failed')
245        for point,new_point in map(None,lofl,new_lofl):
246            self.failUnless(point[0]+x==new_point[0], ' failed')
247            self.failUnless(point[1]+y==new_point[1], ' failed')
248
249           
250        g = Geo_reference()
251        lofl = [[3.0,34.0], [64.0,6.0]]
252        new_lofl = g.get_absolute(lofl)
253        #print "lofl",lofl
254        #print "new_lofl",new_lofl
255
256        self.failUnless(type(new_lofl) == types.ListType, ' failed')
257        self.failUnless(type(new_lofl) == type(lofl), ' failed')
258        for point,new_point in map(None,lofl,new_lofl):
259            self.failUnless(point[0]==new_point[0], ' failed')
260            self.failUnless(point[1]==new_point[1], ' failed')
261           
262    def test_is_absolute(self):
263       
264        g = Geo_reference(34,0,0)
265        points = [[3.0,34.0], [64.0,6.0]]
266
267        assert g.is_absolute()
268
269        g = Geo_reference(34,7,-6)
270        assert not g.is_absolute()       
271
272                       
273    def test___cmp__(self):
274        g = Geo_reference(56,1.9,1.9,)
275        new_g = Geo_reference(56,1.9,1.9)
276     
277        self.failUnless(g == new_g, 'test___cmp__ failed')   
278
279
280    def test_reconcile(self):
281        g1 = Geo_reference(56,2,5)
282        g2 = Geo_reference(50,4,5)
283        g3 = Geo_reference(50,66,6)       
284        g_default = Geo_reference()               
285
286
287        g2.reconcile_zones(g3)
288        assert g2.get_zone() == g3.get_zone()
289
290        g_default.reconcile_zones(g3)
291        assert g_default.get_zone() == g3.get_zone()
292
293        g_default = Geo_reference()               
294        g3.reconcile_zones(g_default)
295        assert g_default.get_zone() == g3.get_zone()               
296
297        try:
298            g1.reconcile_zones(g2)
299        except:
300            pass
301        else:
302            msg = 'Should have raised an exception'
303            raise msg
304 
305    def test_bad_ASCII_title(self):     
306 # create an text file
307        point_file = tempfile.mktemp(".xxx")
308        fd = open(point_file,'w')
309        fd.write("# hey! \n")
310        fd.close()
311       
312        fd = open(point_file,'r')
313        #
314        #new_g = Geo_reference(ASCIIFile=fd)
315        try:
316            new_g = Geo_reference(ASCIIFile=fd)
317            fd.close()
318            os.remove(point_file)
319        except TitleError:
320            fd.close()
321            os.remove(point_file)
322        else:
323            self.failUnless(0 ==1,
324                        'bad text file did not raise error!')
325            os.remove(point_file)
326
327    def test_read_write_ASCII_test_and_fail(self):
328        from Scientific.IO.NetCDF import NetCDFFile
329
330        # This is to test a fail
331        g = Geo_reference(56,1.9,1.9)
332        file_name = tempfile.mktemp(".geo_referenceTest")
333        fd = open(file_name,'w')
334        g.write_ASCII(fd)
335        fd.close()       
336        fd = open(file_name,'r')
337        line = fd.readline()
338        line = " #Geo"
339        try:
340            new_g = Geo_reference(ASCIIFile=fd, read_title=line)
341            fd.close()
342            os.remove(file_name)
343        except TitleError:
344            fd.close()
345            os.remove(file_name)
346        else:
347            self.failUnless(0 ==1,
348                        'bad text file did not raise error!')
349
350        # this tests a pass
351        g = Geo_reference(56,1.9,1.9)
352        file_name = tempfile.mktemp(".geo_referenceTest")
353        fd = open(file_name,'w')
354        g.write_ASCII(fd)
355        fd.close()
356       
357        fd = open(file_name,'r')
358        line = fd.readline()
359        line = "#geo_yeah"
360        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
361        fd.close()
362        os.remove(file_name)
363
364        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
365       
366        # this tests a pass
367        g = Geo_reference(56,1.9,1.9)
368        file_name = tempfile.mktemp(".geo_referenceTest")
369        fd = open(file_name,'w')
370        g.write_ASCII(fd)
371        fd.close()
372       
373        fd = open(file_name,'r')
374        line = fd.readline()
375        line = "#geo crap"
376        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
377        fd.close()
378        os.remove(file_name)
379
380        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
381       
382    def test_good_title(self):     
383 # create an .xxx file
384        point_file = tempfile.mktemp(".xxx")
385        fd = open(point_file,'w')
386        fd.write("#Geo crap \n 56\n ")
387        fd.close()
388       
389        fd = open(point_file,'r')
390        #
391        #new_g = Geo_reference(ASCIIFile=fd)
392        try:
393            new_g = Geo_reference(ASCIIFile=fd)
394            fd.close()
395            os.remove(point_file)
396        except ValueError:
397            fd.close()
398            os.remove(point_file)
399        else:
400            self.failUnless(0 ==1,
401                        'bad text file did not raise error!')
402            os.remove(point_file)
403
404    def test_error_message_ShapeError(self):
405       
406        new_g = Geo_reference()
407        try:
408            new_g.get_absolute((8.9, 7.8, 9.0)) 
409        except ShapeError:
410            pass
411        else:
412            self.failUnless(0 ==1,
413                        'bad shape did not raise error!')
414            os.remove(point_file)
415           
416        new_g = Geo_reference()
417        try:
418            new_g.get_absolute(((8.9, 7.8, 9.0))) 
419        except ShapeError:
420            pass
421        else:
422            self.failUnless(0 ==1,
423                        'bad shape did not raise error!')
424            os.remove(point_file)
425
426    def test_functionality_get_absolute(self):
427        x0 = 1000.0
428        y0 = 2000.0
429        geo = Geo_reference(56, x0, y0)
430
431        # iterable points (*not* num.array())
432        points = ((2,3), (3,1), (5,2))
433        abs_points = geo.get_absolute(points)
434        # check we haven't changed 'points' itself
435        self.failIf(num.alltrue(abs_points == points))
436        new_points = abs_points.copy()
437        new_points[:,0] -= x0
438        new_points[:,1] -= y0
439        self.failUnless(num.alltrue(new_points == points))
440
441        # points in num.array()
442        points = num.array(((2,3), (3,1), (5,2)), num.Float)
443        abs_points = geo.get_absolute(points)
444        # check we haven't changed 'points' itself
445        self.failIf(num.alltrue(abs_points == points))
446        new_points = abs_points.copy()
447        new_points[:,0] -= x0
448        new_points[:,1] -= y0
449        self.failUnless(num.alltrue(new_points == points))
450
451    def test_georef_types(self):
452        '''Ensure that attributes of a georeference are of correct type.
453       
454        zone            int
455        false_easting   int
456        false_northing  int
457        xllcorner       float
458        yllcorner       float
459        '''
460
461        from Scientific.IO.NetCDF import NetCDFFile
462
463        # ensure that basic instance attributes are correct
464        g = Geo_reference(56, 1.8, 1.8)
465        self.failUnless(isinstance(g.zone, int),
466                        "geo_ref .zone should be 'int' type, "
467                        "was '%s' type" % type(g.zone)) 
468        self.failUnless(isinstance(g.false_easting, int),
469                        "geo_ref .false_easting should be int type, "
470                        "was '%s' type" % type(g.false_easting)) 
471        self.failUnless(isinstance(g.false_northing, int),
472                        "geo_ref .false_northing should be int type, "
473                        "was '%s' type" % type(g.false_northing))
474        self.failUnless(isinstance(g.xllcorner, float),
475                        "geo_ref .xllcorner should be float type, "
476                        "was '%s' type" % type(g.xllcorner))
477        self.failUnless(isinstance(g.yllcorner, float),
478                        "geo_ref .yllcorner should be float type, "
479                        "was '%s' type" % type(g.yllcorner))
480
481        # now write fikle, read back and check types again
482        file_name = tempfile.mktemp(".geo_referenceTest")
483       
484        out_file = NetCDFFile(file_name, 'w')
485        g.write_NetCDF(out_file)
486        out_file.close()
487       
488        in_file = NetCDFFile(file_name, 'r')
489        new_g = Geo_reference(NetCDFObject=in_file)
490        in_file.close()
491        os.remove(file_name)
492
493        self.failUnless(isinstance(new_g.zone, int),
494                        "geo_ref .zone should be 'int' type, "
495                        "was '%s' type" % type(new_g.zone)) 
496        self.failUnless(isinstance(new_g.false_easting, int),
497                        "geo_ref .false_easting should be int type, "
498                        "was '%s' type" % type(new_g.false_easting)) 
499        self.failUnless(isinstance(new_g.false_northing, int),
500                        "geo_ref .false_northing should be int type, "
501                        "was '%s' type" % type(new_g.false_northing))
502        self.failUnless(isinstance(new_g.xllcorner, float),
503                        "geo_ref .xllcorner should be float type, "
504                        "was '%s' type" % type(new_g.xllcorner))
505        self.failUnless(isinstance(new_g.yllcorner, float),
506                        "geo_ref .yllcorner should be float type, "
507                        "was '%s' type" % type(new_g.yllcorner))
508       
509    def test_georef_types_coerceable(self):
510        '''Ensure that attributes of a georeference are of correct type.
511       
512        zone            int
513        false_easting   int
514        false_northing  int
515        xllcorner       float
516        yllcorner       float
517        '''
518
519        # now provide wrong types but coerceable
520        g = Geo_reference(56.0, '1.8', '1.8')
521        self.failUnless(isinstance(g.zone, int),
522                        "geo_ref .zone should be 'int' type, "
523                        "was '%s' type" % type(g.zone)) 
524        self.failUnless(isinstance(g.false_easting, int),
525                        "geo_ref .false_easting should be int type, "
526                        "was '%s' type" % type(g.false_easting)) 
527        self.failUnless(isinstance(g.false_northing, int),
528                        "geo_ref .false_northing should be int type, "
529                        "was '%s' type" % type(g.false_northing))
530        self.failUnless(isinstance(g.xllcorner, float),
531                        "geo_ref .xllcorner should be float type, "
532                        "was '%s' type" % type(g.xllcorner))
533        self.failUnless(isinstance(g.yllcorner, float),
534                        "geo_ref .yllcorner should be float type, "
535                        "was '%s' type" % type(g.yllcorner))
536
537       
538#-------------------------------------------------------------
539if __name__ == "__main__":
540
541    suite = unittest.makeSuite(geo_referenceTestCase,'test')
542    runner = unittest.TextTestRunner() #verbosity=2)
543    runner.run(suite)
544   
Note: See TracBrowser for help on using the repository browser.