source: anuga_core/source/anuga/geospatial_data/test_geospatial_data.py @ 3745

Last change on this file since 3745 was 3738, checked in by ole, 17 years ago

Added attributes into clip functions and tested
File size: 57.0 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 anuga.geospatial_data.geospatial_data import *
11from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
12from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
13
14
15class Test_Geospatial_data(unittest.TestCase):
16    def setUp(self):
17        pass
18
19    def tearDown(self):
20        pass
21
22
23    def test_0(self):
24        """test_0(self):
25       
26        Basic points
27        """
28        from anuga.coordinate_transforms.geo_reference import Geo_reference
29       
30        points = [[1.0, 2.1], [3.0, 5.3]]
31        G = Geospatial_data(points)
32
33        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
34
35        #Check __repr__
36
37        rep = `G`
38        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
39        assert rep == ref
40
41
42        #Check getter
43        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
44       
45        #Check defaults
46        assert G.attributes is None
47       
48        assert G.geo_reference.zone == Geo_reference().zone
49        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
50        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
51       
52
53    def test_1(self):
54        points = [[1.0, 2.1], [3.0, 5.3]]
55        attributes = [2, 4]
56        G = Geospatial_data(points, attributes)       
57
58        assert G.attributes.keys()[0] == 'attribute'
59        assert allclose(G.attributes.values()[0], [2, 4])
60       
61
62    def test_2(self):
63        from anuga.coordinate_transforms.geo_reference import Geo_reference
64        points = [[1.0, 2.1], [3.0, 5.3]]
65        attributes = [2, 4]
66        G = Geospatial_data(points, attributes,
67                            geo_reference=Geo_reference(56, 100, 200))
68
69        assert G.geo_reference.zone == 56
70        assert G.geo_reference.xllcorner == 100
71        assert G.geo_reference.yllcorner == 200
72
73
74    def test_get_attributes_1(self):
75        from anuga.coordinate_transforms.geo_reference import Geo_reference
76        points = [[1.0, 2.1], [3.0, 5.3]]
77        attributes = [2, 4]
78        G = Geospatial_data(points, attributes,
79                            geo_reference=Geo_reference(56, 100, 200))
80
81
82        P = G.get_data_points(absolute=False)
83        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
84
85        P = G.get_data_points(absolute=True)
86        assert allclose(P, [[101.0, 202.1], [103.0, 205.3]])       
87
88        V = G.get_attributes() #Simply get them
89        assert allclose(V, [2, 4])
90
91        V = G.get_attributes('attribute') #Get by name
92        assert allclose(V, [2, 4])
93
94    def test_get_attributes_2(self):
95        """Multiple attributes
96        """
97       
98        from anuga.coordinate_transforms.geo_reference import Geo_reference
99        points = [[1.0, 2.1], [3.0, 5.3]]
100        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
101        G = Geospatial_data(points, attributes,
102                            geo_reference=Geo_reference(56, 100, 200),
103                            default_attribute_name='a1')
104
105
106        P = G.get_data_points(absolute=False)
107        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
108       
109        V = G.get_attributes() #Get default attribute
110        assert allclose(V, [2, 4])
111
112        V = G.get_attributes('a0') #Get by name
113        assert allclose(V, [0, 0])
114
115        V = G.get_attributes('a1') #Get by name
116        assert allclose(V, [2, 4])
117
118        V = G.get_attributes('a2') #Get by name
119        assert allclose(V, [79.4, -7])
120
121        try:
122            V = G.get_attributes('hdnoatedu') #Invalid
123        except AssertionError:
124            pass
125        else:
126            raise 'Should have raised exception' 
127
128    def test_get_data_points(self):
129        points_ab = [[12.5,34.7],[-4.5,-60.0]]
130        x_p = -10
131        y_p = -40
132        geo_ref = Geo_reference(56, x_p, y_p)
133        points_rel = geo_ref.change_points_geo_ref(points_ab)
134       
135        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
136
137        results = spatial.get_data_points(absolute=False)
138       
139        assert allclose(results, points_rel)
140       
141        x_p = -1770
142        y_p = 4.01
143        geo_ref = Geo_reference(56, x_p, y_p)
144        points_rel = geo_ref.change_points_geo_ref(points_ab)
145        results = spatial.get_data_points \
146                  ( geo_reference=geo_ref)
147       
148        assert allclose(results, points_rel)
149
150       
151    def test_set_geo_reference(self):
152        points_ab = [[12.5,34.7],[-4.5,-60.0]]
153        x_p = -10
154        y_p = -40
155        geo_ref = Geo_reference(56, x_p, y_p)
156        points_rel = geo_ref.change_points_geo_ref(points_ab)
157        spatial = Geospatial_data(points_rel)
158        # since the geo_ref wasn't set
159        assert not allclose( points_ab, spatial.get_data_points(absolute=True))
160       
161        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
162        assert allclose( points_ab, spatial.get_data_points(absolute=True))
163       
164        x_p = 10
165        y_p = 400
166        new_geo_ref = Geo_reference(56, x_p, y_p)
167        spatial.set_geo_reference(new_geo_ref)
168        assert allclose( points_ab, spatial.get_data_points(absolute=True))
169       
170       
171       
172    def test_conversions_to_points_dict(self):
173        """test conversions to points_dict
174        """
175       
176        from anuga.coordinate_transforms.geo_reference import Geo_reference
177        points = [[1.0, 2.1], [3.0, 5.3]]
178        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
179        G = Geospatial_data(points, attributes,
180                            geo_reference=Geo_reference(56, 100, 200),
181                            default_attribute_name='a1')
182
183
184        points_dict = geospatial_data2points_dictionary(G)
185
186        assert points_dict.has_key('pointlist')
187        assert points_dict.has_key('attributelist')       
188        assert points_dict.has_key('geo_reference')
189
190        assert allclose( points_dict['pointlist'], points )
191
192        A = points_dict['attributelist']
193        assert A.has_key('a0')
194        assert A.has_key('a1')
195        assert A.has_key('a2')       
196
197        assert allclose( A['a0'], [0, 0] )
198        assert allclose( A['a1'], [2, 4] )       
199        assert allclose( A['a2'], [79.4, -7] )
200
201
202        geo = points_dict['geo_reference']
203        assert geo is G.geo_reference
204
205
206    def test_conversions_from_points_dict(self):
207        """test conversions from points_dict
208        """
209
210        from anuga.coordinate_transforms.geo_reference import Geo_reference
211       
212        points = [[1.0, 2.1], [3.0, 5.3]]
213        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
214
215        points_dict = {}
216        points_dict['pointlist'] = points
217        points_dict['attributelist'] = attributes
218        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
219       
220
221        G = points_dictionary2geospatial_data(points_dict)
222
223        P = G.get_data_points(absolute=False)
224        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
225       
226        #V = G.get_attribute_values() #Get default attribute
227        #assert allclose(V, [2, 4])
228
229        V = G.get_attributes('a0') #Get by name
230        assert allclose(V, [0, 0])
231
232        V = G.get_attributes('a1') #Get by name
233        assert allclose(V, [2, 4])
234
235        V = G.get_attributes('a2') #Get by name
236        assert allclose(V, [79.4, -7])
237
238    def test_add(self):
239        """ test the addition of two geospatical objects
240            no geo_reference see next test
241        """
242        points = [[1.0, 2.1], [3.0, 5.3]]
243        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
244        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
245        G1 = Geospatial_data(points, attributes)       
246        G2 = Geospatial_data(points, attributes1) 
247       
248#        g3 = geospatial_data2points_dictionary(G1)
249#        print 'g3=', g3
250       
251        G = G1 + G2
252
253        assert G.attributes.has_key('depth')
254        assert G.attributes.has_key('elevation')
255        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
256        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
257        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
258                                              [1.0, 2.1], [3.0, 5.3]])
259       
260    def test_addII(self):
261        """ test the addition of two geospatical objects
262            no geo_reference see next test
263        """
264        points = [[1.0, 2.1], [3.0, 5.3]]
265        attributes = {'depth':[2, 4]}
266        G1 = Geospatial_data(points, attributes) 
267       
268        points = [[5.0, 2.1], [3.0, 50.3]]
269        attributes = {'depth':[200, 400]}
270        G2 = Geospatial_data(points, attributes)
271       
272#        g3 = geospatial_data2points_dictionary(G1)
273#        print 'g3=', g3
274       
275        G = G1 + G2
276
277        assert G.attributes.has_key('depth') 
278        assert G.attributes.keys(), ['depth']
279        assert allclose(G.attributes['depth'], [2, 4, 200, 400])
280        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
281                                              [5.0, 2.1], [3.0, 50.3]])
282    def test_add_with_geo (self):
283        """
284        Difference in Geo_reference resolved
285        """
286        points1 = [[1.0, 2.1], [3.0, 5.3]]
287        points2 = [[5.0, 6.1], [6.0, 3.3]]
288        attributes1 = [2, 4]
289        attributes2 = [5, 76]
290        geo_ref1= Geo_reference(55, 1.0, 2.0)
291        geo_ref2 = Geo_reference(zone=55,
292                                 xllcorner=0.1,
293                                 yllcorner=3.0,
294                                 datum='wgs84',
295                                 projection='UTM',
296                                 units='m')
297                               
298        G1 = Geospatial_data(points1, attributes1, geo_ref1)
299        G2 = Geospatial_data(points2, attributes2, geo_ref2)
300
301        #Check that absolute values are as expected
302        P1 = G1.get_data_points(absolute=True)
303        assert allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
304
305        P2 = G2.get_data_points(absolute=True)
306        assert allclose(P2, [[5.1, 9.1], [6.1, 6.3]])       
307       
308        G = G1 + G2
309       
310        assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
311        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
312
313        P = G.get_data_points(absolute=True)
314
315        P_relative = G.get_data_points(absolute=False)
316       
317        assert allclose(P_relative, P - [0.1, 2.0])
318
319        assert allclose(P, concatenate( (P1,P2) ))
320        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
321                            [5.1, 9.1], [6.1, 6.3]])
322       
323
324
325       
326
327    def test_add_with_geo_absolute (self):
328        """
329        Difference in Geo_reference resolved
330        """
331        points1 = array([[2.0, 4.1], [4.0, 7.3]])
332        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
333        attributes1 = [2, 4]
334        attributes2 = [5, 76]
335        geo_ref1= Geo_reference(55, 1.0, 2.0)
336        geo_ref2 = Geo_reference(55, 2.0, 3.0)
337
338       
339                               
340        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
341                             attributes1, geo_ref1)
342       
343        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
344                             attributes2, geo_ref2)
345
346        #Check that absolute values are as expected
347        P1 = G1.get_data_points(absolute=True)
348        assert allclose(P1, points1)
349
350        P1 = G1.get_data_points(absolute=False)
351        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
352
353        P2 = G2.get_data_points(absolute=True)
354        assert allclose(P2, points2)
355
356        P2 = G2.get_data_points(absolute=False)
357        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
358       
359        G = G1 + G2
360       
361        assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
362        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
363
364        P = G.get_data_points(absolute=True)
365
366        P_relative = G.get_data_points(absolute=False)
367       
368        assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
369
370        assert allclose(P, concatenate( (points1,points2) ))
371                           
372       
373    def test_clip0(self):
374        """test_clip0(self):
375       
376        Test that point sets can be clipped by a polygon
377        """
378       
379        from anuga.coordinate_transforms.geo_reference import Geo_reference
380       
381        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
382                  [0, 0], [2.4, 3.3]]
383        G = Geospatial_data(points)
384
385        # First try the unit square   
386        U = [[0,0], [1,0], [1,1], [0,1]] 
387        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
388
389        # Then a more complex polygon
390        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
391        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
392        G = Geospatial_data(points)
393
394        assert allclose(G.clip(polygon).get_data_points(),
395                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
396
397    def test_clip0_with_attributes(self):
398        """test_clip0_with_attributes(self):
399       
400        Test that point sets with attributes can be clipped by a polygon
401        """
402       
403        from anuga.coordinate_transforms.geo_reference import Geo_reference
404       
405        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
406                  [0, 0], [2.4, 3.3]]
407
408        attributes = [2, -4, 5, 76, -2, 0.1, 3]
409        att_dict = {'att1': attributes,
410                    'att2': array(attributes)+1}
411       
412        G = Geospatial_data(points, att_dict)
413
414        # First try the unit square   
415        U = [[0,0], [1,0], [1,1], [0,1]] 
416        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
417        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
418        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])               
419
420        # Then a more complex polygon
421        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
422        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
423
424        # This time just one attribute
425        attributes = [2, -4, 5, 76, -2, 0.1]
426        G = Geospatial_data(points, attributes)
427
428        assert allclose(G.clip(polygon).get_data_points(),
429                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
430        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
431       
432
433    def test_clip1(self):
434        """test_clip1(self):
435       
436        Test that point sets can be clipped by a polygon given as
437        another Geospatial dataset
438        """
439       
440        from anuga.coordinate_transforms.geo_reference import Geo_reference
441       
442        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
443                  [0, 0], [2.4, 3.3]]
444        attributes = [2, -4, 5, 76, -2, 0.1, 3]
445        att_dict = {'att1': attributes,
446                    'att2': array(attributes)+1}
447        G = Geospatial_data(points, att_dict)
448       
449        # First try the unit square   
450        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
451        assert allclose(G.clip(U).get_data_points(),
452                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
453
454        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
455        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                       
456       
457        # Then a more complex polygon
458        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
459        attributes = [2, -4, 5, 76, -2, 0.1]       
460        G = Geospatial_data(points, attributes)
461        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
462       
463
464        assert allclose(G.clip(polygon).get_data_points(),
465                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
466        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
467       
468
469    def test_clip0_outside(self):
470        """test_clip0_outside(self):
471       
472        Test that point sets can be clipped outside of a polygon
473        """
474       
475        from anuga.coordinate_transforms.geo_reference import Geo_reference
476       
477        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
478                  [0, 0], [2.4, 3.3]]
479        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
480        G = Geospatial_data(points, attributes)
481
482        # First try the unit square   
483        U = [[0,0], [1,0], [1,1], [0,1]]
484        assert allclose(G.clip_outside(U).get_data_points(),
485                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
486        #print G.clip_outside(U).get_attributes()
487        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])       
488       
489
490        # Then a more complex polygon
491        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
492        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
493        attributes = [2, -4, 5, 76, -2, 0.1]       
494        G = Geospatial_data(points, attributes)
495
496        assert allclose(G.clip_outside(polygon).get_data_points(),
497                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
498        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])               
499
500
501    def test_clip1_outside(self):
502        """test_clip1_outside(self):
503       
504        Test that point sets can be clipped outside of a polygon given as
505        another Geospatial dataset
506        """
507       
508        from anuga.coordinate_transforms.geo_reference import Geo_reference
509       
510        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
511                  [0, 0], [2.4, 3.3]]
512        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
513        G = Geospatial_data(points, attributes)       
514
515        # First try the unit square   
516        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
517        assert allclose(G.clip_outside(U).get_data_points(),
518                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
519        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])       
520
521        # Then a more complex polygon
522        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
523        attributes = [2, -4, 5, 76, -2, 0.1]       
524        G = Geospatial_data(points, attributes)
525
526        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
527       
528
529        assert allclose(G.clip_outside(polygon).get_data_points(),
530                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
531        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
532       
533
534
535    def test_clip1_inside_outside(self):
536        """test_clip1_inside_outside(self):
537       
538        Test that point sets can be clipped outside of a polygon given as
539        another Geospatial dataset
540        """
541       
542        from anuga.coordinate_transforms.geo_reference import Geo_reference
543       
544        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
545                  [0, 0], [2.4, 3.3]]
546        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
547        G = Geospatial_data(points, attributes)
548
549        # First try the unit square   
550        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
551        G1 = G.clip(U)
552        assert allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
553        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
554       
555        G2 = G.clip_outside(U)
556        assert allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
557                                              [3.0, 5.3], [2.4, 3.3]])
558        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])               
559
560       
561        # New ordering
562        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\
563                     [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
564
565        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
566       
567        assert allclose((G1+G2).get_data_points(), new_points)
568        assert allclose((G1+G2).get_attributes(), new_attributes)
569
570        G = G1+G2
571        FN = 'test_combine.pts'
572        G.export_points_file(FN)
573
574
575        # Read it back in
576        G3 = Geospatial_data(FN)
577
578
579        # Check result
580        assert allclose(G3.get_data_points(), new_points)       
581        assert allclose(G3.get_attributes(), new_attributes)       
582       
583        os.remove(FN)
584
585       
586    def test_create_from_xya_file(self):
587        """Check that object can be created from a points file (.pts and .xya)
588        """
589
590        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
591        attributes = [2, 4, 5, 76]
592        '''
593        # Use old pointsdict format
594        pointsdict = {'pointlist': points,
595                      'attributelist': {'att1': attributes,
596                                        'att2': array(attributes) + 1}}
597        '''
598        att_dict = {'att1': attributes,
599                    'att2': array(attributes) +1}
600                   
601        # Create points as an xya file
602        FN = 'test_points.xya'
603        G1 = Geospatial_data(points, att_dict)
604        G1.export_points_file(FN)
605#        G1.export_points_file(ofile)
606
607        #Create object from file
608        G = Geospatial_data(file_name = FN)
609       
610        assert allclose(G.get_data_points(), points)
611        assert allclose(G.get_attributes('att1'), attributes)
612        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
613       
614        os.remove(FN)
615
616    def test_create_from_xya_file1(self):
617        """
618        Check that object can be created from an Absolute xya file
619        """
620
621        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
622        attributes = [2, 4, 5, 76]
623
624        att_dict = {'att1': attributes,
625                    'att2': array(attributes) +1}
626
627        geo_ref = Geo_reference(56, 10, 5)
628                   
629        # Create points as an xya file
630        FN = 'test_points.xya'
631        G1 = Geospatial_data(points, att_dict, geo_ref)
632
633        G1.export_points_file(FN, absolute=True)
634
635        #Create object from file
636        G = Geospatial_data(file_name = FN)
637       
638        assert allclose(G.get_data_points(absolute=True), 
639                        G1.get_data_points(absolute=True))
640        assert allclose(G.get_attributes('att1'), attributes)
641        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
642       
643        os.remove(FN)
644       
645    def test_loadxya(self):
646        """
647        comma delimited
648        """
649        fileName = tempfile.mktemp(".xya")
650        file = open(fileName,"w")
651        file.write("elevation  , speed \n\
6521.0, 0.0, 10.0, 0.0\n\
6530.0, 1.0, 0.0, 10.0\n\
6541.0, 0.0, 10.4, 40.0\n")
655        file.close()
656        results = Geospatial_data(fileName, delimiter=',')
657        os.remove(fileName)
658#        print 'data', results.get_data_points()
659        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
660        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
661        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
662
663    def test_loadxya2(self):
664        """
665        space delimited
666        """
667        import os
668       
669        fileName = tempfile.mktemp(".xya")
670        file = open(fileName,"w")
671        file.write("  elevation   speed \n\
6721.0 0.0 10.0 0.0\n\
6730.0 1.0 0.0 10.0\n\
6741.0 0.0 10.4 40.0\n")
675        file.close()
676
677        results = Geospatial_data(fileName, delimiter=' ')
678
679        os.remove(fileName)
680
681        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
682        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
683        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
684     
685    def test_loadxya3(self):
686        """
687        space delimited
688        """
689        import os
690       
691        fileName = tempfile.mktemp(".xya")
692        file = open(fileName,"w")
693        file.write("  elevation   speed \n\
6941.0 0.0 10.0 0.0\n\
6950.0 1.0 0.0 10.0\n\
6961.0 0.0 10.4 40.0\n\
697#geocrap\n\
69856\n\
69956.6\n\
7003\n")
701        file.close()
702
703        results = Geospatial_data(fileName, delimiter=' ')
704
705        os.remove(fileName)
706        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
707        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
708        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
709
710    def BADtest_loadxya4(self):
711        """
712        comma delimited
713        """
714        fileName = tempfile.mktemp(".xya")
715        file = open(fileName,"w")
716        file.write("elevation  , speed \n\
7171.0, 0.0, splat, 0.0\n\
7180.0, 1.0, 0.0, 10.0\n\
7191.0, 0.0, 10.4, 40.0\n")
720        file.close()
721        results = Geospatial_data(fileName, delimiter=',')
722        os.remove(fileName)
723#        print 'data', results.get_data_points()
724        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
725        assert allclose(results.get_attributes(attribute_name='elevation'), ["splat", 0.0, 10.4])
726        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
727       
728    def test_read_write_points_file_bad2(self):
729        att_dict = {}
730        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
731        att_dict['elevation'] = array([10.0, 0.0, 10.4])
732        att_dict['brightness'] = array([10.0, 0.0, 10.4])
733        geo_reference=Geo_reference(56,1.9,1.9)
734       
735        G = Geospatial_data(pointlist, att_dict, geo_reference)
736       
737        try:
738            G.export_points_file("_???/yeah.xya")
739           
740        except IOError:
741            pass
742        else:
743            msg = 'bad points file extension did not raise error!'
744            raise msg
745#            self.failUnless(0 == 1,
746#                        'bad points file extension did not raise error!')
747                   
748    def test_loadxy_bad(self):
749        import os
750       
751        fileName = tempfile.mktemp(".xya")
752        file = open(fileName,"w")
753        file.write("  elevation   \n\
7541.0 0.0 10.0 0.0\n\
7550.0 1.0 0.0 10.0\n\
7561.0 0.0 10.4 40.0\n")
757        file.close()
758        #print fileName
759        try:
760            results = Geospatial_data(fileName, delimiter=' ')
761        except IOError:
762            pass
763        else:
764            msg = 'bad xya file did not raise error!'
765            raise msg
766#            self.failUnless(0 == 1,
767#                        'bad xya file did not raise error!')
768        os.remove(fileName)
769       
770    def test_loadxy_bad2(self):
771        import os
772       
773        fileName = tempfile.mktemp(".xya")
774        file = open(fileName,"w")
775        file.write("elevation\n\
7761.0 0.0 10.0 \n\
7770.0 1.0\n\
7781.0 \n")
779        file.close()
780        #print fileName
781        try:
782            results = Geospatial_data(fileName, delimiter=' ')
783        except IOError:
784            pass
785        else:
786            msg = 'bad xya file did not raise error!'
787            raise msg
788        os.remove(fileName)
789   
790    def test_loadxy_bad3(self):
791        """
792        specifying wrong delimiter
793        """
794        import os
795       
796        fileName = tempfile.mktemp(".xya")
797        file = open(fileName,"w")
798        file.write("  elevation  , speed \n\
7991.0, 0.0, 10.0, 0.0\n\
8000.0, 1.0, 0.0, 10.0\n\
8011.0, 0.0, 10.4, 40.0\n")
802        file.close()
803        try:
804            results = Geospatial_data(fileName, delimiter=' ')
805        except IOError:
806            pass
807        else:
808            msg = 'bad xya file did not raise error!'
809            raise msg
810        os.remove(fileName)
811     
812    def test_loadxy_bad4(self):
813        """
814         specifying wrong delimiter
815        """
816        import os
817        fileName = tempfile.mktemp(".xya")
818        file = open(fileName,"w")
819        file.write("  elevation   speed \n\
8201.0 0.0 10.0 0.0\n\
8210.0 1.0 0.0 10.0\n\
8221.0 0.0 10.4 40.0\n\
823#geocrap\n\
82456\n\
82556.6\n\
8263\n"
827)
828        file.close()
829        try:
830            results = Geospatial_data(fileName, delimiter=',')
831        except IOError:
832            pass
833        else:
834            msg = 'bad xya file did not raise error!'
835            raise msg
836
837        os.remove(fileName)
838
839    def test_loadxy_bad5(self):
840        """
841        specifying wrong delimiter
842        """
843        import os
844       
845        fileName = tempfile.mktemp(".xya")
846        file = open(fileName,"w")
847        file.write("  elevation   speed \n\
8481.0 0.0 10.0 0.0\n\
8490.0 1.0 0.0 10.0\n\
8501.0 0.0 10.4 40.0\n\
851#geocrap\n\
852crap")
853        file.close()
854        try:
855#            dict = import_points_file(fileName,delimiter=' ')
856#            results = Geospatial_data()
857            results = Geospatial_data(fileName, delimiter=' ', verbose=True)
858#            results.import_points_file(fileName, delimiter=' ')
859        except IOError:
860            pass
861        else:
862            msg = 'bad xya file did not raise error!'
863            raise msg
864
865#            self.failUnless(0 ==1,
866#                        'bad xya file did not raise error!')   
867        os.remove(fileName)             
868
869    def test_loadxy_bad_no_file_xya(self):
870        import os
871       
872        fileName = tempfile.mktemp(".xya")
873        try:
874            results = Geospatial_data(fileName, delimiter=' ')
875        except IOError:
876            pass
877        else:
878            msg = 'imaginary file did not raise error!'
879            raise msg
880
881#        except IOError:
882#            pass
883#        else:
884#            self.failUnless(0 == 1,
885#                        'imaginary file did not raise error!')
886
887                       
888  ###################### .XYA ##############################
889       
890    def test_export_xya_file(self):
891#        dict = {}
892        att_dict = {}
893        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
894        att_dict['elevation'] = array([10.0, 0.0, 10.4])
895        att_dict['brightness'] = array([10.0, 0.0, 10.4])
896#        dict['attributelist'] = att_dict
897        geo_reference=Geo_reference(56,1.9,1.9)
898       
899       
900        fileName = tempfile.mktemp(".xya")
901        G = Geospatial_data(pointlist, att_dict, geo_reference)
902        G.export_points_file(fileName, False)
903
904#        dict2 = import_points_file(fileName)
905        results = Geospatial_data(file_name = fileName)
906        #print "fileName",fileName
907        os.remove(fileName)
908       
909        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
910        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
911        answer = [10.0, 0.0, 10.4]
912        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
913        #print "dict2['geo_reference']",dict2['geo_reference']
914        self.failUnless(results.get_geo_reference() == geo_reference,
915                         'test_writepts failed. Test geo_reference')
916
917    def test_export_xya_file2(self):
918        """test absolute xya file
919        """
920        att_dict = {}
921        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
922        att_dict['elevation'] = array([10.0, 0.0, 10.4])
923        att_dict['brightness'] = array([10.0, 0.0, 10.4])
924       
925        fileName = tempfile.mktemp(".xya")
926        G = Geospatial_data(pointlist, att_dict)
927        G.export_points_file(fileName)
928        results = Geospatial_data(file_name = fileName)
929#        dict2 = import_points_file(fileName)
930        os.remove(fileName)
931       
932        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
933        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
934        answer = [10.0, 0.0, 10.4]
935        assert allclose(results.get_attributes('brightness'), answer)
936
937    def test_export_xya_file3(self):
938        """test absolute xya file with geo_ref
939        """
940        att_dict = {}
941        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
942        att_dict['elevation'] = array([10.0, 0.0, 10.4])
943        att_dict['brightness'] = array([10.0, 0.0, 10.4])
944        geo_reference=Geo_reference(56,1.9,1.9)
945       
946       
947        fileName = tempfile.mktemp(".xya")
948        G = Geospatial_data(pointlist, att_dict, geo_reference)
949       
950        G.export_points_file(fileName, absolute=True)
951       
952        results = Geospatial_data(file_name = fileName)
953        os.remove(fileName)
954
955        assert allclose(results.get_data_points(),
956                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
957        assert allclose(results.get_attributes(attribute_name='elevation'),
958                         [10.0, 0.0, 10.4])
959        answer = [10.0, 0.0, 10.4]
960        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
961        self.failUnless(results.get_geo_reference() == geo_reference,
962                         'test_writepts failed. Test geo_reference')                         
963                       
964                       
965                       
966    def test_new_export_pts_file(self):
967        att_dict = {}
968        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
969        att_dict['elevation'] = array([10.1, 0.0, 10.4])
970        att_dict['brightness'] = array([10.0, 1.0, 10.4])
971       
972        fileName = tempfile.mktemp(".pts")
973       
974        G = Geospatial_data(pointlist, att_dict)
975       
976        G.export_points_file(fileName)
977
978        results = Geospatial_data(file_name = fileName)
979
980        os.remove(fileName)
981       
982        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
983        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
984        answer = [10.0, 1.0, 10.4]
985        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
986
987    def test_new_export_absolute_pts_file(self):
988        att_dict = {}
989        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
990        att_dict['elevation'] = array([10.1, 0.0, 10.4])
991        att_dict['brightness'] = array([10.0, 1.0, 10.4])
992        geo_ref = Geo_reference(50, 25, 55)
993       
994        fileName = tempfile.mktemp(".pts")
995       
996        G = Geospatial_data(pointlist, att_dict, geo_ref)
997       
998        G.export_points_file(fileName, absolute=True)
999
1000        results = Geospatial_data(file_name = fileName)
1001
1002        os.remove(fileName)
1003       
1004        assert allclose(results.get_data_points(), G.get_data_points(True))
1005        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
1006        answer = [10.0, 1.0, 10.4]
1007        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
1008
1009    def test_loadpts(self):
1010       
1011        from Scientific.IO.NetCDF import NetCDFFile
1012
1013        fileName = tempfile.mktemp(".pts")
1014        # NetCDF file definition
1015        outfile = NetCDFFile(fileName, 'w')
1016       
1017        # dimension definitions
1018        outfile.createDimension('number_of_points', 3)   
1019        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1020   
1021        # variable definitions
1022        outfile.createVariable('points', Float, ('number_of_points',
1023                                                 'number_of_dimensions'))
1024        outfile.createVariable('elevation', Float, ('number_of_points',))
1025   
1026        # Get handles to the variables
1027        points = outfile.variables['points']
1028        elevation = outfile.variables['elevation']
1029 
1030        points[0, :] = [1.0,0.0]
1031        elevation[0] = 10.0 
1032        points[1, :] = [0.0,1.0]
1033        elevation[1] = 0.0 
1034        points[2, :] = [1.0,0.0]
1035        elevation[2] = 10.4   
1036
1037        outfile.close()
1038       
1039        results = Geospatial_data(file_name = fileName)
1040        os.remove(fileName)
1041        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
1042        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1043        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
1044       
1045    def test_writepts(self):
1046        """test_writepts: Test that storage of x,y,attributes works
1047        """
1048        att_dict = {}
1049        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1050        att_dict['elevation'] = array([10.0, 0.0, 10.4])
1051        att_dict['brightness'] = array([10.0, 0.0, 10.4])
1052        geo_reference=Geo_reference(56,1.9,1.9)
1053
1054        # Test pts format
1055        fileName = tempfile.mktemp(".pts")
1056        G = Geospatial_data(pointlist, att_dict, geo_reference)
1057        G.export_points_file(fileName, False)
1058        results = Geospatial_data(file_name=fileName)
1059        os.remove(fileName)
1060
1061        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1062        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1063        answer = [10.0, 0.0, 10.4]
1064        assert allclose(results.get_attributes('brightness'), answer)
1065        self.failUnless(geo_reference == geo_reference,
1066                         'test_writepts failed. Test geo_reference')
1067
1068        # Test xya format
1069        fileName = tempfile.mktemp(".xya")
1070        G = Geospatial_data(pointlist, att_dict, geo_reference)
1071        G.export_points_file(fileName, False)
1072        results = Geospatial_data(file_name=fileName)
1073        os.remove(fileName)
1074        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1075        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1076        answer = [10.0, 0.0, 10.4]
1077        assert allclose(results.get_attributes('brightness'), answer)
1078        self.failUnless(geo_reference == geo_reference,
1079                         'test_writepts failed. Test geo_reference')
1080
1081    def test_writepts_no_attributes(self):
1082        """test_writepts_no_attributes: Test that storage of x,y alone works
1083        """
1084        att_dict = {}
1085        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1086        geo_reference=Geo_reference(56,1.9,1.9)
1087
1088        # Test pts format
1089        fileName = tempfile.mktemp(".pts")
1090        G = Geospatial_data(pointlist, None, geo_reference)
1091        G.export_points_file(fileName, False)
1092        results = Geospatial_data(file_name=fileName)
1093        os.remove(fileName)
1094
1095        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1096        self.failUnless(geo_reference == geo_reference,
1097                         'test_writepts failed. Test geo_reference')
1098
1099        # Test xya format
1100        fileName = tempfile.mktemp(".xya")
1101        G = Geospatial_data(pointlist, None, geo_reference)
1102        G.export_points_file(fileName, False)
1103        results = Geospatial_data(file_name=fileName)
1104        os.remove(fileName)
1105        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1106        self.failUnless(geo_reference == geo_reference,
1107                         'test_writepts failed. Test geo_reference')
1108
1109       
1110       
1111 ########################## BAD .PTS ##########################         
1112
1113    def test_load_bad_no_file_pts(self):
1114        import os
1115        import tempfile
1116       
1117        fileName = tempfile.mktemp(".pts")
1118        #print fileName
1119        try:
1120            results = Geospatial_data(file_name = fileName)
1121#            dict = import_points_file(fileName)
1122        except IOError:
1123            pass
1124        else:
1125            msg = 'imaginary file did not raise error!'
1126            raise msg
1127#            self.failUnless(0 == 1,
1128#                        'imaginary file did not raise error!')
1129
1130
1131    def test_create_from_pts_file(self):
1132       
1133        from Scientific.IO.NetCDF import NetCDFFile
1134
1135#        fileName = tempfile.mktemp(".pts")
1136        FN = 'test_points.pts'
1137        # NetCDF file definition
1138        outfile = NetCDFFile(FN, 'w')
1139       
1140        # dimension definitions
1141        outfile.createDimension('number_of_points', 3)   
1142        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1143   
1144        # variable definitions
1145        outfile.createVariable('points', Float, ('number_of_points',
1146                                                 'number_of_dimensions'))
1147        outfile.createVariable('elevation', Float, ('number_of_points',))
1148   
1149        # Get handles to the variables
1150        points = outfile.variables['points']
1151        elevation = outfile.variables['elevation']
1152 
1153        points[0, :] = [1.0,0.0]
1154        elevation[0] = 10.0 
1155        points[1, :] = [0.0,1.0]
1156        elevation[1] = 0.0 
1157        points[2, :] = [1.0,0.0]
1158        elevation[2] = 10.4   
1159
1160        outfile.close()
1161
1162        G = Geospatial_data(file_name = FN)
1163
1164        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
1165        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
1166
1167        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1168        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1169        os.remove(FN)
1170
1171    def test_create_from_pts_file_with_geo(self):
1172        """This test reveals if Geospatial data is correctly instantiated from a pts file.
1173        """
1174       
1175        from Scientific.IO.NetCDF import NetCDFFile
1176
1177        FN = 'test_points.pts'
1178        # NetCDF file definition
1179        outfile = NetCDFFile(FN, 'w')
1180
1181        # Make up an arbitrary georef
1182        xll = 0.1
1183        yll = 20
1184        geo_reference=Geo_reference(56, xll, yll)
1185        geo_reference.write_NetCDF(outfile)
1186
1187        # dimension definitions
1188        outfile.createDimension('number_of_points', 3)   
1189        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1190   
1191        # variable definitions
1192        outfile.createVariable('points', Float, ('number_of_points',
1193                                                 'number_of_dimensions'))
1194        outfile.createVariable('elevation', Float, ('number_of_points',))
1195   
1196        # Get handles to the variables
1197        points = outfile.variables['points']
1198        elevation = outfile.variables['elevation']
1199
1200        points[0, :] = [1.0,0.0]
1201        elevation[0] = 10.0 
1202        points[1, :] = [0.0,1.0]
1203        elevation[1] = 0.0 
1204        points[2, :] = [1.0,0.0]
1205        elevation[2] = 10.4   
1206
1207        outfile.close()
1208
1209        G = Geospatial_data(file_name = FN)
1210
1211        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
1212        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
1213
1214        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
1215                                              [0.0+xll, 1.0+yll],
1216                                              [1.0+xll, 0.0+yll]])
1217       
1218        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1219        os.remove(FN)
1220
1221       
1222    def test_add_(self):
1223        '''adds an xya and pts files, reads the files and adds them
1224           checking results are correct
1225        '''
1226
1227        # create files
1228        att_dict1 = {}
1229        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1230        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
1231        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
1232        geo_reference1 = Geo_reference(56, 2.0, 1.0)
1233       
1234        att_dict2 = {}
1235        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
1236        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
1237        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
1238        geo_reference2 = Geo_reference(56, 1.0, 2.0) 
1239
1240        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1241        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1242       
1243        fileName1 = tempfile.mktemp(".xya")
1244        fileName2 = tempfile.mktemp(".pts")
1245
1246        #makes files
1247        G1.export_points_file(fileName1)
1248        G2.export_points_file(fileName2)
1249       
1250        # add files
1251       
1252        G3 = Geospatial_data(file_name = fileName1)
1253        G4 = Geospatial_data(file_name = fileName2)
1254       
1255        G = G3 + G4
1256
1257       
1258        #read results
1259#        print'res', G.get_data_points()
1260#        print'res1', G.get_data_points(False)
1261        assert allclose(G.get_data_points(),
1262                        [[ 3.0, 1.0], [ 2.0, 2.0],
1263                         [ 3.0, 1.0], [ 3.0, 3.0],
1264                         [ 2.0, 4.0], [ 3.0, 3.0]])
1265                         
1266        assert allclose(G.get_attributes(attribute_name='elevation'),
1267                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
1268       
1269        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
1270        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
1271       
1272        self.failUnless(G.get_geo_reference() == geo_reference1,
1273                         'test_writepts failed. Test geo_reference')
1274                         
1275        os.remove(fileName1)
1276        os.remove(fileName2)
1277       
1278    def test_ensure_absolute(self):
1279        points = [[2.0, 0.0],[1.0, 1.0],
1280                         [2.0, 0.0],[2.0, 2.0],
1281                         [1.0, 3.0],[2.0, 2.0]]
1282        new_points = ensure_absolute(points)
1283       
1284        assert allclose(new_points, points)
1285
1286        points = array([[2.0, 0.0],[1.0, 1.0],
1287                         [2.0, 0.0],[2.0, 2.0],
1288                         [1.0, 3.0],[2.0, 2.0]])
1289        new_points = ensure_absolute(points)
1290       
1291        assert allclose(new_points, points)
1292       
1293        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1294                         [2.0, 0.0],[2.0, 2.0],
1295                         [1.0, 3.0],[2.0, 2.0]])
1296       
1297        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1298
1299        data_points = zeros((ab_points.shape), Float)
1300        #Shift datapoints according to new origins
1301        for k in range(len(ab_points)):
1302            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1303            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1304        #print "data_points",data_points     
1305        new_points = ensure_absolute(data_points,
1306                                             geo_reference=mesh_origin)
1307        #print "new_points",new_points
1308        #print "ab_points",ab_points
1309           
1310        assert allclose(new_points, ab_points)
1311
1312        geo = Geo_reference(56,67,-56)
1313
1314        data_points = geo.change_points_geo_ref(ab_points)   
1315        new_points = ensure_absolute(data_points,
1316                                             geo_reference=geo)
1317        #print "new_points",new_points
1318        #print "ab_points",ab_points
1319           
1320        assert allclose(new_points, ab_points)
1321
1322
1323        geo_reference = Geo_reference(56, 100, 200)
1324        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1325        points = geo_reference.change_points_geo_ref(ab_points)
1326        attributes = [2, 4]
1327        #print "geo in points", points
1328        G = Geospatial_data(points, attributes,
1329                            geo_reference=geo_reference)
1330         
1331        new_points = ensure_absolute(G)
1332        #print "new_points",new_points
1333        #print "ab_points",ab_points
1334           
1335        assert allclose(new_points, ab_points)
1336
1337       
1338        fileName = tempfile.mktemp(".xya")
1339        file = open(fileName,"w")
1340        file.write("  elevation   speed \n\
13411.0 0.0 10.0 0.0\n\
13420.0 1.0 0.0 10.0\n\
13431.0 0.0 10.4 40.0\n\
1344#geocrap\n\
134556\n\
134610\n\
134720\n")
1348        file.close()
1349       
1350        ab_points = ensure_absolute(fileName)
1351        actual =  [[11, 20.0],[10.0, 21.0],[11.0, 20.0]]
1352        assert allclose(ab_points, actual)
1353        os.remove(fileName)
1354
1355       
1356    def test_ensure_geospatial(self):
1357        points = [[2.0, 0.0],[1.0, 1.0],
1358                         [2.0, 0.0],[2.0, 2.0],
1359                         [1.0, 3.0],[2.0, 2.0]]
1360        new_points = ensure_geospatial(points)
1361       
1362        assert allclose(new_points.get_data_points(absolute = True), points)
1363
1364        points = array([[2.0, 0.0],[1.0, 1.0],
1365                         [2.0, 0.0],[2.0, 2.0],
1366                         [1.0, 3.0],[2.0, 2.0]])
1367        new_points = ensure_geospatial(points)
1368       
1369        assert allclose(new_points.get_data_points(absolute = True), points)
1370       
1371        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1372                         [2.0, 0.0],[2.0, 2.0],
1373                         [1.0, 3.0],[2.0, 2.0]])
1374       
1375        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1376
1377        data_points = zeros((ab_points.shape), Float)
1378        #Shift datapoints according to new origins
1379        for k in range(len(ab_points)):
1380            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1381            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1382        #print "data_points",data_points     
1383        new_geospatial = ensure_geospatial(data_points,
1384                                             geo_reference=mesh_origin)
1385        new_points = new_geospatial.get_data_points(absolute=True)
1386        #print "new_points",new_points
1387        #print "ab_points",ab_points
1388           
1389        assert allclose(new_points, ab_points)
1390
1391        geo = Geo_reference(56,67,-56)
1392
1393        data_points = geo.change_points_geo_ref(ab_points)   
1394        new_geospatial = ensure_geospatial(data_points,
1395                                             geo_reference=geo)
1396        new_points = new_geospatial.get_data_points(absolute=True)
1397        #print "new_points",new_points
1398        #print "ab_points",ab_points
1399           
1400        assert allclose(new_points, ab_points)
1401
1402
1403        geo_reference = Geo_reference(56, 100, 200)
1404        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1405        points = geo_reference.change_points_geo_ref(ab_points)
1406        attributes = [2, 4]
1407        #print "geo in points", points
1408        G = Geospatial_data(points, attributes,
1409                            geo_reference=geo_reference)
1410         
1411        new_geospatial  = ensure_geospatial(G)
1412        new_points = new_geospatial.get_data_points(absolute=True)
1413        #print "new_points",new_points
1414        #print "ab_points",ab_points
1415           
1416        assert allclose(new_points, ab_points)
1417       
1418    def test_isinstance(self):
1419
1420        import os
1421       
1422        fileName = tempfile.mktemp(".xya")
1423        file = open(fileName,"w")
1424        file.write("  elevation   speed \n\
14251.0 0.0 10.0 0.0\n\
14260.0 1.0 0.0 10.0\n\
14271.0 0.0 10.4 40.0\n\
1428#geocrap\n\
142956\n\
143056.6\n\
14313\n")
1432        file.close()
1433
1434        results = Geospatial_data(fileName)
1435        assert allclose(results.get_data_points(absolute=False), \
1436                        [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1437        assert allclose(results.get_attributes(attribute_name='elevation'), \
1438                        [10.0, 0.0, 10.4])
1439        assert allclose(results.get_attributes(attribute_name='speed'), \
1440                        [0.0, 10.0, 40.0])
1441
1442        os.remove(fileName)
1443       
1444    def test_delimiter(self):
1445       
1446        try:
1447            G = Geospatial_data(delimiter=',')
1448#            results = Geospatial_data(file_name = fileName)
1449#            dict = import_points_file(fileName)
1450        except ValueError:
1451            pass
1452        else:
1453            msg = 'Instance with No fileName but has a delimiter\
1454                  did not raise error!'
1455            raise msg
1456
1457    def test_no_constructors(self):
1458       
1459        try:
1460            G = Geospatial_data()
1461#            results = Geospatial_data(file_name = fileName)
1462#            dict = import_points_file(fileName)
1463        except ValueError:
1464            pass
1465        else:
1466            msg = 'Instance must have a filename or data points'
1467            raise msg       
1468
1469    def test_check_geo_reference(self):
1470        """
1471        checks geo reference details are OK. eg can be called '#geo reference'
1472        if not throws a clear error message
1473        """
1474        import os
1475        fileName = tempfile.mktemp(".xya")
1476        file = open(fileName,"w")
1477        file.write("  elevation  \n\
14781.0 0.0 10.0\n\
14790.0 1.0 0.0\n\
14801.0 0.0 10.4\n\
1481#ge oreference\n\
148256\n\
14831.1\n\
14841.0\n")
1485
1486        file.close()
1487        results = Geospatial_data(fileName)
1488        assert allclose(results.get_geo_reference().get_xllcorner(), 1.1)
1489        assert allclose(results.get_geo_reference().get_yllcorner(), 1.0)
1490
1491        os.remove(fileName)
1492       
1493        fileName = tempfile.mktemp(".xya")
1494        file = open(fileName,"w")
1495        file.write("  elevation  \n\
14961.0 0.0 10.0\n\
14970.0 1.0 0.0\n\
14981.0 0.0 10.4\n")
1499
1500        file.close()
1501        results = Geospatial_data(fileName)
1502       
1503        os.remove(fileName)
1504       
1505    def test_check_geo_reference1(self):
1506        """
1507        checks geo reference details are OK. eg can be called '#geo reference'
1508        if not throws a clear error message
1509        """
1510        import os
1511        fileName = tempfile.mktemp(".xya")
1512        file = open(fileName,"w")
1513        file.write("  elevation  \n\
15141.0 0.0 10.0\n\
15150.0 1.0 0.0\n\
15161.0 0.0 10.4\n\
1517#geo t t\n\
151856\n\
15191.1\n"
1520)
1521        file.close()
1522
1523        try:
1524            results = Geospatial_data(fileName, delimiter = " ")
1525        except IOError:
1526            pass
1527        else:
1528            msg = 'Geo reference data format is incorrect'
1529            raise msg       
1530
1531
1532        os.remove(fileName)
1533
1534
1535       
1536    def test_lat_long(self):
1537        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1538        lon_gong = degminsec2decimal_degrees(150,55,0.)
1539       
1540        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1541        lon_2 = degminsec2decimal_degrees(150,00,0.)
1542       
1543        lats = [lat_gong, lat_2]
1544        longs = [lon_gong, lon_2]
1545        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1546
1547        points = gsd.get_data_points(absolute=True)
1548       
1549        assert allclose(points[0][0], 308728.009)
1550        assert allclose(points[0][1], 6180432.601)
1551        assert allclose(points[1][0],  222908.705)
1552        assert allclose(points[1][1], 6233785.284)
1553        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1554                        'Bad zone error!')
1555       
1556        try:
1557            results = Geospatial_data(latitudes=lats)
1558        except ValueError:
1559            pass
1560        else:
1561            self.failUnless(0 ==1,  'Error not thrown error!')
1562        try:
1563            results = Geospatial_data(latitudes=lats)
1564        except ValueError:
1565            pass
1566        else:
1567            self.failUnless(0 ==1,  'Error not thrown error!')
1568        try:
1569            results = Geospatial_data(longitudes=lats)
1570        except ValueError:
1571            pass
1572        else:
1573            self.failUnless(0 ==1, 'Error not thrown error!')
1574        try:
1575            results = Geospatial_data(latitudes=lats, longitudes=longs,
1576                                      geo_reference="p")
1577        except ValueError:
1578            pass
1579        else:
1580            self.failUnless(0 ==1,  'Error not thrown error!')
1581           
1582        try:
1583            results = Geospatial_data(latitudes=lats, longitudes=longs,
1584                                      data_points=12)
1585        except ValueError:
1586            pass
1587        else:
1588            self.failUnless(0 ==1,  'Error not thrown error!')
1589
1590    def test_lat_long2(self):
1591        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1592        lon_gong = degminsec2decimal_degrees(150,55,0.)
1593       
1594        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1595        lon_2 = degminsec2decimal_degrees(150,00,0.)
1596       
1597        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1598        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1599
1600        points = gsd.get_data_points(absolute=True)
1601       
1602        assert allclose(points[0][0], 308728.009)
1603        assert allclose(points[0][1], 6180432.601)
1604        assert allclose(points[1][0],  222908.705)
1605        assert allclose(points[1][1], 6233785.284)
1606        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1607                        'Bad zone error!')
1608
1609        try:
1610            results = Geospatial_data(points_are_lats_longs=True)
1611        except ValueError:
1612            pass
1613        else:
1614            self.failUnless(0 ==1,  'Error not thrown error!')
1615
1616    def test_len(self):
1617       
1618        points = [[1.0, 2.1], [3.0, 5.3]]
1619        G = Geospatial_data(points)
1620        self.failUnless(2 ==len(G),  'Len error!')
1621       
1622        points = [[1.0, 2.1]]
1623        G = Geospatial_data(points)
1624        self.failUnless(1 ==len(G),  'Len error!')
1625
1626        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1627        G = Geospatial_data(points)
1628        self.failUnless(4 ==len(G),  'Len error!')
1629         
1630if __name__ == "__main__":
1631
1632    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_ensure_geospatial')
1633    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1634    runner = unittest.TextTestRunner()
1635    runner.run(suite)
1636
1637   
Note: See TracBrowser for help on using the repository browser.