source: inundation/geospatial_data/test_geospatial_data.py @ 2570

Last change on this file since 2570 was 2570, checked in by nick, 17 years ago

Added functionality to geospatial_data.py including import and exporting pts and xya files
File size: 22.5 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import os
6from Numeric import zeros, array, allclose, concatenate
7from math import sqrt, pi
8import tempfile
9
10from geospatial_data import *
11
12from coordinate_transforms.geo_reference import Geo_reference
13
14class Test_Geospatial_data(unittest.TestCase):
15    def setUp(self):
16        pass
17
18    def tearDown(self):
19        pass
20
21
22    def test_0(self):
23        """Basic points
24        """
25        from coordinate_transforms.geo_reference import Geo_reference
26       
27        points = [[1.0, 2.1], [3.0, 5.3]]
28        G = Geospatial_data(points)
29
30        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
31
32        #Check defaults
33        assert G.attributes is None
34       
35        assert G.geo_reference.zone == Geo_reference().zone
36        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
37        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
38       
39
40    def test_1(self):
41        points = [[1.0, 2.1], [3.0, 5.3]]
42        attributes = [2, 4]
43        G = Geospatial_data(points, attributes)       
44
45        assert G.attributes.keys()[0] == 'attribute'
46        assert allclose(G.attributes.values()[0], [2, 4])
47       
48
49    def test_2(self):
50        from coordinate_transforms.geo_reference import Geo_reference
51        points = [[1.0, 2.1], [3.0, 5.3]]
52        attributes = [2, 4]
53        G = Geospatial_data(points, attributes,
54                            geo_reference = Geo_reference(56, 100, 200))
55
56        assert G.geo_reference.zone == 56
57        assert G.geo_reference.xllcorner == 100
58        assert G.geo_reference.yllcorner == 200
59
60
61    def test_get_attributes_1(self):
62        from coordinate_transforms.geo_reference import Geo_reference
63        points = [[1.0, 2.1], [3.0, 5.3]]
64        attributes = [2, 4]
65        G = Geospatial_data(points, attributes,
66                            geo_reference = Geo_reference(56, 100, 200))
67
68
69        P = G.get_data_points(absolute = False)
70        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
71
72        P = G.get_data_points(absolute = True)
73        assert allclose(P, [[101.0, 202.1], [103.0, 205.3]])       
74
75        V = G.get_attributes() #Simply get them
76        assert allclose(V, [2, 4])
77
78        V = G.get_attributes('attribute') #Get by name
79        assert allclose(V, [2, 4])
80
81    def test_get_attributes_2(self):
82        """Multiple attributes
83        """
84       
85        from coordinate_transforms.geo_reference import Geo_reference
86        points = [[1.0, 2.1], [3.0, 5.3]]
87        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
88        G = Geospatial_data(points, attributes,
89                            geo_reference = Geo_reference(56, 100, 200),
90                            default_attribute_name = 'a1')
91
92
93        P = G.get_data_points(absolute = False)
94        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
95       
96        V = G.get_attributes() #Get default attribute
97        assert allclose(V, [2, 4])
98
99        V = G.get_attributes('a0') #Get by name
100        assert allclose(V, [0, 0])
101
102        V = G.get_attributes('a1') #Get by name
103        assert allclose(V, [2, 4])
104
105        V = G.get_attributes('a2') #Get by name
106        assert allclose(V, [79.4, -7])
107
108        try:
109            V = G.get_attributes('hdnoatedu') #Invalid
110        except AssertionError:
111            pass
112        else:
113            raise 'Should have raised exception' 
114
115
116    def test_conversions_to_points_dict(self):
117        """test conversions to points_dict
118        """
119       
120        from coordinate_transforms.geo_reference import Geo_reference
121        points = [[1.0, 2.1], [3.0, 5.3]]
122        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
123        G = Geospatial_data(points, attributes,
124                            geo_reference = Geo_reference(56, 100, 200),
125                            default_attribute_name = 'a1')
126
127
128        points_dict = geospatial_data2points_dictionary(G)
129
130        assert points_dict.has_key('pointlist')
131        assert points_dict.has_key('attributelist')       
132        assert points_dict.has_key('geo_reference')
133
134        assert allclose( points_dict['pointlist'], points )
135
136        A = points_dict['attributelist']
137        assert A.has_key('a0')
138        assert A.has_key('a1')
139        assert A.has_key('a2')       
140
141        assert allclose( A['a0'], [0, 0] )
142        assert allclose( A['a1'], [2, 4] )       
143        assert allclose( A['a2'], [79.4, -7] )
144
145
146        geo = points_dict['geo_reference']
147        assert geo is G.geo_reference
148
149
150    def test_conversions_from_points_dict(self):
151        """test conversions from points_dict
152        """
153
154        from coordinate_transforms.geo_reference import Geo_reference
155       
156        points = [[1.0, 2.1], [3.0, 5.3]]
157        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
158
159        points_dict = {}
160        points_dict['pointlist'] = points
161        points_dict['attributelist'] = attributes
162        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
163       
164
165        G = points_dictionary2geospatial_data(points_dict)
166
167        P = G.get_data_points(absolute = False)
168        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
169       
170        #V = G.get_attribute_values() #Get default attribute
171        #assert allclose(V, [2, 4])
172
173        V = G.get_attributes('a0') #Get by name
174        assert allclose(V, [0, 0])
175
176        V = G.get_attributes('a1') #Get by name
177        assert allclose(V, [2, 4])
178
179        V = G.get_attributes('a2') #Get by name
180        assert allclose(V, [79.4, -7])
181
182    def test_add(self):
183        """ test the addition of two geospatical objects
184            no geo_reference see next test
185        """
186        points = [[1.0, 2.1], [3.0, 5.3]]
187        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
188        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
189        G1 = Geospatial_data(points, attributes)       
190        G2 = Geospatial_data(points, attributes1) 
191       
192#        print 'G1 depth=', G1.get_attributes('depth')
193#        print 'G1 attrib keys=', G1.attributes.keys()
194        g3 = geospatial_data2points_dictionary(G1)
195#        print 'g3=', g3
196       
197        G = G1 + G2
198
199#        g = geospatial_data2points_dictionary(G)
200#        print 'G=', g
201#        print 'G points =', G.get_data_points()
202#        print 'G attrib keys =', G.attributes.keys()
203#        print 'G test =', G.get_attributes('elevation')
204        assert G.attributes.has_key('depth')
205        assert G.attributes.has_key('elevation')
206        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
207        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
208        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
209                                              [1.0, 2.1], [3.0, 5.3]])
210       
211    def test_add1 (self):
212        """
213        difference in Geo_reference resolved
214        """
215        points = [[1.0, 2.1], [3.0, 5.3]]
216        points1 = [[5.0, 6.1], [6.0, 3.3]]
217        attributes = [2, 4]
218        attributes1 = [5, 76]
219        geo_ref = Geo_reference(55, 1.0, 2.0)
220        geo_ref1 = Geo_reference(zone = 55, xllcorner = 0.1,
221                                yllcorner = 3.0, datum = 'wgs84',
222                                projection = 'UTM', units = 'm')
223                               
224        G1 = Geospatial_data(points, attributes, geo_ref)
225        G2 = Geospatial_data(points1, attributes1, geo_ref1)
226
227        G = G1 + G2
228       
229        assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
230        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
231#        print 'G data points=', G.get_data_points()
232        assert allclose(G.get_data_points(), [[2.0, 4.1], [4.0, 7.3], [5.1, 9.1], [6.1, 6.3]])                             
233
234
235
236
237    def test_create_from_xya_file(self):
238        """Check that object can be created from a points file (.pts and .xya)
239        """
240
241        from load_mesh.loadASCII import export_points_file
242       
243
244        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
245        attributes = [2, 4, 5, 76]
246
247        # Use old pointsdict format
248        pointsdict = {'pointlist': points,
249                      'attributelist': {'att1': attributes,
250                                        'att2': array(attributes) + 1}} 
251     
252        # Create points as an xya file
253        FN = 'test_points.xya'
254        export_points_file(FN, pointsdict)
255
256
257        #Create object from file
258        G = Geospatial_data(file_name = FN)
259       
260        assert allclose(G.get_data_points(), points)
261        assert allclose(G.get_attributes('att1'), attributes)
262        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
263       
264        os.remove(FN)
265       
266    def test_loadxya(self):
267        """
268        comma delimited
269        """
270        fileName = tempfile.mktemp(".xya")
271        file = open(fileName,"w")
272        file.write("elevation  , speed \n\
2731.0, 0.0, 10.0, 0.0\n\
2740.0, 1.0, 0.0, 10.0\n\
2751.0, 0.0, 10.4, 40.0\n")
276        file.close()
277        #print fileName
278        dict = import_points_file(fileName,delimiter = ',')
279        os.remove(fileName)
280        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
281        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
282        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
283
284    def test_loadxya2(self):
285        """
286        space delimited
287        """
288        import os
289       
290        fileName = tempfile.mktemp(".xya")
291        file = open(fileName,"w")
292        file.write("  elevation   speed \n\
2931.0 0.0 10.0 0.0\n\
2940.0 1.0 0.0 10.0\n\
2951.0 0.0 10.4 40.0\n")
296        file.close()
297        #print fileName
298#        dict = import_points_file(fileName,delimiter = ' ')
299        dict = import_points_file(fileName,delimiter = ' ')
300        os.remove(fileName)
301        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
302        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
303        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
304     
305    def test_loadxya3(self):
306        """
307        space delimited
308        """
309        import os
310       
311        fileName = tempfile.mktemp(".xya")
312        file = open(fileName,"w")
313        file.write("  elevation   speed \n\
3141.0 0.0 10.0 0.0\n\
3150.0 1.0 0.0 10.0\n\
3161.0 0.0 10.4 40.0\n\
317#geocrap\n\
31856\n\
31956.6\n\
3203\n")
321        file.close()
322        #print fileName
323        dict = import_points_file(fileName,delimiter = ' ')
324        os.remove(fileName)
325        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
326        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
327        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
328
329    def test_loadxy_bad_no_file_xya(self):
330        import os
331       
332        fileName = tempfile.mktemp(".xya")
333        #print fileName
334        try:
335            dict = import_points_file(fileName,delimiter = ' ')
336        except IOError:
337            pass
338        else:
339            self.failUnless(0 ==1,
340                        'imaginary file did not raise error!')
341 
342    def Xtest_read_write_points_file_bad(self):
343        #not used yet as i'm uncertain if "tri_dict" is nesscessary to geospatial_data
344        dict = self.tri_dict.copy()
345        fileName = tempfile.mktemp(".xxx")
346        try:
347            export_points_file(fileName,dict)
348        except IOError:
349            pass
350        else:
351            self.failUnless(0 ==1,
352                        'bad points file extension did not raise error!')
353
354    def Xtest_read_write_points_file_bad2(self):
355        dict = {}
356        att_dict = {}
357        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
358        att_dict['elevation'] = array([10.0, 0.0, 10.4])
359        att_dict['brightness'] = array([10.0, 0.0, 10.4])
360        dict['attributelist'] = att_dict
361        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
362        try:
363            export_points_file("_???/yeah.xya",dict)
364        except IOError:
365            pass
366        else:
367            self.failUnless(0 ==1,
368                        'bad points file extension did not raise error!')
369                   
370    def test_loadxy_bad(self):
371        import os
372       
373        fileName = tempfile.mktemp(".xya")
374        file = open(fileName,"w")
375        file.write("  elevation   \n\
3761.0 0.0 10.0 0.0\n\
3770.0 1.0 0.0 10.0\n\
3781.0 0.0 10.4 40.0\n")
379        file.close()
380        #print fileName
381        try:
382            dict = import_points_file(fileName,delimiter = ' ')
383        except IOError:
384            pass
385        else:
386            self.failUnless(0 ==1,
387                        'bad xya file did not raise error!')
388        os.remove(fileName)
389       
390    def test_loadxy_bad2(self):
391        import os
392       
393        fileName = tempfile.mktemp(".xya")
394        file = open(fileName,"w")
395        file.write("elevation\n\
3961.0 0.0 10.0 \n\
3970.0 1.0\n\
3981.0 \n")
399        file.close()
400        #print fileName
401        try:
402            dict = import_points_file(fileName,delimiter = ' ')
403        except IOError:
404            pass
405        else:
406            self.failUnless(0 ==1,
407                        'bad xya file did not raise error!')   
408        os.remove(fileName)
409   
410    def test_loadxy_bad3(self):
411        """
412        specifying wrong delimiter
413        """
414        import os
415       
416        fileName = tempfile.mktemp(".xya")
417        file = open(fileName,"w")
418        file.write("  elevation  , speed \n\
4191.0, 0.0, 10.0, 0.0\n\
4200.0, 1.0, 0.0, 10.0\n\
4211.0, 0.0, 10.4, 40.0\n")
422        file.close()
423        try:
424            dict = import_points_file(fileName,delimiter = ' ')
425        except IOError:
426            pass
427        else:
428            self.failUnless(0 ==1,
429                        'bad xya file did not raise error!')   
430        os.remove(fileName)
431     
432    def test_loadxy_bad4(self):
433        """
434        specifying wrong delimiter
435        """
436        import os
437       
438        fileName = tempfile.mktemp(".xya")
439        file = open(fileName,"w")
440        file.write("  elevation   speed \n\
4411.0 0.0 10.0 0.0\n\
4420.0 1.0 0.0 10.0\n\
4431.0 0.0 10.4 40.0\n\
444yeah")
445        file.close()
446        try:
447            dict = import_points_file(fileName,delimiter = ' ')
448        except IOError:
449            pass
450        else:
451            self.failUnless(0 ==1,
452                        'bad xya file did not raise error!')   
453        os.remove(fileName)
454 
455    def test_loadxy_bad4(self):
456        """
457        specifying wrong delimiter
458        """
459        import os
460       
461        fileName = tempfile.mktemp(".xya")
462        file = open(fileName,"w")
463        file.write("  elevation   speed \n\
4641.0 0.0 10.0 0.0\n\
4650.0 1.0 0.0 10.0\n\
4661.0 0.0 10.4 40.0\n\
467#geocrap")
468        file.close()
469        try:
470            dict = import_points_file(fileName,delimiter = ' ')
471        except IOError:
472            pass
473        else:
474            self.failUnless(0 ==1,
475                        'bad xya file did not raise error!')   
476        os.remove(fileName)
477
478    def test_loadxy_bad5(self):
479        """
480        specifying wrong delimiter
481        """
482        import os
483       
484        fileName = tempfile.mktemp(".xya")
485        file = open(fileName,"w")
486        file.write("  elevation   speed \n\
4871.0 0.0 10.0 0.0\n\
4880.0 1.0 0.0 10.0\n\
4891.0 0.0 10.4 40.0\n\
490#geocrap\n\
491crap")
492        file.close()
493        try:
494            dict = import_points_file(fileName,delimiter = ' ')
495        except IOError:
496            pass
497        else:
498            self.failUnless(0 ==1,
499                        'bad xya file did not raise error!')   
500        os.remove(fileName)             
501
502    def test_loadxy_bad_no_file_xya(self):
503        import os
504       
505        fileName = tempfile.mktemp(".xya")
506        #print fileName
507        try:
508            dict = import_points_file(fileName,delimiter = ' ')
509        except IOError:
510            pass
511        else:
512            self.failUnless(0 ==1,
513                        'imaginary file did not raise error!')
514                       
515  ###################### .XYA ##############################
516       
517    def test_export_xya_file(self):
518        dict = {}
519        att_dict = {}
520        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
521        att_dict['elevation'] = array([10.0, 0.0, 10.4])
522        att_dict['brightness'] = array([10.0, 0.0, 10.4])
523        dict['attributelist'] = att_dict
524        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
525       
526       
527        fileName = tempfile.mktemp(".xya")
528        export_points_file(fileName, dict)
529        dict2 = import_points_file(fileName)
530        #print "fileName",fileName
531        os.remove(fileName)
532        #print "dict2",dict2
533       
534        assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
535        assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
536        answer = [10.0, 0.0, 10.4]
537        assert allclose(dict2['attributelist']['brightness'], answer)
538        #print "dict2['geo_reference']",dict2['geo_reference']
539        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
540                         'test_writepts failed. Test geo_reference')
541
542    def test_export_xya_file2(self):
543        dict = {}
544        att_dict = {}
545        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
546        att_dict['elevation'] = array([10.0, 0.0, 10.4])
547        att_dict['brightness'] = array([10.0, 0.0, 10.4])
548        dict['attributelist'] = att_dict
549       
550       
551        fileName = tempfile.mktemp(".xya")
552        export_points_file(fileName, dict)
553        dict2 = import_points_file(fileName)
554        #print "fileName",fileName
555        os.remove(fileName)
556        #print "dict2",dict2
557       
558        assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
559        assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
560        answer = [10.0, 0.0, 10.4]
561        assert allclose(dict2['attributelist']['brightness'], answer)
562
563
564    def test_loadpts(self):
565       
566        from Scientific.IO.NetCDF import NetCDFFile
567
568        fileName = tempfile.mktemp(".pts")
569        # NetCDF file definition
570        outfile = NetCDFFile(fileName, 'w')
571       
572        # dimension definitions
573        outfile.createDimension('number_of_points', 3)   
574        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
575   
576        # variable definitions
577        outfile.createVariable('points', Float, ('number_of_points',
578                                                 'number_of_dimensions'))
579        outfile.createVariable('elevation', Float, ('number_of_points',))
580   
581        # Get handles to the variables
582        points = outfile.variables['points']
583        elevation = outfile.variables['elevation']
584 
585        points[0, :] = [1.0,0.0]
586        elevation[0] = 10.0 
587        points[1, :] = [0.0,1.0]
588        elevation[1] = 0.0 
589        points[2, :] = [1.0,0.0]
590        elevation[2] = 10.4   
591
592        outfile.close()
593       
594        dict = import_points_file(fileName)
595        os.remove(fileName)
596        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
597        #print "dict['pointlist']",dict['pointlist']
598        #print "answer",answer
599        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
600        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
601       
602    def test_writepts(self):
603        dict = {}
604        att_dict = {}
605        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
606        att_dict['elevation'] = array([10.0, 0.0, 10.4])
607        att_dict['brightness'] = array([10.0, 0.0, 10.4])
608        dict['attributelist'] = att_dict
609        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
610       
611       
612        fileName = tempfile.mktemp(".pts")
613        export_points_file(fileName, dict)
614        dict2 = import_points_file(fileName)
615        #print "fileName",fileName
616        os.remove(fileName)
617        #print "dict2",dict2
618       
619        assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
620        assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
621        answer = [10.0, 0.0, 10.4]
622        assert allclose(dict2['attributelist']['brightness'], answer)
623
624        #print "dict['geo_reference'] ",dict['geo_reference']
625        #print "dict2['geo_reference']",dict2['geo_reference']
626       
627        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
628                         'test_writepts failed. Test geo_reference')
629       
630 ########################## BAD .PTS ##########################         
631
632    def test_load_bad_no_file_pts(self):
633        import os
634        import tempfile
635       
636        fileName = tempfile.mktemp(".pts")
637        #print fileName
638        try:
639            dict = import_points_file(fileName)
640        except IOError:
641            pass
642        else:
643            self.failUnless(0 ==1,
644                        'imaginary file did not raise error!')
645
646
647    def test_create_from_pts_file(self):
648       
649        from Scientific.IO.NetCDF import NetCDFFile
650
651
652#        fileName = tempfile.mktemp(".pts")
653        FN = 'test_points.pts'
654        # NetCDF file definition
655        outfile = NetCDFFile(FN, 'w')
656       
657        # dimension definitions
658        outfile.createDimension('number_of_points', 3)   
659        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
660   
661        # variable definitions
662        outfile.createVariable('points', Float, ('number_of_points',
663                                                 'number_of_dimensions'))
664        outfile.createVariable('elevation', Float, ('number_of_points',))
665   
666        # Get handles to the variables
667        points = outfile.variables['points']
668        elevation = outfile.variables['elevation']
669 
670        points[0, :] = [1.0,0.0]
671        elevation[0] = 10.0 
672        points[1, :] = [0.0,1.0]
673        elevation[1] = 0.0 
674        points[2, :] = [1.0,0.0]
675        elevation[2] = 10.4   
676
677        outfile.close()
678
679        G = Geospatial_data(file_name = FN)
680
681        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
682        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
683        os.remove(FN)
684
685if __name__ == "__main__":
686    suite = unittest.makeSuite(Test_Geospatial_data,'test')
687    runner = unittest.TextTestRunner()
688    runner.run(suite)
689
690   
Note: See TracBrowser for help on using the repository browser.