source: inundation/geospatial_data/test_geospatial_data.py @ 3158

Last change on this file since 3158 was 3154, checked in by duncan, 19 years ago

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