source: inundation/geospatial_data/test_geospatial_data.py @ 2950

Last change on this file since 2950 was 2942, checked in by duncan, 19 years ago

changing errors thrown by geo_ref. adding tests to geospatial
File size: 37.1 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, TitleError
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        fileName = tempfile.mktemp(".xya")
511        file = open(fileName,"w")
512        file.write("  elevation   speed \n\
5131.0 0.0 10.0 0.0\n\
5140.0 1.0 0.0 10.0\n\
5151.0 0.0 10.4 40.0\n\
516#geo\n\
51756\n\
5181.1\n\
5191.0\n")
520
521        file.close()
522        try:
523            results = Geospatial_data(fileName, delimiter=',')
524        except IOError:
525            pass
526        else:
527            msg = 'bad xya file did not raise error!'
528            raise msg
529
530        os.remove(fileName)
531
532    def test_loadxy_bad5(self):
533        """
534        specifying wrong delimiter
535        """
536        import os
537       
538        fileName = tempfile.mktemp(".xya")
539        file = open(fileName,"w")
540        file.write("  elevation   speed \n\
5411.0 0.0 10.0 0.0\n\
5420.0 1.0 0.0 10.0\n\
5431.0 0.0 10.4 40.0\n\
544#geocrap\n\
545crap")
546        file.close()
547        try:
548#            dict = import_points_file(fileName,delimiter=' ')
549#            results = Geospatial_data()
550            results = Geospatial_data(fileName, delimiter=' ', verbose=True)
551#            results.import_points_file(fileName, delimiter=' ')
552        except IOError:
553            pass
554        else:
555            msg = 'bad xya file did not raise error!'
556            raise msg
557
558#            self.failUnless(0 ==1,
559#                        'bad xya file did not raise error!')   
560        os.remove(fileName)             
561
562    def test_loadxy_bad_no_file_xya(self):
563        import os
564       
565        fileName = tempfile.mktemp(".xya")
566        try:
567            results = Geospatial_data(fileName, delimiter=' ')
568        except IOError:
569            pass
570        else:
571            msg = 'imaginary file did not raise error!'
572            raise msg
573
574#        except IOError:
575#            pass
576#        else:
577#            self.failUnless(0 == 1,
578#                        'imaginary file did not raise error!')
579
580                       
581  ###################### .XYA ##############################
582       
583    def test_export_xya_file(self):
584#        dict = {}
585        att_dict = {}
586        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
587        att_dict['elevation'] = array([10.0, 0.0, 10.4])
588        att_dict['brightness'] = array([10.0, 0.0, 10.4])
589#        dict['attributelist'] = att_dict
590        geo_reference=Geo_reference(56,1.9,1.9)
591       
592       
593        fileName = tempfile.mktemp(".xya")
594        G = Geospatial_data(pointlist, att_dict, geo_reference)
595        G.export_points_file(fileName, False)
596
597#        dict2 = import_points_file(fileName)
598        results = Geospatial_data(file_name = fileName)
599        #print "fileName",fileName
600        os.remove(fileName)
601       
602        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
603        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
604        answer = [10.0, 0.0, 10.4]
605        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
606        #print "dict2['geo_reference']",dict2['geo_reference']
607        self.failUnless(results.get_geo_reference() == geo_reference,
608                         'test_writepts failed. Test geo_reference')
609
610    def test_export_xya_file2(self):
611        """test absolute xya file
612        """
613        att_dict = {}
614        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
615        att_dict['elevation'] = array([10.0, 0.0, 10.4])
616        att_dict['brightness'] = array([10.0, 0.0, 10.4])
617       
618        fileName = tempfile.mktemp(".xya")
619        G = Geospatial_data(pointlist, att_dict)
620        G.export_points_file(fileName)
621        results = Geospatial_data(file_name = fileName)
622#        dict2 = import_points_file(fileName)
623        os.remove(fileName)
624       
625        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
626        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
627        answer = [10.0, 0.0, 10.4]
628        assert allclose(results.get_attributes('brightness'), answer)
629
630    def test_export_xya_file3(self):
631        """test absolute xya file with geo_ref
632        """
633        att_dict = {}
634        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
635        att_dict['elevation'] = array([10.0, 0.0, 10.4])
636        att_dict['brightness'] = array([10.0, 0.0, 10.4])
637        geo_reference=Geo_reference(56,1.9,1.9)
638       
639       
640        fileName = tempfile.mktemp(".xya")
641        G = Geospatial_data(pointlist, att_dict, geo_reference)
642       
643        G.export_points_file(fileName, absolute=True)
644       
645        results = Geospatial_data(file_name = fileName)
646        os.remove(fileName)
647
648        assert allclose(results.get_data_points(),
649                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
650        assert allclose(results.get_attributes(attribute_name='elevation'),
651                         [10.0, 0.0, 10.4])
652        answer = [10.0, 0.0, 10.4]
653        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
654        self.failUnless(results.get_geo_reference() == geo_reference,
655                         'test_writepts failed. Test geo_reference')                         
656                       
657                       
658                       
659    def test_new_export_pts_file(self):
660        att_dict = {}
661        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
662        att_dict['elevation'] = array([10.1, 0.0, 10.4])
663        att_dict['brightness'] = array([10.0, 1.0, 10.4])
664       
665        fileName = tempfile.mktemp(".pts")
666       
667        G = Geospatial_data(pointlist, att_dict)
668       
669        G.export_points_file(fileName)
670
671        results = Geospatial_data(file_name = fileName)
672
673        os.remove(fileName)
674       
675        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
676        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
677        answer = [10.0, 1.0, 10.4]
678        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
679
680    def test_new_export_absolute_pts_file(self):
681        att_dict = {}
682        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
683        att_dict['elevation'] = array([10.1, 0.0, 10.4])
684        att_dict['brightness'] = array([10.0, 1.0, 10.4])
685        geo_ref = Geo_reference(50, 25, 55)
686       
687        fileName = tempfile.mktemp(".pts")
688       
689        G = Geospatial_data(pointlist, att_dict, geo_ref)
690       
691        G.export_points_file(fileName, absolute=True)
692
693        results = Geospatial_data(file_name = fileName)
694
695        os.remove(fileName)
696       
697        assert allclose(results.get_data_points(), G.get_data_points(True))
698        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
699        answer = [10.0, 1.0, 10.4]
700        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
701
702    def test_loadpts(self):
703       
704        from Scientific.IO.NetCDF import NetCDFFile
705
706        fileName = tempfile.mktemp(".pts")
707        # NetCDF file definition
708        outfile = NetCDFFile(fileName, 'w')
709       
710        # dimension definitions
711        outfile.createDimension('number_of_points', 3)   
712        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
713   
714        # variable definitions
715        outfile.createVariable('points', Float, ('number_of_points',
716                                                 'number_of_dimensions'))
717        outfile.createVariable('elevation', Float, ('number_of_points',))
718   
719        # Get handles to the variables
720        points = outfile.variables['points']
721        elevation = outfile.variables['elevation']
722 
723        points[0, :] = [1.0,0.0]
724        elevation[0] = 10.0 
725        points[1, :] = [0.0,1.0]
726        elevation[1] = 0.0 
727        points[2, :] = [1.0,0.0]
728        elevation[2] = 10.4   
729
730        outfile.close()
731       
732        results = Geospatial_data(file_name = fileName)
733        os.remove(fileName)
734        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
735        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
736        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
737       
738    def test_writepts(self):
739        att_dict = {}
740        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
741        att_dict['elevation'] = array([10.0, 0.0, 10.4])
742        att_dict['brightness'] = array([10.0, 0.0, 10.4])
743        geo_reference=Geo_reference(56,1.9,1.9)
744       
745        fileName = tempfile.mktemp(".pts")
746       
747        G = Geospatial_data(pointlist, att_dict, geo_reference)
748       
749        G.export_points_file(fileName, False)
750       
751        results = Geospatial_data(file_name = fileName)
752
753        os.remove(fileName)
754
755        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
756        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
757        answer = [10.0, 0.0, 10.4]
758        assert allclose(results.get_attributes('brightness'), answer)
759
760       
761        self.failUnless(geo_reference == geo_reference,
762                         'test_writepts failed. Test geo_reference')
763       
764 ########################## BAD .PTS ##########################         
765
766    def test_load_bad_no_file_pts(self):
767        import os
768        import tempfile
769       
770        fileName = tempfile.mktemp(".pts")
771        #print fileName
772        try:
773            results = Geospatial_data(file_name = fileName)
774#            dict = import_points_file(fileName)
775        except IOError:
776            pass
777        else:
778            msg = 'imaginary file did not raise error!'
779            raise msg
780#            self.failUnless(0 == 1,
781#                        'imaginary file did not raise error!')
782
783
784    def test_create_from_pts_file(self):
785       
786        from Scientific.IO.NetCDF import NetCDFFile
787
788#        fileName = tempfile.mktemp(".pts")
789        FN = 'test_points.pts'
790        # NetCDF file definition
791        outfile = NetCDFFile(FN, 'w')
792       
793        # dimension definitions
794        outfile.createDimension('number_of_points', 3)   
795        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
796   
797        # variable definitions
798        outfile.createVariable('points', Float, ('number_of_points',
799                                                 'number_of_dimensions'))
800        outfile.createVariable('elevation', Float, ('number_of_points',))
801   
802        # Get handles to the variables
803        points = outfile.variables['points']
804        elevation = outfile.variables['elevation']
805 
806        points[0, :] = [1.0,0.0]
807        elevation[0] = 10.0 
808        points[1, :] = [0.0,1.0]
809        elevation[1] = 0.0 
810        points[2, :] = [1.0,0.0]
811        elevation[2] = 10.4   
812
813        outfile.close()
814
815        G = Geospatial_data(file_name = FN)
816
817        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
818        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
819
820        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
821        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
822        os.remove(FN)
823
824    def test_create_from_pts_file_with_geo(self):
825        """This test reveals if Geospatial data is correctly instantiated from a pts file.
826        """
827       
828        from Scientific.IO.NetCDF import NetCDFFile
829
830        FN = 'test_points.pts'
831        # NetCDF file definition
832        outfile = NetCDFFile(FN, 'w')
833
834        # Make up an arbitrary georef
835        xll = 0.1
836        yll = 20
837        geo_reference=Geo_reference(56, xll, yll)
838        geo_reference.write_NetCDF(outfile)
839
840        # dimension definitions
841        outfile.createDimension('number_of_points', 3)   
842        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
843   
844        # variable definitions
845        outfile.createVariable('points', Float, ('number_of_points',
846                                                 'number_of_dimensions'))
847        outfile.createVariable('elevation', Float, ('number_of_points',))
848   
849        # Get handles to the variables
850        points = outfile.variables['points']
851        elevation = outfile.variables['elevation']
852
853        points[0, :] = [1.0,0.0]
854        elevation[0] = 10.0 
855        points[1, :] = [0.0,1.0]
856        elevation[1] = 0.0 
857        points[2, :] = [1.0,0.0]
858        elevation[2] = 10.4   
859
860        outfile.close()
861
862        G = Geospatial_data(file_name = FN)
863
864        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
865        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
866
867        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
868                                              [0.0+xll, 1.0+yll],
869                                              [1.0+xll, 0.0+yll]])
870       
871        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
872        os.remove(FN)
873
874       
875    def test_add_(self):
876        '''adds an xya and pts files, reads the files and adds them
877           checking results are correct
878        '''
879
880        #FIXME (Ole): I changed the georeference arbitrarily, but the test still passes.
881       
882        # create files
883        att_dict1 = {}
884        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
885        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
886        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
887        geo_reference1 = Geo_reference(56, 2.0, 1.0)
888       
889        att_dict2 = {}
890        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
891        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
892        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
893        geo_reference2 = Geo_reference(56, 1.0, 2.0) #FIXME (Ole): I changed this, why does it still pass
894        #OK - it fails now, after the fix revealed by test_create_from_pts_file_with_geo')
895
896        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
897        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
898       
899        fileName1 = tempfile.mktemp(".xya")
900        fileName2 = tempfile.mktemp(".pts")
901
902        #makes files
903        G1.export_points_file(fileName1)
904        G2.export_points_file(fileName2)
905       
906        # add files
907       
908        G3 = Geospatial_data(file_name = fileName1)
909        G4 = Geospatial_data(file_name = fileName2)
910       
911        G = G3 + G4
912
913       
914        #read results
915#        print'res', G.get_data_points()
916#        print'res1', G.get_data_points(False)
917        assert allclose(G.get_data_points(),
918                        [[ 3.0, 1.0], [ 2.0, 2.0],
919                         [ 3.0, 1.0], [ 3.0, 3.0],
920                         [ 2.0, 4.0], [ 3.0, 3.0]])
921                         
922        assert allclose(G.get_attributes(attribute_name='elevation'),
923                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
924       
925        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
926        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
927       
928        self.failUnless(G.get_geo_reference() == geo_reference1,
929                         'test_writepts failed. Test geo_reference')
930                         
931        os.remove(fileName1)
932        os.remove(fileName2)
933       
934    def test_ensure_absolute(self):
935        points = [[2.0, 0.0],[1.0, 1.0],
936                         [2.0, 0.0],[2.0, 2.0],
937                         [1.0, 3.0],[2.0, 2.0]]
938        new_points = ensure_absolute(points)
939       
940        assert allclose(new_points, points)
941
942        points = array([[2.0, 0.0],[1.0, 1.0],
943                         [2.0, 0.0],[2.0, 2.0],
944                         [1.0, 3.0],[2.0, 2.0]])
945        new_points = ensure_absolute(points)
946       
947        assert allclose(new_points, points)
948       
949        ab_points = array([[2.0, 0.0],[1.0, 1.0],
950                         [2.0, 0.0],[2.0, 2.0],
951                         [1.0, 3.0],[2.0, 2.0]])
952       
953        mesh_origin = (56, 290000, 618000) #zone, easting, northing
954
955        data_points = zeros((ab_points.shape), Float)
956        #Shift datapoints according to new origins
957        for k in range(len(ab_points)):
958            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
959            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
960        #print "data_points",data_points     
961        new_points = ensure_absolute(data_points,
962                                             geo_reference=mesh_origin)
963        #print "new_points",new_points
964        #print "ab_points",ab_points
965           
966        assert allclose(new_points, ab_points)
967
968        geo = Geo_reference(56,67,-56)
969
970        data_points = geo.change_points_geo_ref(ab_points)   
971        new_points = ensure_absolute(data_points,
972                                             geo_reference=geo)
973        #print "new_points",new_points
974        #print "ab_points",ab_points
975           
976        assert allclose(new_points, ab_points)
977
978
979        geo_reference = Geo_reference(56, 100, 200)
980        ab_points = [[1.0, 2.1], [3.0, 5.3]]
981        points = geo_reference.change_points_geo_ref(ab_points)
982        attributes = [2, 4]
983        #print "geo in points", points
984        G = Geospatial_data(points, attributes,
985                            geo_reference=geo_reference)
986         
987        new_points = ensure_absolute(G)
988        #print "new_points",new_points
989        #print "ab_points",ab_points
990           
991        assert allclose(new_points, ab_points)
992       
993    def test_isinstance(self):
994
995        import os
996       
997        fileName = tempfile.mktemp(".xya")
998        file = open(fileName,"w")
999        file.write("  elevation   speed \n\
10001.0 0.0 10.0 0.0\n\
10010.0 1.0 0.0 10.0\n\
10021.0 0.0 10.4 40.0\n\
1003#geocrap\n\
100456\n\
100556.6\n\
10063\n")
1007        file.close()
1008
1009        results = Geospatial_data(fileName)
1010       
1011        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1012        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
1013        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
1014
1015        os.remove(fileName)
1016       
1017    def test_delimiter(self):
1018       
1019        try:
1020            G = Geospatial_data(delimiter=',')
1021#            results = Geospatial_data(file_name = fileName)
1022#            dict = import_points_file(fileName)
1023        except ValueError:
1024            pass
1025        else:
1026            msg = 'Instance with No fileName but has a delimiter\
1027                  did not raise error!'
1028            raise msg
1029
1030    def test_no_constructors(self):
1031       
1032        try:
1033            G = Geospatial_data()
1034#            results = Geospatial_data(file_name = fileName)
1035#            dict = import_points_file(fileName)
1036        except ValueError:
1037            pass
1038        else:
1039            msg = 'Instance must have a filename or data points'
1040            raise msg       
1041
1042    def test_check_geo_reference(self):
1043        """
1044        checks geo reference details are OK. eg can be called '#geo reference'
1045        if not throws a clear error message
1046        """
1047        import os
1048        fileName = tempfile.mktemp(".xya")
1049        file = open(fileName,"w")
1050        file.write("  elevation  \n\
10511.0 0.0 10.0\n\
10520.0 1.0 0.0\n\
10531.0 0.0 10.4\n\
1054#ge oreference\n\
105556\n\
10561.1\n\
10571.0\n")
1058
1059        file.close()
1060        results = Geospatial_data(fileName)
1061        assert allclose(results.get_geo_reference().get_xllcorner(), 1.1)
1062        assert allclose(results.get_geo_reference().get_yllcorner(), 1.0)
1063
1064        fileName = tempfile.mktemp(".xya")
1065        file = open(fileName,"w")
1066        file.write("  elevation  \n\
10671.0 0.0 10.0\n\
10680.0 1.0 0.0\n\
10691.0 0.0 10.4\n")
1070
1071        file.close()
1072        results = Geospatial_data(fileName)
1073       
1074    def test_check_geo_reference1(self):
1075        """
1076        checks geo reference details are OK. eg can be called '#geo reference'
1077        if not throws a clear error message
1078        """
1079        import os
1080        fileName = tempfile.mktemp(".xya")
1081        file = open(fileName,"w")
1082        file.write("  elevation  \n\
10831.0 0.0 10.0\n\
10840.0 1.0 0.0\n\
10851.0 0.0 10.4\n\
1086#geo t t\n\
108756\n\
10881.1\n"
1089)
1090        file.close()
1091
1092        try:
1093            results = Geospatial_data(fileName, delimiter = " ")
1094        except IOError:
1095            pass
1096        else:
1097            msg = 'Geo reference data format is incorrect'
1098            raise msg       
1099
1100
1101
1102
1103if __name__ == "__main__":
1104    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_check_geo_reference1')
1105    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1106    runner = unittest.TextTestRunner()
1107    runner.run(suite)
1108
1109   
Note: See TracBrowser for help on using the repository browser.