source: inundation/geospatial_data/test_geospatial_data.py @ 2928

Last change on this file since 2928 was 2928, checked in by nick, 19 years ago

geospatial_data object updated so not to use dictionaries

File size: 35.9 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#        g3 = geospatial_data2points_dictionary(G1)
193#        print 'g3=', g3
194       
195        G = G1 + G2
196
197        assert G.attributes.has_key('depth')
198        assert G.attributes.has_key('elevation')
199        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
200        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
201        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
202                                              [1.0, 2.1], [3.0, 5.3]])
203       
204    def test_add_with_geo (self):
205        """
206        Difference in Geo_reference resolved
207        """
208        points1 = [[1.0, 2.1], [3.0, 5.3]]
209        points2 = [[5.0, 6.1], [6.0, 3.3]]
210        attributes1 = [2, 4]
211        attributes2 = [5, 76]
212        geo_ref1= Geo_reference(55, 1.0, 2.0)
213        geo_ref2 = Geo_reference(zone=55,
214                                 xllcorner=0.1,
215                                 yllcorner=3.0,
216                                 datum='wgs84',
217                                 projection='UTM',
218                                 units='m')
219                               
220        G1 = Geospatial_data(points1, attributes1, geo_ref1)
221        G2 = Geospatial_data(points2, attributes2, geo_ref2)
222
223        #Check that absolute values are as expected
224        P1 = G1.get_data_points(absolute=True)
225        assert allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
226
227        P2 = G2.get_data_points(absolute=True)
228        assert allclose(P2, [[5.1, 9.1], [6.1, 6.3]])       
229       
230        G = G1 + G2
231       
232        assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
233        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
234
235        P = G.get_data_points(absolute=True)
236
237        P_relative = G.get_data_points(absolute=False)
238       
239        assert allclose(P_relative, P - [0.1, 2.0])
240
241        assert allclose(P, concatenate( (P1,P2) ))
242        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
243                            [5.1, 9.1], [6.1, 6.3]])
244       
245
246
247       
248
249    def test_add_with_geo_absolute (self):
250        """
251        Difference in Geo_reference resolved
252        """
253        points1 = array([[2.0, 4.1], [4.0, 7.3]])
254        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
255        attributes1 = [2, 4]
256        attributes2 = [5, 76]
257        geo_ref1= Geo_reference(55, 1.0, 2.0)
258        geo_ref2 = Geo_reference(55, 2.0, 3.0)
259
260       
261                               
262        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
263                             attributes1, geo_ref1)
264       
265        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
266                             attributes2, geo_ref2)
267
268        #Check that absolute values are as expected
269        P1 = G1.get_data_points(absolute=True)
270        assert allclose(P1, points1)
271
272        P1 = G1.get_data_points(absolute=False)
273        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
274
275        P2 = G2.get_data_points(absolute=True)
276        assert allclose(P2, points2)
277
278        P2 = G2.get_data_points(absolute=False)
279        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
280       
281        G = G1 + G2
282       
283        assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
284        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
285
286        P = G.get_data_points(absolute=True)
287
288        P_relative = G.get_data_points(absolute=False)
289       
290        assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
291
292        assert allclose(P, concatenate( (points1,points2) ))
293                           
294       
295
296
297    def test_create_from_xya_file(self):
298        """Check that object can be created from a points file (.pts and .xya)
299        """
300
301        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
302        attributes = [2, 4, 5, 76]
303        '''
304        # Use old pointsdict format
305        pointsdict = {'pointlist': points,
306                      'attributelist': {'att1': attributes,
307                                        'att2': array(attributes) + 1}}
308        '''
309        att_dict = {'att1': attributes,
310                    'att2': array(attributes) +1}
311                   
312        # Create points as an xya file
313        FN = 'test_points.xya'
314        G1 = Geospatial_data(points, att_dict)
315        G1.export_points_file(FN)
316#        G1.export_points_file(ofile)
317
318        #Create object from file
319        G = Geospatial_data(file_name = FN)
320       
321        assert allclose(G.get_data_points(), points)
322        assert allclose(G.get_attributes('att1'), attributes)
323        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
324       
325        os.remove(FN)
326
327    def test_create_from_xya_file1(self):
328        """
329        Check that object can be created from an Absolute xya file
330        """
331
332        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
333        attributes = [2, 4, 5, 76]
334
335        att_dict = {'att1': attributes,
336                    'att2': array(attributes) +1}
337
338        geo_ref = Geo_reference(56, 10, 5)
339                   
340        # Create points as an xya file
341        FN = 'test_points.xya'
342        G1 = Geospatial_data(points, att_dict, geo_ref)
343
344        G1.export_points_file(FN, absolute = True)
345
346        #Create object from file
347        G = Geospatial_data(file_name = FN)
348       
349        assert allclose(G.get_data_points(absolute = True), 
350                        G1.get_data_points(absolute = True))
351        assert allclose(G.get_attributes('att1'), attributes)
352        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
353       
354        os.remove(FN)
355       
356    def test_loadxya(self):
357        """
358        comma delimited
359        """
360        fileName = tempfile.mktemp(".xya")
361        file = open(fileName,"w")
362        file.write("elevation  , speed \n\
3631.0, 0.0, 10.0, 0.0\n\
3640.0, 1.0, 0.0, 10.0\n\
3651.0, 0.0, 10.4, 40.0\n")
366        file.close()
367        results = Geospatial_data(fileName, delimiter = ',')
368        os.remove(fileName)
369#        print 'data', results.get_data_points()
370        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
371        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
372        assert allclose(results.get_attributes(attribute_name = 'speed'), [0.0, 10.0, 40.0])
373
374    def test_loadxya2(self):
375        """
376        space delimited
377        """
378        import os
379       
380        fileName = tempfile.mktemp(".xya")
381        file = open(fileName,"w")
382        file.write("  elevation   speed \n\
3831.0 0.0 10.0 0.0\n\
3840.0 1.0 0.0 10.0\n\
3851.0 0.0 10.4 40.0\n")
386        file.close()
387
388        results = Geospatial_data(fileName, delimiter = ' ')
389
390        os.remove(fileName)
391
392        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
393        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
394        assert allclose(results.get_attributes(attribute_name = 'speed'), [0.0, 10.0, 40.0])
395     
396    def test_loadxya3(self):
397        """
398        space delimited
399        """
400        import os
401       
402        fileName = tempfile.mktemp(".xya")
403        file = open(fileName,"w")
404        file.write("  elevation   speed \n\
4051.0 0.0 10.0 0.0\n\
4060.0 1.0 0.0 10.0\n\
4071.0 0.0 10.4 40.0\n\
408#geocrap\n\
40956\n\
41056.6\n\
4113\n")
412        file.close()
413
414        results = Geospatial_data(fileName, delimiter = ' ')
415
416        os.remove(fileName)
417        assert allclose(results.get_data_points(absolute = False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
418        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
419        assert allclose(results.get_attributes(attribute_name = 'speed'), [0.0, 10.0, 40.0])
420
421    def test_read_write_points_file_bad2(self):
422        att_dict = {}
423        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
424        att_dict['elevation'] = array([10.0, 0.0, 10.4])
425        att_dict['brightness'] = array([10.0, 0.0, 10.4])
426        geo_reference = Geo_reference(56,1.9,1.9)
427       
428        G = Geospatial_data(pointlist, att_dict, geo_reference)
429       
430        try:
431            G.export_points_file("_???/yeah.xya")
432           
433        except IOError:
434            pass
435        else:
436            msg = 'bad points file extension did not raise error!'
437            raise msg
438#            self.failUnless(0 == 1,
439#                        'bad points file extension did not raise error!')
440                   
441    def test_loadxy_bad(self):
442        import os
443       
444        fileName = tempfile.mktemp(".xya")
445        file = open(fileName,"w")
446        file.write("  elevation   \n\
4471.0 0.0 10.0 0.0\n\
4480.0 1.0 0.0 10.0\n\
4491.0 0.0 10.4 40.0\n")
450        file.close()
451        #print fileName
452        try:
453            results = Geospatial_data(fileName, delimiter = ' ')
454        except IOError:
455            pass
456        else:
457            msg = 'bad xya file did not raise error!'
458            raise msg
459#            self.failUnless(0 == 1,
460#                        'bad xya file did not raise error!')
461        os.remove(fileName)
462       
463    def test_loadxy_bad2(self):
464        import os
465       
466        fileName = tempfile.mktemp(".xya")
467        file = open(fileName,"w")
468        file.write("elevation\n\
4691.0 0.0 10.0 \n\
4700.0 1.0\n\
4711.0 \n")
472        file.close()
473        #print fileName
474        try:
475            results = Geospatial_data(fileName, delimiter = ' ')
476        except IOError:
477            pass
478        else:
479            msg = 'bad xya file did not raise error!'
480            raise msg
481        os.remove(fileName)
482   
483    def test_loadxy_bad3(self):
484        """
485        specifying wrong delimiter
486        """
487        import os
488       
489        fileName = tempfile.mktemp(".xya")
490        file = open(fileName,"w")
491        file.write("  elevation  , speed \n\
4921.0, 0.0, 10.0, 0.0\n\
4930.0, 1.0, 0.0, 10.0\n\
4941.0, 0.0, 10.4, 40.0\n")
495        file.close()
496        try:
497            results = Geospatial_data(fileName, delimiter = ' ')
498        except IOError:
499            pass
500        else:
501            msg = 'bad xya file did not raise error!'
502            raise msg
503        os.remove(fileName)
504     
505    def test_loadxy_bad4(self):
506        """
507        specifying wrong delimiter
508        """
509        import os
510       
511        fileName = tempfile.mktemp(".xya")
512        file = open(fileName,"w")
513        file.write("  elevation   speed \n\
5141.0 0.0 10.0 0.0\n\
5150.0 1.0 0.0 10.0\n\
5161.0 0.0 10.4 40.0\n\
517#geocrap")
518        file.close()
519        try:
520#            dict = import_points_file(fileName,delimiter = ' ')
521#            results = Geospatial_data()
522            results = Geospatial_data(fileName, delimiter = ',')
523#            results.import_points_file(fileName, delimiter = ' ')
524        except IOError:
525            pass
526        else:
527            msg = 'bad xya file did not raise error!'
528            raise msg
529
530#            self.failUnless(0 ==1,
531#                        'bad xya file did not raise error!')   
532        os.remove(fileName)
533
534    def test_loadxy_bad5(self):
535        """
536        specifying wrong delimiter
537        """
538        import os
539       
540        fileName = tempfile.mktemp(".xya")
541        file = open(fileName,"w")
542        file.write("  elevation   speed \n\
5431.0 0.0 10.0 0.0\n\
5440.0 1.0 0.0 10.0\n\
5451.0 0.0 10.4 40.0\n\
546#geocrap\n\
547crap")
548        file.close()
549        try:
550#            dict = import_points_file(fileName,delimiter = ' ')
551#            results = Geospatial_data()
552            results = Geospatial_data(fileName, delimiter = ' ')
553#            results.import_points_file(fileName, delimiter = ' ')
554        except IOError:
555            pass
556        else:
557            msg = 'bad xya file did not raise error!'
558            raise msg
559
560#            self.failUnless(0 ==1,
561#                        'bad xya file did not raise error!')   
562        os.remove(fileName)             
563
564    def test_loadxy_bad_no_file_xya(self):
565        import os
566       
567        fileName = tempfile.mktemp(".xya")
568        try:
569            results = Geospatial_data(fileName, delimiter = ' ')
570        except IOError:
571            pass
572        else:
573            msg = 'imaginary file did not raise error!'
574            raise msg
575
576#        except IOError:
577#            pass
578#        else:
579#            self.failUnless(0 == 1,
580#                        'imaginary file did not raise error!')
581
582                       
583  ###################### .XYA ##############################
584       
585    def test_export_xya_file(self):
586#        dict = {}
587        att_dict = {}
588        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
589        att_dict['elevation'] = array([10.0, 0.0, 10.4])
590        att_dict['brightness'] = array([10.0, 0.0, 10.4])
591#        dict['attributelist'] = att_dict
592        geo_reference = Geo_reference(56,1.9,1.9)
593       
594       
595        fileName = tempfile.mktemp(".xya")
596        G = Geospatial_data(pointlist, att_dict, geo_reference)
597        G.export_points_file(fileName)
598
599#        dict2 = import_points_file(fileName)
600        results = Geospatial_data(file_name = fileName)
601        #print "fileName",fileName
602        os.remove(fileName)
603       
604        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
605        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
606        answer = [10.0, 0.0, 10.4]
607        assert allclose(results.get_attributes(attribute_name = 'brightness'), answer)
608        #print "dict2['geo_reference']",dict2['geo_reference']
609        self.failUnless(results.get_geo_reference() == geo_reference,
610                         'test_writepts failed. Test geo_reference')
611
612    def test_export_xya_file2(self):
613        """test absolute xya file
614        """
615        att_dict = {}
616        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
617        att_dict['elevation'] = array([10.0, 0.0, 10.4])
618        att_dict['brightness'] = array([10.0, 0.0, 10.4])
619       
620        fileName = tempfile.mktemp(".xya")
621        G = Geospatial_data(pointlist, att_dict)
622        G.export_points_file(fileName)
623        results = Geospatial_data(file_name = fileName)
624#        dict2 = import_points_file(fileName)
625        os.remove(fileName)
626       
627        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
628        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
629        answer = [10.0, 0.0, 10.4]
630        assert allclose(results.get_attributes('brightness'), answer)
631
632    def test_export_xya_file3(self):
633        """test absolute xya file with geo_ref
634        """
635        att_dict = {}
636        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
637        att_dict['elevation'] = array([10.0, 0.0, 10.4])
638        att_dict['brightness'] = array([10.0, 0.0, 10.4])
639        geo_reference = Geo_reference(56,1.9,1.9)
640       
641       
642        fileName = tempfile.mktemp(".xya")
643        G = Geospatial_data(pointlist, att_dict, geo_reference)
644       
645        G.export_points_file(fileName, absolute=True)
646       
647        results = Geospatial_data(file_name = fileName)
648        os.remove(fileName)
649
650        assert allclose(results.get_data_points(),
651                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
652        assert allclose(results.get_attributes(attribute_name = 'elevation'),
653                         [10.0, 0.0, 10.4])
654        answer = [10.0, 0.0, 10.4]
655        assert allclose(results.get_attributes(attribute_name = 'brightness'), answer)
656        self.failUnless(results.get_geo_reference() == geo_reference,
657                         'test_writepts failed. Test geo_reference')                         
658                       
659                       
660                       
661    def test_new_export_pts_file(self):
662        att_dict = {}
663        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
664        att_dict['elevation'] = array([10.1, 0.0, 10.4])
665        att_dict['brightness'] = array([10.0, 1.0, 10.4])
666       
667        fileName = tempfile.mktemp(".pts")
668       
669        G = Geospatial_data(pointlist, att_dict)
670       
671        G.export_points_file(fileName)
672
673        results = Geospatial_data(file_name = fileName)
674
675        os.remove(fileName)
676       
677        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
678        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.1, 0.0, 10.4])
679        answer = [10.0, 1.0, 10.4]
680        assert allclose(results.get_attributes(attribute_name = 'brightness'), answer)
681
682    def test_new_export_absolute_pts_file(self):
683        att_dict = {}
684        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
685        att_dict['elevation'] = array([10.1, 0.0, 10.4])
686        att_dict['brightness'] = array([10.0, 1.0, 10.4])
687        geo_ref = Geo_reference(50, 25, 55)
688       
689        fileName = tempfile.mktemp(".pts")
690       
691        G = Geospatial_data(pointlist, att_dict, geo_ref)
692       
693        G.export_points_file(fileName, absolute=True)
694
695        results = Geospatial_data(file_name = fileName)
696
697        os.remove(fileName)
698       
699        assert allclose(results.get_data_points(), G.get_data_points(True))
700        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.1, 0.0, 10.4])
701        answer = [10.0, 1.0, 10.4]
702        assert allclose(results.get_attributes(attribute_name = 'brightness'), answer)
703
704    def test_loadpts(self):
705       
706        from Scientific.IO.NetCDF import NetCDFFile
707
708        fileName = tempfile.mktemp(".pts")
709        # NetCDF file definition
710        outfile = NetCDFFile(fileName, 'w')
711       
712        # dimension definitions
713        outfile.createDimension('number_of_points', 3)   
714        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
715   
716        # variable definitions
717        outfile.createVariable('points', Float, ('number_of_points',
718                                                 'number_of_dimensions'))
719        outfile.createVariable('elevation', Float, ('number_of_points',))
720   
721        # Get handles to the variables
722        points = outfile.variables['points']
723        elevation = outfile.variables['elevation']
724 
725        points[0, :] = [1.0,0.0]
726        elevation[0] = 10.0 
727        points[1, :] = [0.0,1.0]
728        elevation[1] = 0.0 
729        points[2, :] = [1.0,0.0]
730        elevation[2] = 10.4   
731
732        outfile.close()
733       
734        results = Geospatial_data(file_name = fileName)
735        os.remove(fileName)
736        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
737        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
738        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
739       
740    def test_writepts(self):
741        att_dict = {}
742        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
743        att_dict['elevation'] = array([10.0, 0.0, 10.4])
744        att_dict['brightness'] = array([10.0, 0.0, 10.4])
745        geo_reference = Geo_reference(56,1.9,1.9)
746       
747        fileName = tempfile.mktemp(".pts")
748       
749        G = Geospatial_data(pointlist, att_dict, geo_reference)
750       
751        G.export_points_file(fileName)
752       
753        results = Geospatial_data(file_name = fileName)
754
755        os.remove(fileName)
756
757        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
758        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
759        answer = [10.0, 0.0, 10.4]
760        assert allclose(results.get_attributes('brightness'), answer)
761
762       
763        self.failUnless(geo_reference == geo_reference,
764                         'test_writepts failed. Test geo_reference')
765       
766 ########################## BAD .PTS ##########################         
767
768    def test_load_bad_no_file_pts(self):
769        import os
770        import tempfile
771       
772        fileName = tempfile.mktemp(".pts")
773        #print fileName
774        try:
775            results = Geospatial_data(file_name = fileName)
776#            dict = import_points_file(fileName)
777        except IOError:
778            pass
779        else:
780            msg = 'imaginary file did not raise error!'
781            raise msg
782#            self.failUnless(0 == 1,
783#                        'imaginary file did not raise error!')
784
785
786    def test_create_from_pts_file(self):
787       
788        from Scientific.IO.NetCDF import NetCDFFile
789
790#        fileName = tempfile.mktemp(".pts")
791        FN = 'test_points.pts'
792        # NetCDF file definition
793        outfile = NetCDFFile(FN, 'w')
794       
795        # dimension definitions
796        outfile.createDimension('number_of_points', 3)   
797        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
798   
799        # variable definitions
800        outfile.createVariable('points', Float, ('number_of_points',
801                                                 'number_of_dimensions'))
802        outfile.createVariable('elevation', Float, ('number_of_points',))
803   
804        # Get handles to the variables
805        points = outfile.variables['points']
806        elevation = outfile.variables['elevation']
807 
808        points[0, :] = [1.0,0.0]
809        elevation[0] = 10.0 
810        points[1, :] = [0.0,1.0]
811        elevation[1] = 0.0 
812        points[2, :] = [1.0,0.0]
813        elevation[2] = 10.4   
814
815        outfile.close()
816
817        G = Geospatial_data(file_name = FN)
818
819        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
820        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
821
822        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
823        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
824        os.remove(FN)
825
826    def test_create_from_pts_file_with_geo(self):
827        """This test reveals if Geospatial data is correctly instantiated from a pts file.
828        """
829       
830        from Scientific.IO.NetCDF import NetCDFFile
831
832        FN = 'test_points.pts'
833        # NetCDF file definition
834        outfile = NetCDFFile(FN, 'w')
835
836        # Make up an arbitrary georef
837        xll = 0.1
838        yll = 20
839        geo_reference = Geo_reference(56, xll, yll)
840        geo_reference.write_NetCDF(outfile)
841
842        # dimension definitions
843        outfile.createDimension('number_of_points', 3)   
844        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
845   
846        # variable definitions
847        outfile.createVariable('points', Float, ('number_of_points',
848                                                 'number_of_dimensions'))
849        outfile.createVariable('elevation', Float, ('number_of_points',))
850   
851        # Get handles to the variables
852        points = outfile.variables['points']
853        elevation = outfile.variables['elevation']
854
855        points[0, :] = [1.0,0.0]
856        elevation[0] = 10.0 
857        points[1, :] = [0.0,1.0]
858        elevation[1] = 0.0 
859        points[2, :] = [1.0,0.0]
860        elevation[2] = 10.4   
861
862        outfile.close()
863
864        G = Geospatial_data(file_name = FN)
865
866        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
867        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
868
869        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
870                                              [0.0+xll, 1.0+yll],
871                                              [1.0+xll, 0.0+yll]])
872       
873        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
874        os.remove(FN)
875
876       
877    def test_add_points_files(self):
878        '''adds an xya and pts files and checks resulting file is correct
879        '''
880
881        #FIXME (Ole): I changed the georeference arbitrarily, but the test still passes.
882       
883        # create files
884#        dict1 = {}
885        att_dict1 = {}
886        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
887        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
888        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
889#        dict1['attributelist'] = att_dict1
890        geo_reference1 = Geo_reference(56, 2.0, 1.0)
891       
892#        dict2 = {}
893        att_dict2 = {}
894        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
895        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
896        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
897#        dict2['attributelist'] = att_dict2
898        geo_reference2 = Geo_reference(56, 1.0, 2.0) #FIXME (Ole): I changed this, why does it still pass
899        #OK - it fails now, after the fix revealed by test_create_from_pts_file_with_geo')
900
901        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
902        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
903       
904        fileName1 = tempfile.mktemp(".xya")
905        fileName2 = tempfile.mktemp(".pts")
906
907        #makes files
908        G1.export_points_file(fileName1)
909        G2.export_points_file(fileName2)
910       
911        # add files
912       
913        G3 = Geospatial_data(file_name = fileName1)
914        G4 = Geospatial_data(file_name = fileName2)
915       
916        G = G3 + G4
917
918       
919        #read results
920        assert allclose(G.get_data_points(absolute = False),
921                        [[2.0, 0.0],[1.0, 1.0],
922                         [2.0, 0.0],[2.0, 2.0],
923                         [1.0, 3.0],[2.0, 2.0]])
924        assert allclose(G.get_attributes(attribute_name = 'elevation'),
925                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
926       
927        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
928        assert allclose(G.get_attributes(attribute_name = 'brightness'), answer)
929       
930        self.failUnless(G.get_geo_reference() == geo_reference1,
931                         'test_writepts failed. Test geo_reference')
932                         
933        os.remove(fileName1)
934        os.remove(fileName2)
935       
936    def test_ensure_absolute(self):
937        points = [[2.0, 0.0],[1.0, 1.0],
938                         [2.0, 0.0],[2.0, 2.0],
939                         [1.0, 3.0],[2.0, 2.0]]
940        new_points = ensure_absolute(points)
941       
942        assert allclose(new_points, points)
943
944        points = array([[2.0, 0.0],[1.0, 1.0],
945                         [2.0, 0.0],[2.0, 2.0],
946                         [1.0, 3.0],[2.0, 2.0]])
947        new_points = ensure_absolute(points)
948       
949        assert allclose(new_points, points)
950       
951        ab_points = array([[2.0, 0.0],[1.0, 1.0],
952                         [2.0, 0.0],[2.0, 2.0],
953                         [1.0, 3.0],[2.0, 2.0]])
954       
955        mesh_origin = (56, 290000, 618000) #zone, easting, northing
956
957        data_points = zeros((ab_points.shape), Float)
958        #Shift datapoints according to new origins
959        for k in range(len(ab_points)):
960            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
961            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
962        #print "data_points",data_points     
963        new_points = ensure_absolute(data_points,
964                                             geo_reference = mesh_origin)
965        #print "new_points",new_points
966        #print "ab_points",ab_points
967           
968        assert allclose(new_points, ab_points)
969
970        geo = Geo_reference(56,67,-56)
971
972        data_points = geo.change_points_geo_ref(ab_points)   
973        new_points = ensure_absolute(data_points,
974                                             geo_reference = geo)
975        #print "new_points",new_points
976        #print "ab_points",ab_points
977           
978        assert allclose(new_points, ab_points)
979
980
981        geo_reference = Geo_reference(56, 100, 200)
982        ab_points = [[1.0, 2.1], [3.0, 5.3]]
983        points = geo_reference.change_points_geo_ref(ab_points)
984        attributes = [2, 4]
985        #print "geo in points", points
986        G = Geospatial_data(points, attributes,
987                            geo_reference = geo_reference)
988         
989        new_points = ensure_absolute(G)
990        #print "new_points",new_points
991        #print "ab_points",ab_points
992           
993        assert allclose(new_points, ab_points)
994       
995    def test_isinstance(self):
996
997        import os
998       
999        fileName = tempfile.mktemp(".xya")
1000        file = open(fileName,"w")
1001        file.write("  elevation   speed \n\
10021.0 0.0 10.0 0.0\n\
10030.0 1.0 0.0 10.0\n\
10041.0 0.0 10.4 40.0\n\
1005#geocrap\n\
100656\n\
100756.6\n\
10083\n")
1009        file.close()
1010
1011        results = Geospatial_data(fileName)
1012       
1013        assert allclose(results.get_data_points(absolute = False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1014        assert allclose(results.get_attributes(attribute_name = 'elevation'), [10.0, 0.0, 10.4])
1015        assert allclose(results.get_attributes(attribute_name = 'speed'), [0.0, 10.0, 40.0])
1016
1017        os.remove(fileName)
1018       
1019    def test_delimiter(self):
1020       
1021        try:
1022            G = Geospatial_data(delimiter = ',')
1023#            results = Geospatial_data(file_name = fileName)
1024#            dict = import_points_file(fileName)
1025        except ValueError:
1026            pass
1027        else:
1028            msg = 'Instance with No fileName but has a delimiter\
1029                  did not raise error!'
1030            raise msg
1031
1032    def test_no_constructors(self):
1033       
1034        try:
1035            G = Geospatial_data()
1036#            results = Geospatial_data(file_name = fileName)
1037#            dict = import_points_file(fileName)
1038        except ValueError:
1039            pass
1040        else:
1041            msg = 'Instance must have a filename or data points'
1042            raise msg       
1043
1044
1045
1046if __name__ == "__main__":
1047    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_ensure_absolute')
1048    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1049    runner = unittest.TextTestRunner()
1050    runner.run(suite)
1051
1052   
Note: See TracBrowser for help on using the repository browser.