source: branches/numpy/anuga/coordinate_transforms/test_geo_reference.py @ 6410

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

numpy changes.

File size: 17.3 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 numpy 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(isinstance(new_lofl, num.ndarray), ' 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(isinstance(new_lofl, num.ndarray), ' 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(isinstance(new_lofl, num.ndarray), ' 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    def test_get_absolute_list(self):
233        # test with supplied offsets
234        x = 7.0
235        y = 3.0
236       
237        g = Geo_reference(56, x, y)
238        points = [[3.0,34.0], [64.0,6.0]]
239        new_points = g.get_absolute(points)
240
241        self.failUnless(type(new_points) == types.ListType, 'failed')
242        self.failUnless(type(new_points) == type(points), 'failed')
243        for point, new_point in map(None, points, new_points):
244            self.failUnless(point[0]+x == new_point[0], 'failed')
245            self.failUnless(point[1]+y == new_point[1], 'failed')
246
247        # test with no supplied offsets
248        g = Geo_reference()
249        points = [[3.0,34.0], [64.0,6.0]]
250        new_points = g.get_absolute(points)
251
252        self.failUnless(type(new_points) == types.ListType, 'failed')
253        self.failUnless(type(new_points) == type(points), 'failed')
254        for point, new_point in map(None, points, new_points):
255            self.failUnless(point[0] == new_point[0], 'failed')
256            self.failUnless(point[1] == new_point[1], 'failed')
257           
258        # test that calling get_absolute twice does the right thing
259        # first call
260        dx = 10.0
261        dy = 10.0
262        g = Geo_reference(56, dx, dy)
263        points = [[3.0,34.0], [64.0,6.0]]
264        expected_new_points = [[3.0+dx,34.0+dy], [64.0+dx,6.0+dy]]
265        new_points = g.get_absolute(points)
266
267        self.failUnless(type(new_points) == types.ListType, 'failed')
268        self.failUnless(type(new_points) == type(points), 'failed')
269        for point, new_point in map(None, expected_new_points, new_points):
270            self.failUnless(point[0] == new_point[0], 'failed')
271            self.failUnless(point[1] == new_point[1], 'failed')
272
273        # and repeat from 'new_points = g.get_absolute(points)' above
274        # to see if second call with same input gives same results.
275        new_points = g.get_absolute(points)
276
277        self.failUnless(type(new_points) == types.ListType, 'failed')
278        self.failUnless(type(new_points) == type(points), 'failed')
279        for point, new_point in map(None, expected_new_points, new_points):
280            self.failUnless(point[0] == new_point[0], 'failed')
281            self.failUnless(point[1] == new_point[1], 'failed')
282
283    def test_get_absolute_array(self):
284        '''Same test as test_get_absolute_list(), but with numeric arrays.'''
285
286        # test with supplied offsets
287        x = 7.0
288        y = 3.0
289       
290        g = Geo_reference(56, x, y)
291        points = num.array([[3.0,34.0], [64.0,6.0]])
292        new_points = g.get_absolute(points)
293
294        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
295        self.failUnless(type(new_points) == type(points), 'failed')
296        self.failUnless(num.alltrue(points == new_points), 'failed')
297
298        # test with no supplied offsets
299        g = Geo_reference()
300        points = num.array([[3.0,34.0], [64.0,6.0]])
301        new_points = g.get_absolute(points)
302
303        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
304        self.failUnless(type(new_points) == type(points), 'failed')
305        self.failUnless(num.alltrue(points == new_points), 'failed')
306
307        # test that calling get_absolute twice does the right thing
308        # first call
309        dx = 10.0
310        dy = 10.0
311        g = Geo_reference(56, dx, dy)
312        points = num.array([[3.0,34.0], [64.0,6.0]])
313        expected_new_points = num.array([[3.0+dx,34.0+dy], [64.0+dx,6.0+dy]])
314        new_points = g.get_absolute(points)
315
316        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
317        self.failUnless(type(new_points) == type(points), 'failed')
318        msg = ('First call of .get_absolute() returned %s\nexpected %s'
319               % (str(new_points), str(expected_new_points)))
320        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
321
322        # and repeat from 'new_points = g.get_absolute(points)' above
323        # to see if second call with same input gives same results.
324        new_points = g.get_absolute(points)
325
326        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
327        self.failUnless(type(new_points) == type(points), 'failed')
328        msg = ('Second call of .get_absolute() returned\n%s\nexpected\n%s'
329               % (str(new_points), str(expected_new_points)))
330        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
331
332        # and repeat again to see if *third* call with same input
333        # gives same results.
334        new_points = g.get_absolute(points)
335
336        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
337        self.failUnless(type(new_points) == type(points), 'failed')
338        msg = ('Second call of .get_absolute() returned %s\nexpected %s'
339               % (str(new_points), str(expected_new_points)))
340        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
341
342    def test_is_absolute(self):
343       
344        g = Geo_reference(34,0,0)
345        points = [[3.0,34.0], [64.0,6.0]]
346
347        assert g.is_absolute()
348
349        g = Geo_reference(34,7,-6)
350        assert not g.is_absolute()       
351
352                       
353    def test___cmp__(self):
354        g = Geo_reference(56,1.9,1.9,)
355        new_g = Geo_reference(56,1.9,1.9)
356     
357        self.failUnless(g == new_g, 'test___cmp__ failed')   
358
359
360    def test_reconcile(self):
361        g1 = Geo_reference(56,2,5)
362        g2 = Geo_reference(50,4,5)
363        g3 = Geo_reference(50,66,6)       
364        g_default = Geo_reference()               
365
366
367        g2.reconcile_zones(g3)
368        assert g2.get_zone() == g3.get_zone()
369
370        g_default.reconcile_zones(g3)
371        assert g_default.get_zone() == g3.get_zone()
372
373        g_default = Geo_reference()               
374        g3.reconcile_zones(g_default)
375        assert g_default.get_zone() == g3.get_zone()               
376
377        try:
378            g1.reconcile_zones(g2)
379        except:
380            pass
381        else:
382            msg = 'Should have raised an exception'
383            raise msg
384 
385    def test_bad_ASCII_title(self):     
386 # create an text file
387        point_file = tempfile.mktemp(".xxx")
388        fd = open(point_file,'w')
389        fd.write("# hey! \n")
390        fd.close()
391       
392        fd = open(point_file,'r')
393        #
394        #new_g = Geo_reference(ASCIIFile=fd)
395        try:
396            new_g = Geo_reference(ASCIIFile=fd)
397            fd.close()
398            os.remove(point_file)
399        except TitleError:
400            fd.close()
401            os.remove(point_file)
402        else:
403            self.failUnless(0 ==1,
404                        'bad text file did not raise error!')
405            os.remove(point_file)
406
407    def test_read_write_ASCII_test_and_fail(self):
408        from Scientific.IO.NetCDF import NetCDFFile
409
410        # This is to test a fail
411        g = Geo_reference(56,1.9,1.9)
412        file_name = tempfile.mktemp(".geo_referenceTest")
413        fd = open(file_name,'w')
414        g.write_ASCII(fd)
415        fd.close()       
416        fd = open(file_name,'r')
417        line = fd.readline()
418        line = " #Geo"
419        try:
420            new_g = Geo_reference(ASCIIFile=fd, read_title=line)
421            fd.close()
422            os.remove(file_name)
423        except TitleError:
424            fd.close()
425            os.remove(file_name)
426        else:
427            self.failUnless(0 ==1,
428                        'bad text file did not raise error!')
429
430        # this tests a pass
431        g = Geo_reference(56,1.9,1.9)
432        file_name = tempfile.mktemp(".geo_referenceTest")
433        fd = open(file_name,'w')
434        g.write_ASCII(fd)
435        fd.close()
436       
437        fd = open(file_name,'r')
438        line = fd.readline()
439        line = "#geo_yeah"
440        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
441        fd.close()
442        os.remove(file_name)
443
444        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
445       
446        # this tests a pass
447        g = Geo_reference(56,1.9,1.9)
448        file_name = tempfile.mktemp(".geo_referenceTest")
449        fd = open(file_name,'w')
450        g.write_ASCII(fd)
451        fd.close()
452       
453        fd = open(file_name,'r')
454        line = fd.readline()
455        line = "#geo crap"
456        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
457        fd.close()
458        os.remove(file_name)
459
460        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
461       
462    def test_good_title(self):     
463 # create an .xxx file
464        point_file = tempfile.mktemp(".xxx")
465        fd = open(point_file,'w')
466        fd.write("#Geo crap \n 56\n ")
467        fd.close()
468       
469        fd = open(point_file,'r')
470        #
471        #new_g = Geo_reference(ASCIIFile=fd)
472        try:
473            new_g = Geo_reference(ASCIIFile=fd)
474            fd.close()
475            os.remove(point_file)
476        except ValueError:
477            fd.close()
478            os.remove(point_file)
479        else:
480            self.failUnless(0 ==1,
481                        'bad text file did not raise error!')
482            os.remove(point_file)
483
484    def test_error_message_ShapeError(self):
485       
486        new_g = Geo_reference()
487        try:
488            new_g.get_absolute((8.9, 7.8, 9.0)) 
489        except ShapeError:
490            pass
491        else:
492            self.failUnless(0 ==1,
493                        'bad shape did not raise error!')
494            os.remove(point_file)
495           
496        new_g = Geo_reference()
497        try:
498            new_g.get_absolute(((8.9, 7.8, 9.0))) 
499        except ShapeError:
500            pass
501        else:
502            self.failUnless(0 ==1,
503                        'bad shape did not raise error!')
504            os.remove(point_file)
505       
506#-------------------------------------------------------------
507
508if __name__ == "__main__":
509    suite = unittest.makeSuite(geo_referenceTestCase,'test')
510    runner = unittest.TextTestRunner() #verbosity=2)
511    runner.run(suite)
512   
Note: See TracBrowser for help on using the repository browser.