source: inundation/geospatial_data/test_geospatial_data.py @ 3149

Last change on this file since 3149 was 3149, checked in by duncan, 18 years ago

made geospatial.set_geo_reference more flexible.

File size: 40.4 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    def test_set_geo_reference(self):
116        points_ab = [[12.5,34.7],[-4.5,-60.0]]
117        x_p = -10
118        y_p = -40
119        geo_ref = Geo_reference(56, x_p, y_p)
120        points_rel = geo_ref.change_points_geo_ref(points_ab)
121        spatial = Geospatial_data(points_rel)
122        # since the geo_ref wasn't set
123        assert not allclose( points_ab, spatial.get_data_points(absolute=True))
124       
125        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
126        assert allclose( points_ab, spatial.get_data_points(absolute=True))
127       
128        x_p = 10
129        y_p = 400
130        new_geo_ref = Geo_reference(56, x_p, y_p)
131        spatial.set_geo_reference(new_geo_ref)
132        assert allclose( points_ab, spatial.get_data_points(absolute=True))
133       
134       
135       
136    def test_conversions_to_points_dict(self):
137        """test conversions to points_dict
138        """
139       
140        from coordinate_transforms.geo_reference import Geo_reference
141        points = [[1.0, 2.1], [3.0, 5.3]]
142        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
143        G = Geospatial_data(points, attributes,
144                            geo_reference=Geo_reference(56, 100, 200),
145                            default_attribute_name='a1')
146
147
148        points_dict = geospatial_data2points_dictionary(G)
149
150        assert points_dict.has_key('pointlist')
151        assert points_dict.has_key('attributelist')       
152        assert points_dict.has_key('geo_reference')
153
154        assert allclose( points_dict['pointlist'], points )
155
156        A = points_dict['attributelist']
157        assert A.has_key('a0')
158        assert A.has_key('a1')
159        assert A.has_key('a2')       
160
161        assert allclose( A['a0'], [0, 0] )
162        assert allclose( A['a1'], [2, 4] )       
163        assert allclose( A['a2'], [79.4, -7] )
164
165
166        geo = points_dict['geo_reference']
167        assert geo is G.geo_reference
168
169
170    def test_conversions_from_points_dict(self):
171        """test conversions from points_dict
172        """
173
174        from coordinate_transforms.geo_reference import Geo_reference
175       
176        points = [[1.0, 2.1], [3.0, 5.3]]
177        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
178
179        points_dict = {}
180        points_dict['pointlist'] = points
181        points_dict['attributelist'] = attributes
182        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
183       
184
185        G = points_dictionary2geospatial_data(points_dict)
186
187        P = G.get_data_points(absolute=False)
188        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
189       
190        #V = G.get_attribute_values() #Get default attribute
191        #assert allclose(V, [2, 4])
192
193        V = G.get_attributes('a0') #Get by name
194        assert allclose(V, [0, 0])
195
196        V = G.get_attributes('a1') #Get by name
197        assert allclose(V, [2, 4])
198
199        V = G.get_attributes('a2') #Get by name
200        assert allclose(V, [79.4, -7])
201
202    def test_add(self):
203        """ test the addition of two geospatical objects
204            no geo_reference see next test
205        """
206        points = [[1.0, 2.1], [3.0, 5.3]]
207        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
208        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
209        G1 = Geospatial_data(points, attributes)       
210        G2 = Geospatial_data(points, attributes1) 
211       
212#        g3 = geospatial_data2points_dictionary(G1)
213#        print 'g3=', g3
214       
215        G = G1 + G2
216
217        assert G.attributes.has_key('depth')
218        assert G.attributes.has_key('elevation')
219        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
220        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
221        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
222                                              [1.0, 2.1], [3.0, 5.3]])
223       
224    def test_add_with_geo (self):
225        """
226        Difference in Geo_reference resolved
227        """
228        points1 = [[1.0, 2.1], [3.0, 5.3]]
229        points2 = [[5.0, 6.1], [6.0, 3.3]]
230        attributes1 = [2, 4]
231        attributes2 = [5, 76]
232        geo_ref1= Geo_reference(55, 1.0, 2.0)
233        geo_ref2 = Geo_reference(zone=55,
234                                 xllcorner=0.1,
235                                 yllcorner=3.0,
236                                 datum='wgs84',
237                                 projection='UTM',
238                                 units='m')
239                               
240        G1 = Geospatial_data(points1, attributes1, geo_ref1)
241        G2 = Geospatial_data(points2, attributes2, geo_ref2)
242
243        #Check that absolute values are as expected
244        P1 = G1.get_data_points(absolute=True)
245        assert allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
246
247        P2 = G2.get_data_points(absolute=True)
248        assert allclose(P2, [[5.1, 9.1], [6.1, 6.3]])       
249       
250        G = G1 + G2
251       
252        assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
253        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
254
255        P = G.get_data_points(absolute=True)
256
257        P_relative = G.get_data_points(absolute=False)
258       
259        assert allclose(P_relative, P - [0.1, 2.0])
260
261        assert allclose(P, concatenate( (P1,P2) ))
262        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
263                            [5.1, 9.1], [6.1, 6.3]])
264       
265
266
267       
268
269    def test_add_with_geo_absolute (self):
270        """
271        Difference in Geo_reference resolved
272        """
273        points1 = array([[2.0, 4.1], [4.0, 7.3]])
274        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
275        attributes1 = [2, 4]
276        attributes2 = [5, 76]
277        geo_ref1= Geo_reference(55, 1.0, 2.0)
278        geo_ref2 = Geo_reference(55, 2.0, 3.0)
279
280       
281                               
282        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
283                             attributes1, geo_ref1)
284       
285        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
286                             attributes2, geo_ref2)
287
288        #Check that absolute values are as expected
289        P1 = G1.get_data_points(absolute=True)
290        assert allclose(P1, points1)
291
292        P1 = G1.get_data_points(absolute=False)
293        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
294
295        P2 = G2.get_data_points(absolute=True)
296        assert allclose(P2, points2)
297
298        P2 = G2.get_data_points(absolute=False)
299        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
300       
301        G = G1 + G2
302       
303        assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
304        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
305
306        P = G.get_data_points(absolute=True)
307
308        P_relative = G.get_data_points(absolute=False)
309       
310        assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
311
312        assert allclose(P, concatenate( (points1,points2) ))
313                           
314       
315
316
317    def test_create_from_xya_file(self):
318        """Check that object can be created from a points file (.pts and .xya)
319        """
320
321        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
322        attributes = [2, 4, 5, 76]
323        '''
324        # Use old pointsdict format
325        pointsdict = {'pointlist': points,
326                      'attributelist': {'att1': attributes,
327                                        'att2': array(attributes) + 1}}
328        '''
329        att_dict = {'att1': attributes,
330                    'att2': array(attributes) +1}
331                   
332        # Create points as an xya file
333        FN = 'test_points.xya'
334        G1 = Geospatial_data(points, att_dict)
335        G1.export_points_file(FN)
336#        G1.export_points_file(ofile)
337
338        #Create object from file
339        G = Geospatial_data(file_name = FN)
340       
341        assert allclose(G.get_data_points(), points)
342        assert allclose(G.get_attributes('att1'), attributes)
343        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
344       
345        os.remove(FN)
346
347    def test_create_from_xya_file1(self):
348        """
349        Check that object can be created from an Absolute xya file
350        """
351
352        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
353        attributes = [2, 4, 5, 76]
354
355        att_dict = {'att1': attributes,
356                    'att2': array(attributes) +1}
357
358        geo_ref = Geo_reference(56, 10, 5)
359                   
360        # Create points as an xya file
361        FN = 'test_points.xya'
362        G1 = Geospatial_data(points, att_dict, geo_ref)
363
364        G1.export_points_file(FN, absolute=True)
365
366        #Create object from file
367        G = Geospatial_data(file_name = FN)
368       
369        assert allclose(G.get_data_points(absolute=True), 
370                        G1.get_data_points(absolute=True))
371        assert allclose(G.get_attributes('att1'), attributes)
372        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
373       
374        os.remove(FN)
375       
376    def test_loadxya(self):
377        """
378        comma delimited
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        results = Geospatial_data(fileName, delimiter=',')
388        os.remove(fileName)
389#        print 'data', results.get_data_points()
390        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
391        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
392        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
393
394    def test_loadxya2(self):
395        """
396        space delimited
397        """
398        import os
399       
400        fileName = tempfile.mktemp(".xya")
401        file = open(fileName,"w")
402        file.write("  elevation   speed \n\
4031.0 0.0 10.0 0.0\n\
4040.0 1.0 0.0 10.0\n\
4051.0 0.0 10.4 40.0\n")
406        file.close()
407
408        results = Geospatial_data(fileName, delimiter=' ')
409
410        os.remove(fileName)
411
412        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
413        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
414        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
415     
416    def test_loadxya3(self):
417        """
418        space delimited
419        """
420        import os
421       
422        fileName = tempfile.mktemp(".xya")
423        file = open(fileName,"w")
424        file.write("  elevation   speed \n\
4251.0 0.0 10.0 0.0\n\
4260.0 1.0 0.0 10.0\n\
4271.0 0.0 10.4 40.0\n\
428#geocrap\n\
42956\n\
43056.6\n\
4313\n")
432        file.close()
433
434        results = Geospatial_data(fileName, delimiter=' ')
435
436        os.remove(fileName)
437        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
438        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
439        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
440
441    def test_read_write_points_file_bad2(self):
442        att_dict = {}
443        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
444        att_dict['elevation'] = array([10.0, 0.0, 10.4])
445        att_dict['brightness'] = array([10.0, 0.0, 10.4])
446        geo_reference=Geo_reference(56,1.9,1.9)
447       
448        G = Geospatial_data(pointlist, att_dict, geo_reference)
449       
450        try:
451            G.export_points_file("_???/yeah.xya")
452           
453        except IOError:
454            pass
455        else:
456            msg = 'bad points file extension did not raise error!'
457            raise msg
458#            self.failUnless(0 == 1,
459#                        'bad points file extension did not raise error!')
460                   
461    def test_loadxy_bad(self):
462        import os
463       
464        fileName = tempfile.mktemp(".xya")
465        file = open(fileName,"w")
466        file.write("  elevation   \n\
4671.0 0.0 10.0 0.0\n\
4680.0 1.0 0.0 10.0\n\
4691.0 0.0 10.4 40.0\n")
470        file.close()
471        #print fileName
472        try:
473            results = Geospatial_data(fileName, delimiter=' ')
474        except IOError:
475            pass
476        else:
477            msg = 'bad xya file did not raise error!'
478            raise msg
479#            self.failUnless(0 == 1,
480#                        'bad xya file did not raise error!')
481        os.remove(fileName)
482       
483    def test_loadxy_bad2(self):
484        import os
485       
486        fileName = tempfile.mktemp(".xya")
487        file = open(fileName,"w")
488        file.write("elevation\n\
4891.0 0.0 10.0 \n\
4900.0 1.0\n\
4911.0 \n")
492        file.close()
493        #print fileName
494        try:
495            results = Geospatial_data(fileName, delimiter=' ')
496        except IOError:
497            pass
498        else:
499            msg = 'bad xya file did not raise error!'
500            raise msg
501        os.remove(fileName)
502   
503    def test_loadxy_bad3(self):
504        """
505        specifying wrong delimiter
506        """
507        import os
508       
509        fileName = tempfile.mktemp(".xya")
510        file = open(fileName,"w")
511        file.write("  elevation  , speed \n\
5121.0, 0.0, 10.0, 0.0\n\
5130.0, 1.0, 0.0, 10.0\n\
5141.0, 0.0, 10.4, 40.0\n")
515        file.close()
516        try:
517            results = Geospatial_data(fileName, delimiter=' ')
518        except IOError:
519            pass
520        else:
521            msg = 'bad xya file did not raise error!'
522            raise msg
523        os.remove(fileName)
524     
525    def test_loadxy_bad4(self):
526        """
527         specifying wrong delimiter
528        """
529        import os
530        fileName = tempfile.mktemp(".xya")
531        file = open(fileName,"w")
532        file.write("  elevation   speed \n\
5331.0 0.0 10.0 0.0\n\
5340.0 1.0 0.0 10.0\n\
5351.0 0.0 10.4 40.0\n\
536#geocrap\n\
53756\n\
53856.6\n\
5393\n"
540)
541        file.close()
542        try:
543            results = Geospatial_data(fileName, delimiter=',')
544        except IOError:
545            pass
546        else:
547            msg = 'bad xya file did not raise error!'
548            raise msg
549
550        os.remove(fileName)
551
552    def test_loadxy_bad5(self):
553        """
554        specifying wrong delimiter
555        """
556        import os
557       
558        fileName = tempfile.mktemp(".xya")
559        file = open(fileName,"w")
560        file.write("  elevation   speed \n\
5611.0 0.0 10.0 0.0\n\
5620.0 1.0 0.0 10.0\n\
5631.0 0.0 10.4 40.0\n\
564#geocrap\n\
565crap")
566        file.close()
567        try:
568#            dict = import_points_file(fileName,delimiter=' ')
569#            results = Geospatial_data()
570            results = Geospatial_data(fileName, delimiter=' ', verbose=True)
571#            results.import_points_file(fileName, delimiter=' ')
572        except IOError:
573            pass
574        else:
575            msg = 'bad xya file did not raise error!'
576            raise msg
577
578#            self.failUnless(0 ==1,
579#                        'bad xya file did not raise error!')   
580        os.remove(fileName)             
581
582    def test_loadxy_bad_no_file_xya(self):
583        import os
584       
585        fileName = tempfile.mktemp(".xya")
586        try:
587            results = Geospatial_data(fileName, delimiter=' ')
588        except IOError:
589            pass
590        else:
591            msg = 'imaginary file did not raise error!'
592            raise msg
593
594#        except IOError:
595#            pass
596#        else:
597#            self.failUnless(0 == 1,
598#                        'imaginary file did not raise error!')
599
600                       
601  ###################### .XYA ##############################
602       
603    def test_export_xya_file(self):
604#        dict = {}
605        att_dict = {}
606        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
607        att_dict['elevation'] = array([10.0, 0.0, 10.4])
608        att_dict['brightness'] = array([10.0, 0.0, 10.4])
609#        dict['attributelist'] = att_dict
610        geo_reference=Geo_reference(56,1.9,1.9)
611       
612       
613        fileName = tempfile.mktemp(".xya")
614        G = Geospatial_data(pointlist, att_dict, geo_reference)
615        G.export_points_file(fileName, False)
616
617#        dict2 = import_points_file(fileName)
618        results = Geospatial_data(file_name = fileName)
619        #print "fileName",fileName
620        os.remove(fileName)
621       
622        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
623        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
624        answer = [10.0, 0.0, 10.4]
625        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
626        #print "dict2['geo_reference']",dict2['geo_reference']
627        self.failUnless(results.get_geo_reference() == geo_reference,
628                         'test_writepts failed. Test geo_reference')
629
630    def test_export_xya_file2(self):
631        """test absolute xya file
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       
638        fileName = tempfile.mktemp(".xya")
639        G = Geospatial_data(pointlist, att_dict)
640        G.export_points_file(fileName)
641        results = Geospatial_data(file_name = fileName)
642#        dict2 = import_points_file(fileName)
643        os.remove(fileName)
644       
645        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
646        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
647        answer = [10.0, 0.0, 10.4]
648        assert allclose(results.get_attributes('brightness'), answer)
649
650    def test_export_xya_file3(self):
651        """test absolute xya file with geo_ref
652        """
653        att_dict = {}
654        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
655        att_dict['elevation'] = array([10.0, 0.0, 10.4])
656        att_dict['brightness'] = array([10.0, 0.0, 10.4])
657        geo_reference=Geo_reference(56,1.9,1.9)
658       
659       
660        fileName = tempfile.mktemp(".xya")
661        G = Geospatial_data(pointlist, att_dict, geo_reference)
662       
663        G.export_points_file(fileName, absolute=True)
664       
665        results = Geospatial_data(file_name = fileName)
666        os.remove(fileName)
667
668        assert allclose(results.get_data_points(),
669                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
670        assert allclose(results.get_attributes(attribute_name='elevation'),
671                         [10.0, 0.0, 10.4])
672        answer = [10.0, 0.0, 10.4]
673        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
674        self.failUnless(results.get_geo_reference() == geo_reference,
675                         'test_writepts failed. Test geo_reference')                         
676                       
677                       
678                       
679    def test_new_export_pts_file(self):
680        att_dict = {}
681        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
682        att_dict['elevation'] = array([10.1, 0.0, 10.4])
683        att_dict['brightness'] = array([10.0, 1.0, 10.4])
684       
685        fileName = tempfile.mktemp(".pts")
686       
687        G = Geospatial_data(pointlist, att_dict)
688       
689        G.export_points_file(fileName)
690
691        results = Geospatial_data(file_name = fileName)
692
693        os.remove(fileName)
694       
695        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
696        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
697        answer = [10.0, 1.0, 10.4]
698        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
699
700    def test_new_export_absolute_pts_file(self):
701        att_dict = {}
702        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
703        att_dict['elevation'] = array([10.1, 0.0, 10.4])
704        att_dict['brightness'] = array([10.0, 1.0, 10.4])
705        geo_ref = Geo_reference(50, 25, 55)
706       
707        fileName = tempfile.mktemp(".pts")
708       
709        G = Geospatial_data(pointlist, att_dict, geo_ref)
710       
711        G.export_points_file(fileName, absolute=True)
712
713        results = Geospatial_data(file_name = fileName)
714
715        os.remove(fileName)
716       
717        assert allclose(results.get_data_points(), G.get_data_points(True))
718        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
719        answer = [10.0, 1.0, 10.4]
720        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
721
722    def test_loadpts(self):
723       
724        from Scientific.IO.NetCDF import NetCDFFile
725
726        fileName = tempfile.mktemp(".pts")
727        # NetCDF file definition
728        outfile = NetCDFFile(fileName, 'w')
729       
730        # dimension definitions
731        outfile.createDimension('number_of_points', 3)   
732        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
733   
734        # variable definitions
735        outfile.createVariable('points', Float, ('number_of_points',
736                                                 'number_of_dimensions'))
737        outfile.createVariable('elevation', Float, ('number_of_points',))
738   
739        # Get handles to the variables
740        points = outfile.variables['points']
741        elevation = outfile.variables['elevation']
742 
743        points[0, :] = [1.0,0.0]
744        elevation[0] = 10.0 
745        points[1, :] = [0.0,1.0]
746        elevation[1] = 0.0 
747        points[2, :] = [1.0,0.0]
748        elevation[2] = 10.4   
749
750        outfile.close()
751       
752        results = Geospatial_data(file_name = fileName)
753        os.remove(fileName)
754        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
755        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
756        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
757       
758    def test_writepts(self):
759        att_dict = {}
760        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
761        att_dict['elevation'] = array([10.0, 0.0, 10.4])
762        att_dict['brightness'] = array([10.0, 0.0, 10.4])
763        geo_reference=Geo_reference(56,1.9,1.9)
764       
765        fileName = tempfile.mktemp(".pts")
766       
767        G = Geospatial_data(pointlist, att_dict, geo_reference)
768       
769        G.export_points_file(fileName, False)
770       
771        results = Geospatial_data(file_name = fileName)
772
773        os.remove(fileName)
774
775        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
776        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
777        answer = [10.0, 0.0, 10.4]
778        assert allclose(results.get_attributes('brightness'), answer)
779
780       
781        self.failUnless(geo_reference == geo_reference,
782                         'test_writepts failed. Test geo_reference')
783       
784 ########################## BAD .PTS ##########################         
785
786    def test_load_bad_no_file_pts(self):
787        import os
788        import tempfile
789       
790        fileName = tempfile.mktemp(".pts")
791        #print fileName
792        try:
793            results = Geospatial_data(file_name = fileName)
794#            dict = import_points_file(fileName)
795        except IOError:
796            pass
797        else:
798            msg = 'imaginary file did not raise error!'
799            raise msg
800#            self.failUnless(0 == 1,
801#                        'imaginary file did not raise error!')
802
803
804    def test_create_from_pts_file(self):
805       
806        from Scientific.IO.NetCDF import NetCDFFile
807
808#        fileName = tempfile.mktemp(".pts")
809        FN = 'test_points.pts'
810        # NetCDF file definition
811        outfile = NetCDFFile(FN, 'w')
812       
813        # dimension definitions
814        outfile.createDimension('number_of_points', 3)   
815        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
816   
817        # variable definitions
818        outfile.createVariable('points', Float, ('number_of_points',
819                                                 'number_of_dimensions'))
820        outfile.createVariable('elevation', Float, ('number_of_points',))
821   
822        # Get handles to the variables
823        points = outfile.variables['points']
824        elevation = outfile.variables['elevation']
825 
826        points[0, :] = [1.0,0.0]
827        elevation[0] = 10.0 
828        points[1, :] = [0.0,1.0]
829        elevation[1] = 0.0 
830        points[2, :] = [1.0,0.0]
831        elevation[2] = 10.4   
832
833        outfile.close()
834
835        G = Geospatial_data(file_name = FN)
836
837        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
838        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
839
840        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
841        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
842        os.remove(FN)
843
844    def test_create_from_pts_file_with_geo(self):
845        """This test reveals if Geospatial data is correctly instantiated from a pts file.
846        """
847       
848        from Scientific.IO.NetCDF import NetCDFFile
849
850        FN = 'test_points.pts'
851        # NetCDF file definition
852        outfile = NetCDFFile(FN, 'w')
853
854        # Make up an arbitrary georef
855        xll = 0.1
856        yll = 20
857        geo_reference=Geo_reference(56, xll, yll)
858        geo_reference.write_NetCDF(outfile)
859
860        # dimension definitions
861        outfile.createDimension('number_of_points', 3)   
862        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
863   
864        # variable definitions
865        outfile.createVariable('points', Float, ('number_of_points',
866                                                 'number_of_dimensions'))
867        outfile.createVariable('elevation', Float, ('number_of_points',))
868   
869        # Get handles to the variables
870        points = outfile.variables['points']
871        elevation = outfile.variables['elevation']
872
873        points[0, :] = [1.0,0.0]
874        elevation[0] = 10.0 
875        points[1, :] = [0.0,1.0]
876        elevation[1] = 0.0 
877        points[2, :] = [1.0,0.0]
878        elevation[2] = 10.4   
879
880        outfile.close()
881
882        G = Geospatial_data(file_name = FN)
883
884        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
885        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
886
887        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
888                                              [0.0+xll, 1.0+yll],
889                                              [1.0+xll, 0.0+yll]])
890       
891        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
892        os.remove(FN)
893
894       
895    def test_add_(self):
896        '''adds an xya and pts files, reads the files and adds them
897           checking results are correct
898        '''
899
900        #FIXME (Ole): I changed the georeference arbitrarily, but the test still passes.
901       
902        # create files
903        att_dict1 = {}
904        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
905        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
906        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
907        geo_reference1 = Geo_reference(56, 2.0, 1.0)
908       
909        att_dict2 = {}
910        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
911        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
912        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
913        geo_reference2 = Geo_reference(56, 1.0, 2.0) #FIXME (Ole): I changed this, why does it still pass
914        #OK - it fails now, after the fix revealed by test_create_from_pts_file_with_geo')
915
916        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
917        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
918       
919        fileName1 = tempfile.mktemp(".xya")
920        fileName2 = tempfile.mktemp(".pts")
921
922        #makes files
923        G1.export_points_file(fileName1)
924        G2.export_points_file(fileName2)
925       
926        # add files
927       
928        G3 = Geospatial_data(file_name = fileName1)
929        G4 = Geospatial_data(file_name = fileName2)
930       
931        G = G3 + G4
932
933       
934        #read results
935#        print'res', G.get_data_points()
936#        print'res1', G.get_data_points(False)
937        assert allclose(G.get_data_points(),
938                        [[ 3.0, 1.0], [ 2.0, 2.0],
939                         [ 3.0, 1.0], [ 3.0, 3.0],
940                         [ 2.0, 4.0], [ 3.0, 3.0]])
941                         
942        assert allclose(G.get_attributes(attribute_name='elevation'),
943                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
944       
945        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
946        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
947       
948        self.failUnless(G.get_geo_reference() == geo_reference1,
949                         'test_writepts failed. Test geo_reference')
950                         
951        os.remove(fileName1)
952        os.remove(fileName2)
953       
954    def test_ensure_absolute(self):
955        points = [[2.0, 0.0],[1.0, 1.0],
956                         [2.0, 0.0],[2.0, 2.0],
957                         [1.0, 3.0],[2.0, 2.0]]
958        new_points = ensure_absolute(points)
959       
960        assert allclose(new_points, points)
961
962        points = array([[2.0, 0.0],[1.0, 1.0],
963                         [2.0, 0.0],[2.0, 2.0],
964                         [1.0, 3.0],[2.0, 2.0]])
965        new_points = ensure_absolute(points)
966       
967        assert allclose(new_points, points)
968       
969        ab_points = array([[2.0, 0.0],[1.0, 1.0],
970                         [2.0, 0.0],[2.0, 2.0],
971                         [1.0, 3.0],[2.0, 2.0]])
972       
973        mesh_origin = (56, 290000, 618000) #zone, easting, northing
974
975        data_points = zeros((ab_points.shape), Float)
976        #Shift datapoints according to new origins
977        for k in range(len(ab_points)):
978            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
979            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
980        #print "data_points",data_points     
981        new_points = ensure_absolute(data_points,
982                                             geo_reference=mesh_origin)
983        #print "new_points",new_points
984        #print "ab_points",ab_points
985           
986        assert allclose(new_points, ab_points)
987
988        geo = Geo_reference(56,67,-56)
989
990        data_points = geo.change_points_geo_ref(ab_points)   
991        new_points = ensure_absolute(data_points,
992                                             geo_reference=geo)
993        #print "new_points",new_points
994        #print "ab_points",ab_points
995           
996        assert allclose(new_points, ab_points)
997
998
999        geo_reference = Geo_reference(56, 100, 200)
1000        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1001        points = geo_reference.change_points_geo_ref(ab_points)
1002        attributes = [2, 4]
1003        #print "geo in points", points
1004        G = Geospatial_data(points, attributes,
1005                            geo_reference=geo_reference)
1006         
1007        new_points = ensure_absolute(G)
1008        #print "new_points",new_points
1009        #print "ab_points",ab_points
1010           
1011        assert allclose(new_points, ab_points)
1012
1013       
1014    def test_ensure_geospatial(self):
1015        points = [[2.0, 0.0],[1.0, 1.0],
1016                         [2.0, 0.0],[2.0, 2.0],
1017                         [1.0, 3.0],[2.0, 2.0]]
1018        new_points = ensure_absolute(points)
1019       
1020        assert allclose(new_points, points)
1021
1022        points = array([[2.0, 0.0],[1.0, 1.0],
1023                         [2.0, 0.0],[2.0, 2.0],
1024                         [1.0, 3.0],[2.0, 2.0]])
1025        new_points = ensure_absolute(points)
1026       
1027        assert allclose(new_points, points)
1028       
1029        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1030                         [2.0, 0.0],[2.0, 2.0],
1031                         [1.0, 3.0],[2.0, 2.0]])
1032       
1033        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1034
1035        data_points = zeros((ab_points.shape), Float)
1036        #Shift datapoints according to new origins
1037        for k in range(len(ab_points)):
1038            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1039            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1040        #print "data_points",data_points     
1041        new_geospatial = ensure_geospatial(data_points,
1042                                             geo_reference=mesh_origin)
1043        new_points = new_geospatial.get_data_points(absolute=True)
1044        #print "new_points",new_points
1045        #print "ab_points",ab_points
1046           
1047        assert allclose(new_points, ab_points)
1048
1049        geo = Geo_reference(56,67,-56)
1050
1051        data_points = geo.change_points_geo_ref(ab_points)   
1052        new_geospatial = ensure_geospatial(data_points,
1053                                             geo_reference=geo)
1054        new_points = new_geospatial.get_data_points(absolute=True)
1055        #print "new_points",new_points
1056        #print "ab_points",ab_points
1057           
1058        assert allclose(new_points, ab_points)
1059
1060
1061        geo_reference = Geo_reference(56, 100, 200)
1062        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1063        points = geo_reference.change_points_geo_ref(ab_points)
1064        attributes = [2, 4]
1065        #print "geo in points", points
1066        G = Geospatial_data(points, attributes,
1067                            geo_reference=geo_reference)
1068         
1069        new_geospatial  = ensure_geospatial(G)
1070        new_points = new_geospatial.get_data_points(absolute=True)
1071        #print "new_points",new_points
1072        #print "ab_points",ab_points
1073           
1074        assert allclose(new_points, ab_points)
1075       
1076    def test_isinstance(self):
1077
1078        import os
1079       
1080        fileName = tempfile.mktemp(".xya")
1081        file = open(fileName,"w")
1082        file.write("  elevation   speed \n\
10831.0 0.0 10.0 0.0\n\
10840.0 1.0 0.0 10.0\n\
10851.0 0.0 10.4 40.0\n\
1086#geocrap\n\
108756\n\
108856.6\n\
10893\n")
1090        file.close()
1091
1092        results = Geospatial_data(fileName)
1093       
1094        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1095        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
1096        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
1097
1098        os.remove(fileName)
1099       
1100    def test_delimiter(self):
1101       
1102        try:
1103            G = Geospatial_data(delimiter=',')
1104#            results = Geospatial_data(file_name = fileName)
1105#            dict = import_points_file(fileName)
1106        except ValueError:
1107            pass
1108        else:
1109            msg = 'Instance with No fileName but has a delimiter\
1110                  did not raise error!'
1111            raise msg
1112
1113    def test_no_constructors(self):
1114       
1115        try:
1116            G = Geospatial_data()
1117#            results = Geospatial_data(file_name = fileName)
1118#            dict = import_points_file(fileName)
1119        except ValueError:
1120            pass
1121        else:
1122            msg = 'Instance must have a filename or data points'
1123            raise msg       
1124
1125    def test_check_geo_reference(self):
1126        """
1127        checks geo reference details are OK. eg can be called '#geo reference'
1128        if not throws a clear error message
1129        """
1130        import os
1131        fileName = tempfile.mktemp(".xya")
1132        file = open(fileName,"w")
1133        file.write("  elevation  \n\
11341.0 0.0 10.0\n\
11350.0 1.0 0.0\n\
11361.0 0.0 10.4\n\
1137#ge oreference\n\
113856\n\
11391.1\n\
11401.0\n")
1141
1142        file.close()
1143        results = Geospatial_data(fileName)
1144        assert allclose(results.get_geo_reference().get_xllcorner(), 1.1)
1145        assert allclose(results.get_geo_reference().get_yllcorner(), 1.0)
1146
1147        os.remove(fileName)
1148       
1149        fileName = tempfile.mktemp(".xya")
1150        file = open(fileName,"w")
1151        file.write("  elevation  \n\
11521.0 0.0 10.0\n\
11530.0 1.0 0.0\n\
11541.0 0.0 10.4\n")
1155
1156        file.close()
1157        results = Geospatial_data(fileName)
1158       
1159        os.remove(fileName)
1160       
1161    def test_check_geo_reference1(self):
1162        """
1163        checks geo reference details are OK. eg can be called '#geo reference'
1164        if not throws a clear error message
1165        """
1166        import os
1167        fileName = tempfile.mktemp(".xya")
1168        file = open(fileName,"w")
1169        file.write("  elevation  \n\
11701.0 0.0 10.0\n\
11710.0 1.0 0.0\n\
11721.0 0.0 10.4\n\
1173#geo t t\n\
117456\n\
11751.1\n"
1176)
1177        file.close()
1178
1179        try:
1180            results = Geospatial_data(fileName, delimiter = " ")
1181        except IOError:
1182            pass
1183        else:
1184            msg = 'Geo reference data format is incorrect'
1185            raise msg       
1186
1187
1188        os.remove(fileName)
1189
1190
1191if __name__ == "__main__":
1192    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_set_geo_reference')
1193    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1194    runner = unittest.TextTestRunner()
1195    runner.run(suite)
1196
1197   
Note: See TracBrowser for help on using the repository browser.