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

Last change on this file since 3736 was 3736, checked in by sexton, 18 years ago

incorporating test to catch get_attributes from clipped geospatial data set (no_test_clip1)
File size: 54.3 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 no_test_clip1(self):
398        """test_clip1(self):
399       
400        Test that point sets can be clipped by a polygon given as
401        another Geospatial dataset
402        """
403       
404        from anuga.coordinate_transforms.geo_reference import Geo_reference
405       
406        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
407                  [0, 0], [2.4, 3.3]]
408        attributes = [2, -4, 5, 76, -2, 0.1, 3]
409        att_dict = {'att1': attributes,
410                    'att2': array(attributes) +1}
411        G = Geospatial_data(points, att_dict)
412       
413        # First try the unit square   
414        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
415        assert allclose(G.clip(U).get_data_points(),
416                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
417
418        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
419        #assert allclose(G.clip(U).get_attributes('att2'), ???)
420       
421        # Then a more complex polygon
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        G = Geospatial_data(points)
424        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
425       
426
427        assert allclose(G.clip(polygon).get_data_points(),
428                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
429
430    def test_clip0_outside(self):
431        """test_clip0_outside(self):
432       
433        Test that point sets can be clipped outside of a polygon
434        """
435       
436        from anuga.coordinate_transforms.geo_reference import Geo_reference
437       
438        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
439                  [0, 0], [2.4, 3.3]]
440        G = Geospatial_data(points)
441
442        # First try the unit square   
443        U = [[0,0], [1,0], [1,1], [0,1]]
444        assert allclose(G.clip_outside(U).get_data_points(),
445                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
446
447        # Then a more complex polygon
448        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
449        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
450        G = Geospatial_data(points)
451
452        assert allclose(G.clip_outside(polygon).get_data_points(),
453                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])   
454
455
456    def test_clip1_outside(self):
457        """test_clip1_outside(self):
458       
459        Test that point sets can be clipped outside of a polygon given as
460        another Geospatial dataset
461        """
462       
463        from anuga.coordinate_transforms.geo_reference import Geo_reference
464       
465        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
466                  [0, 0], [2.4, 3.3]]
467        G = Geospatial_data(points)
468
469        # First try the unit square   
470        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
471        assert allclose(G.clip_outside(U).get_data_points(),
472                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
473
474        # Then a more complex polygon
475        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
476        G = Geospatial_data(points)
477        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
478       
479
480        assert allclose(G.clip_outside(polygon).get_data_points(),
481                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
482
483
484    def test_clip1_inside_outside(self):
485        """test_clip1_inside_outside(self):
486       
487        Test that point sets can be clipped outside of a polygon given as
488        another Geospatial dataset
489        """
490       
491        from anuga.coordinate_transforms.geo_reference import Geo_reference
492       
493        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
494                  [0, 0], [2.4, 3.3]]
495
496        G = Geospatial_data(points)
497
498        # First try the unit square   
499        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
500        G1 = G.clip(U)
501        assert allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
502       
503        G2 = G.clip_outside(U)
504        assert allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
505                                              [3.0, 5.3], [2.4, 3.3]])
506
507       
508        # New ordering
509        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\
510                     [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
511
512        assert allclose((G1+G2).get_data_points(), new_points)
513
514        G = G1+G2
515        FN = 'test_combine.pts'
516        G.export_points_file(FN)
517
518
519        # Read it back in
520        G3 = Geospatial_data(FN)
521
522
523        # Check result
524        assert allclose(G3.get_data_points(), new_points)       
525       
526       
527        os.remove(FN)
528
529       
530    def test_create_from_xya_file(self):
531        """Check that object can be created from a points file (.pts and .xya)
532        """
533
534        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
535        attributes = [2, 4, 5, 76]
536        '''
537        # Use old pointsdict format
538        pointsdict = {'pointlist': points,
539                      'attributelist': {'att1': attributes,
540                                        'att2': array(attributes) + 1}}
541        '''
542        att_dict = {'att1': attributes,
543                    'att2': array(attributes) +1}
544                   
545        # Create points as an xya file
546        FN = 'test_points.xya'
547        G1 = Geospatial_data(points, att_dict)
548        G1.export_points_file(FN)
549#        G1.export_points_file(ofile)
550
551        #Create object from file
552        G = Geospatial_data(file_name = FN)
553       
554        assert allclose(G.get_data_points(), points)
555        assert allclose(G.get_attributes('att1'), attributes)
556        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
557       
558        os.remove(FN)
559
560    def test_create_from_xya_file1(self):
561        """
562        Check that object can be created from an Absolute xya file
563        """
564
565        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
566        attributes = [2, 4, 5, 76]
567
568        att_dict = {'att1': attributes,
569                    'att2': array(attributes) +1}
570
571        geo_ref = Geo_reference(56, 10, 5)
572                   
573        # Create points as an xya file
574        FN = 'test_points.xya'
575        G1 = Geospatial_data(points, att_dict, geo_ref)
576
577        G1.export_points_file(FN, absolute=True)
578
579        #Create object from file
580        G = Geospatial_data(file_name = FN)
581       
582        assert allclose(G.get_data_points(absolute=True), 
583                        G1.get_data_points(absolute=True))
584        assert allclose(G.get_attributes('att1'), attributes)
585        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
586       
587        os.remove(FN)
588       
589    def test_loadxya(self):
590        """
591        comma delimited
592        """
593        fileName = tempfile.mktemp(".xya")
594        file = open(fileName,"w")
595        file.write("elevation  , speed \n\
5961.0, 0.0, 10.0, 0.0\n\
5970.0, 1.0, 0.0, 10.0\n\
5981.0, 0.0, 10.4, 40.0\n")
599        file.close()
600        results = Geospatial_data(fileName, delimiter=',')
601        os.remove(fileName)
602#        print 'data', results.get_data_points()
603        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
604        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
605        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
606
607    def test_loadxya2(self):
608        """
609        space delimited
610        """
611        import os
612       
613        fileName = tempfile.mktemp(".xya")
614        file = open(fileName,"w")
615        file.write("  elevation   speed \n\
6161.0 0.0 10.0 0.0\n\
6170.0 1.0 0.0 10.0\n\
6181.0 0.0 10.4 40.0\n")
619        file.close()
620
621        results = Geospatial_data(fileName, delimiter=' ')
622
623        os.remove(fileName)
624
625        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
626        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
627        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
628     
629    def test_loadxya3(self):
630        """
631        space delimited
632        """
633        import os
634       
635        fileName = tempfile.mktemp(".xya")
636        file = open(fileName,"w")
637        file.write("  elevation   speed \n\
6381.0 0.0 10.0 0.0\n\
6390.0 1.0 0.0 10.0\n\
6401.0 0.0 10.4 40.0\n\
641#geocrap\n\
64256\n\
64356.6\n\
6443\n")
645        file.close()
646
647        results = Geospatial_data(fileName, delimiter=' ')
648
649        os.remove(fileName)
650        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
651        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
652        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
653
654    def BADtest_loadxya4(self):
655        """
656        comma delimited
657        """
658        fileName = tempfile.mktemp(".xya")
659        file = open(fileName,"w")
660        file.write("elevation  , speed \n\
6611.0, 0.0, splat, 0.0\n\
6620.0, 1.0, 0.0, 10.0\n\
6631.0, 0.0, 10.4, 40.0\n")
664        file.close()
665        results = Geospatial_data(fileName, delimiter=',')
666        os.remove(fileName)
667#        print 'data', results.get_data_points()
668        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
669        assert allclose(results.get_attributes(attribute_name='elevation'), ["splat", 0.0, 10.4])
670        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
671       
672    def test_read_write_points_file_bad2(self):
673        att_dict = {}
674        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
675        att_dict['elevation'] = array([10.0, 0.0, 10.4])
676        att_dict['brightness'] = array([10.0, 0.0, 10.4])
677        geo_reference=Geo_reference(56,1.9,1.9)
678       
679        G = Geospatial_data(pointlist, att_dict, geo_reference)
680       
681        try:
682            G.export_points_file("_???/yeah.xya")
683           
684        except IOError:
685            pass
686        else:
687            msg = 'bad points file extension did not raise error!'
688            raise msg
689#            self.failUnless(0 == 1,
690#                        'bad points file extension did not raise error!')
691                   
692    def test_loadxy_bad(self):
693        import os
694       
695        fileName = tempfile.mktemp(".xya")
696        file = open(fileName,"w")
697        file.write("  elevation   \n\
6981.0 0.0 10.0 0.0\n\
6990.0 1.0 0.0 10.0\n\
7001.0 0.0 10.4 40.0\n")
701        file.close()
702        #print fileName
703        try:
704            results = Geospatial_data(fileName, delimiter=' ')
705        except IOError:
706            pass
707        else:
708            msg = 'bad xya file did not raise error!'
709            raise msg
710#            self.failUnless(0 == 1,
711#                        'bad xya file did not raise error!')
712        os.remove(fileName)
713       
714    def test_loadxy_bad2(self):
715        import os
716       
717        fileName = tempfile.mktemp(".xya")
718        file = open(fileName,"w")
719        file.write("elevation\n\
7201.0 0.0 10.0 \n\
7210.0 1.0\n\
7221.0 \n")
723        file.close()
724        #print fileName
725        try:
726            results = Geospatial_data(fileName, delimiter=' ')
727        except IOError:
728            pass
729        else:
730            msg = 'bad xya file did not raise error!'
731            raise msg
732        os.remove(fileName)
733   
734    def test_loadxy_bad3(self):
735        """
736        specifying wrong delimiter
737        """
738        import os
739       
740        fileName = tempfile.mktemp(".xya")
741        file = open(fileName,"w")
742        file.write("  elevation  , speed \n\
7431.0, 0.0, 10.0, 0.0\n\
7440.0, 1.0, 0.0, 10.0\n\
7451.0, 0.0, 10.4, 40.0\n")
746        file.close()
747        try:
748            results = Geospatial_data(fileName, delimiter=' ')
749        except IOError:
750            pass
751        else:
752            msg = 'bad xya file did not raise error!'
753            raise msg
754        os.remove(fileName)
755     
756    def test_loadxy_bad4(self):
757        """
758         specifying wrong delimiter
759        """
760        import os
761        fileName = tempfile.mktemp(".xya")
762        file = open(fileName,"w")
763        file.write("  elevation   speed \n\
7641.0 0.0 10.0 0.0\n\
7650.0 1.0 0.0 10.0\n\
7661.0 0.0 10.4 40.0\n\
767#geocrap\n\
76856\n\
76956.6\n\
7703\n"
771)
772        file.close()
773        try:
774            results = Geospatial_data(fileName, delimiter=',')
775        except IOError:
776            pass
777        else:
778            msg = 'bad xya file did not raise error!'
779            raise msg
780
781        os.remove(fileName)
782
783    def test_loadxy_bad5(self):
784        """
785        specifying wrong delimiter
786        """
787        import os
788       
789        fileName = tempfile.mktemp(".xya")
790        file = open(fileName,"w")
791        file.write("  elevation   speed \n\
7921.0 0.0 10.0 0.0\n\
7930.0 1.0 0.0 10.0\n\
7941.0 0.0 10.4 40.0\n\
795#geocrap\n\
796crap")
797        file.close()
798        try:
799#            dict = import_points_file(fileName,delimiter=' ')
800#            results = Geospatial_data()
801            results = Geospatial_data(fileName, delimiter=' ', verbose=True)
802#            results.import_points_file(fileName, delimiter=' ')
803        except IOError:
804            pass
805        else:
806            msg = 'bad xya file did not raise error!'
807            raise msg
808
809#            self.failUnless(0 ==1,
810#                        'bad xya file did not raise error!')   
811        os.remove(fileName)             
812
813    def test_loadxy_bad_no_file_xya(self):
814        import os
815       
816        fileName = tempfile.mktemp(".xya")
817        try:
818            results = Geospatial_data(fileName, delimiter=' ')
819        except IOError:
820            pass
821        else:
822            msg = 'imaginary file did not raise error!'
823            raise msg
824
825#        except IOError:
826#            pass
827#        else:
828#            self.failUnless(0 == 1,
829#                        'imaginary file did not raise error!')
830
831                       
832  ###################### .XYA ##############################
833       
834    def test_export_xya_file(self):
835#        dict = {}
836        att_dict = {}
837        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
838        att_dict['elevation'] = array([10.0, 0.0, 10.4])
839        att_dict['brightness'] = array([10.0, 0.0, 10.4])
840#        dict['attributelist'] = att_dict
841        geo_reference=Geo_reference(56,1.9,1.9)
842       
843       
844        fileName = tempfile.mktemp(".xya")
845        G = Geospatial_data(pointlist, att_dict, geo_reference)
846        G.export_points_file(fileName, False)
847
848#        dict2 = import_points_file(fileName)
849        results = Geospatial_data(file_name = fileName)
850        #print "fileName",fileName
851        os.remove(fileName)
852       
853        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
854        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
855        answer = [10.0, 0.0, 10.4]
856        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
857        #print "dict2['geo_reference']",dict2['geo_reference']
858        self.failUnless(results.get_geo_reference() == geo_reference,
859                         'test_writepts failed. Test geo_reference')
860
861    def test_export_xya_file2(self):
862        """test absolute xya file
863        """
864        att_dict = {}
865        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
866        att_dict['elevation'] = array([10.0, 0.0, 10.4])
867        att_dict['brightness'] = array([10.0, 0.0, 10.4])
868       
869        fileName = tempfile.mktemp(".xya")
870        G = Geospatial_data(pointlist, att_dict)
871        G.export_points_file(fileName)
872        results = Geospatial_data(file_name = fileName)
873#        dict2 = import_points_file(fileName)
874        os.remove(fileName)
875       
876        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
877        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
878        answer = [10.0, 0.0, 10.4]
879        assert allclose(results.get_attributes('brightness'), answer)
880
881    def test_export_xya_file3(self):
882        """test absolute xya file with geo_ref
883        """
884        att_dict = {}
885        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
886        att_dict['elevation'] = array([10.0, 0.0, 10.4])
887        att_dict['brightness'] = array([10.0, 0.0, 10.4])
888        geo_reference=Geo_reference(56,1.9,1.9)
889       
890       
891        fileName = tempfile.mktemp(".xya")
892        G = Geospatial_data(pointlist, att_dict, geo_reference)
893       
894        G.export_points_file(fileName, absolute=True)
895       
896        results = Geospatial_data(file_name = fileName)
897        os.remove(fileName)
898
899        assert allclose(results.get_data_points(),
900                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
901        assert allclose(results.get_attributes(attribute_name='elevation'),
902                         [10.0, 0.0, 10.4])
903        answer = [10.0, 0.0, 10.4]
904        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
905        self.failUnless(results.get_geo_reference() == geo_reference,
906                         'test_writepts failed. Test geo_reference')                         
907                       
908                       
909                       
910    def test_new_export_pts_file(self):
911        att_dict = {}
912        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
913        att_dict['elevation'] = array([10.1, 0.0, 10.4])
914        att_dict['brightness'] = array([10.0, 1.0, 10.4])
915       
916        fileName = tempfile.mktemp(".pts")
917       
918        G = Geospatial_data(pointlist, att_dict)
919       
920        G.export_points_file(fileName)
921
922        results = Geospatial_data(file_name = fileName)
923
924        os.remove(fileName)
925       
926        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
927        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
928        answer = [10.0, 1.0, 10.4]
929        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
930
931    def test_new_export_absolute_pts_file(self):
932        att_dict = {}
933        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
934        att_dict['elevation'] = array([10.1, 0.0, 10.4])
935        att_dict['brightness'] = array([10.0, 1.0, 10.4])
936        geo_ref = Geo_reference(50, 25, 55)
937       
938        fileName = tempfile.mktemp(".pts")
939       
940        G = Geospatial_data(pointlist, att_dict, geo_ref)
941       
942        G.export_points_file(fileName, absolute=True)
943
944        results = Geospatial_data(file_name = fileName)
945
946        os.remove(fileName)
947       
948        assert allclose(results.get_data_points(), G.get_data_points(True))
949        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
950        answer = [10.0, 1.0, 10.4]
951        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
952
953    def test_loadpts(self):
954       
955        from Scientific.IO.NetCDF import NetCDFFile
956
957        fileName = tempfile.mktemp(".pts")
958        # NetCDF file definition
959        outfile = NetCDFFile(fileName, 'w')
960       
961        # dimension definitions
962        outfile.createDimension('number_of_points', 3)   
963        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
964   
965        # variable definitions
966        outfile.createVariable('points', Float, ('number_of_points',
967                                                 'number_of_dimensions'))
968        outfile.createVariable('elevation', Float, ('number_of_points',))
969   
970        # Get handles to the variables
971        points = outfile.variables['points']
972        elevation = outfile.variables['elevation']
973 
974        points[0, :] = [1.0,0.0]
975        elevation[0] = 10.0 
976        points[1, :] = [0.0,1.0]
977        elevation[1] = 0.0 
978        points[2, :] = [1.0,0.0]
979        elevation[2] = 10.4   
980
981        outfile.close()
982       
983        results = Geospatial_data(file_name = fileName)
984        os.remove(fileName)
985        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
986        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
987        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
988       
989    def test_writepts(self):
990        """test_writepts: Test that storage of x,y,attributes works
991        """
992        att_dict = {}
993        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
994        att_dict['elevation'] = array([10.0, 0.0, 10.4])
995        att_dict['brightness'] = array([10.0, 0.0, 10.4])
996        geo_reference=Geo_reference(56,1.9,1.9)
997
998        # Test pts format
999        fileName = tempfile.mktemp(".pts")
1000        G = Geospatial_data(pointlist, att_dict, geo_reference)
1001        G.export_points_file(fileName, False)
1002        results = Geospatial_data(file_name=fileName)
1003        os.remove(fileName)
1004
1005        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1006        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1007        answer = [10.0, 0.0, 10.4]
1008        assert allclose(results.get_attributes('brightness'), answer)
1009        self.failUnless(geo_reference == geo_reference,
1010                         'test_writepts failed. Test geo_reference')
1011
1012        # Test xya format
1013        fileName = tempfile.mktemp(".xya")
1014        G = Geospatial_data(pointlist, att_dict, geo_reference)
1015        G.export_points_file(fileName, False)
1016        results = Geospatial_data(file_name=fileName)
1017        os.remove(fileName)
1018        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1019        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1020        answer = [10.0, 0.0, 10.4]
1021        assert allclose(results.get_attributes('brightness'), answer)
1022        self.failUnless(geo_reference == geo_reference,
1023                         'test_writepts failed. Test geo_reference')
1024
1025    def test_writepts_no_attributes(self):
1026        """test_writepts_no_attributes: Test that storage of x,y alone works
1027        """
1028        att_dict = {}
1029        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1030        geo_reference=Geo_reference(56,1.9,1.9)
1031
1032        # Test pts format
1033        fileName = tempfile.mktemp(".pts")
1034        G = Geospatial_data(pointlist, None, geo_reference)
1035        G.export_points_file(fileName, False)
1036        results = Geospatial_data(file_name=fileName)
1037        os.remove(fileName)
1038
1039        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1040        self.failUnless(geo_reference == geo_reference,
1041                         'test_writepts failed. Test geo_reference')
1042
1043        # Test xya format
1044        fileName = tempfile.mktemp(".xya")
1045        G = Geospatial_data(pointlist, None, geo_reference)
1046        G.export_points_file(fileName, False)
1047        results = Geospatial_data(file_name=fileName)
1048        os.remove(fileName)
1049        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1050        self.failUnless(geo_reference == geo_reference,
1051                         'test_writepts failed. Test geo_reference')
1052
1053       
1054       
1055 ########################## BAD .PTS ##########################         
1056
1057    def test_load_bad_no_file_pts(self):
1058        import os
1059        import tempfile
1060       
1061        fileName = tempfile.mktemp(".pts")
1062        #print fileName
1063        try:
1064            results = Geospatial_data(file_name = fileName)
1065#            dict = import_points_file(fileName)
1066        except IOError:
1067            pass
1068        else:
1069            msg = 'imaginary file did not raise error!'
1070            raise msg
1071#            self.failUnless(0 == 1,
1072#                        'imaginary file did not raise error!')
1073
1074
1075    def test_create_from_pts_file(self):
1076       
1077        from Scientific.IO.NetCDF import NetCDFFile
1078
1079#        fileName = tempfile.mktemp(".pts")
1080        FN = 'test_points.pts'
1081        # NetCDF file definition
1082        outfile = NetCDFFile(FN, 'w')
1083       
1084        # dimension definitions
1085        outfile.createDimension('number_of_points', 3)   
1086        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1087   
1088        # variable definitions
1089        outfile.createVariable('points', Float, ('number_of_points',
1090                                                 'number_of_dimensions'))
1091        outfile.createVariable('elevation', Float, ('number_of_points',))
1092   
1093        # Get handles to the variables
1094        points = outfile.variables['points']
1095        elevation = outfile.variables['elevation']
1096 
1097        points[0, :] = [1.0,0.0]
1098        elevation[0] = 10.0 
1099        points[1, :] = [0.0,1.0]
1100        elevation[1] = 0.0 
1101        points[2, :] = [1.0,0.0]
1102        elevation[2] = 10.4   
1103
1104        outfile.close()
1105
1106        G = Geospatial_data(file_name = FN)
1107
1108        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
1109        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
1110
1111        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1112        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1113        os.remove(FN)
1114
1115    def test_create_from_pts_file_with_geo(self):
1116        """This test reveals if Geospatial data is correctly instantiated from a pts file.
1117        """
1118       
1119        from Scientific.IO.NetCDF import NetCDFFile
1120
1121        FN = 'test_points.pts'
1122        # NetCDF file definition
1123        outfile = NetCDFFile(FN, 'w')
1124
1125        # Make up an arbitrary georef
1126        xll = 0.1
1127        yll = 20
1128        geo_reference=Geo_reference(56, xll, yll)
1129        geo_reference.write_NetCDF(outfile)
1130
1131        # dimension definitions
1132        outfile.createDimension('number_of_points', 3)   
1133        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1134   
1135        # variable definitions
1136        outfile.createVariable('points', Float, ('number_of_points',
1137                                                 'number_of_dimensions'))
1138        outfile.createVariable('elevation', Float, ('number_of_points',))
1139   
1140        # Get handles to the variables
1141        points = outfile.variables['points']
1142        elevation = outfile.variables['elevation']
1143
1144        points[0, :] = [1.0,0.0]
1145        elevation[0] = 10.0 
1146        points[1, :] = [0.0,1.0]
1147        elevation[1] = 0.0 
1148        points[2, :] = [1.0,0.0]
1149        elevation[2] = 10.4   
1150
1151        outfile.close()
1152
1153        G = Geospatial_data(file_name = FN)
1154
1155        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
1156        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
1157
1158        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
1159                                              [0.0+xll, 1.0+yll],
1160                                              [1.0+xll, 0.0+yll]])
1161       
1162        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1163        os.remove(FN)
1164
1165       
1166    def test_add_(self):
1167        '''adds an xya and pts files, reads the files and adds them
1168           checking results are correct
1169        '''
1170
1171        # create files
1172        att_dict1 = {}
1173        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1174        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
1175        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
1176        geo_reference1 = Geo_reference(56, 2.0, 1.0)
1177       
1178        att_dict2 = {}
1179        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
1180        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
1181        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
1182        geo_reference2 = Geo_reference(56, 1.0, 2.0) 
1183
1184        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1185        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1186       
1187        fileName1 = tempfile.mktemp(".xya")
1188        fileName2 = tempfile.mktemp(".pts")
1189
1190        #makes files
1191        G1.export_points_file(fileName1)
1192        G2.export_points_file(fileName2)
1193       
1194        # add files
1195       
1196        G3 = Geospatial_data(file_name = fileName1)
1197        G4 = Geospatial_data(file_name = fileName2)
1198       
1199        G = G3 + G4
1200
1201       
1202        #read results
1203#        print'res', G.get_data_points()
1204#        print'res1', G.get_data_points(False)
1205        assert allclose(G.get_data_points(),
1206                        [[ 3.0, 1.0], [ 2.0, 2.0],
1207                         [ 3.0, 1.0], [ 3.0, 3.0],
1208                         [ 2.0, 4.0], [ 3.0, 3.0]])
1209                         
1210        assert allclose(G.get_attributes(attribute_name='elevation'),
1211                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
1212       
1213        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
1214        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
1215       
1216        self.failUnless(G.get_geo_reference() == geo_reference1,
1217                         'test_writepts failed. Test geo_reference')
1218                         
1219        os.remove(fileName1)
1220        os.remove(fileName2)
1221       
1222    def test_ensure_absolute(self):
1223        points = [[2.0, 0.0],[1.0, 1.0],
1224                         [2.0, 0.0],[2.0, 2.0],
1225                         [1.0, 3.0],[2.0, 2.0]]
1226        new_points = ensure_absolute(points)
1227       
1228        assert allclose(new_points, points)
1229
1230        points = array([[2.0, 0.0],[1.0, 1.0],
1231                         [2.0, 0.0],[2.0, 2.0],
1232                         [1.0, 3.0],[2.0, 2.0]])
1233        new_points = ensure_absolute(points)
1234       
1235        assert allclose(new_points, points)
1236       
1237        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1238                         [2.0, 0.0],[2.0, 2.0],
1239                         [1.0, 3.0],[2.0, 2.0]])
1240       
1241        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1242
1243        data_points = zeros((ab_points.shape), Float)
1244        #Shift datapoints according to new origins
1245        for k in range(len(ab_points)):
1246            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1247            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1248        #print "data_points",data_points     
1249        new_points = ensure_absolute(data_points,
1250                                             geo_reference=mesh_origin)
1251        #print "new_points",new_points
1252        #print "ab_points",ab_points
1253           
1254        assert allclose(new_points, ab_points)
1255
1256        geo = Geo_reference(56,67,-56)
1257
1258        data_points = geo.change_points_geo_ref(ab_points)   
1259        new_points = ensure_absolute(data_points,
1260                                             geo_reference=geo)
1261        #print "new_points",new_points
1262        #print "ab_points",ab_points
1263           
1264        assert allclose(new_points, ab_points)
1265
1266
1267        geo_reference = Geo_reference(56, 100, 200)
1268        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1269        points = geo_reference.change_points_geo_ref(ab_points)
1270        attributes = [2, 4]
1271        #print "geo in points", points
1272        G = Geospatial_data(points, attributes,
1273                            geo_reference=geo_reference)
1274         
1275        new_points = ensure_absolute(G)
1276        #print "new_points",new_points
1277        #print "ab_points",ab_points
1278           
1279        assert allclose(new_points, ab_points)
1280
1281       
1282        fileName = tempfile.mktemp(".xya")
1283        file = open(fileName,"w")
1284        file.write("  elevation   speed \n\
12851.0 0.0 10.0 0.0\n\
12860.0 1.0 0.0 10.0\n\
12871.0 0.0 10.4 40.0\n\
1288#geocrap\n\
128956\n\
129010\n\
129120\n")
1292        file.close()
1293       
1294        ab_points = ensure_absolute(fileName)
1295        actual =  [[11, 20.0],[10.0, 21.0],[11.0, 20.0]]
1296        assert allclose(ab_points, actual)
1297        os.remove(fileName)
1298
1299       
1300    def test_ensure_geospatial(self):
1301        points = [[2.0, 0.0],[1.0, 1.0],
1302                         [2.0, 0.0],[2.0, 2.0],
1303                         [1.0, 3.0],[2.0, 2.0]]
1304        new_points = ensure_geospatial(points)
1305       
1306        assert allclose(new_points.get_data_points(absolute = True), points)
1307
1308        points = array([[2.0, 0.0],[1.0, 1.0],
1309                         [2.0, 0.0],[2.0, 2.0],
1310                         [1.0, 3.0],[2.0, 2.0]])
1311        new_points = ensure_geospatial(points)
1312       
1313        assert allclose(new_points.get_data_points(absolute = True), points)
1314       
1315        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1316                         [2.0, 0.0],[2.0, 2.0],
1317                         [1.0, 3.0],[2.0, 2.0]])
1318       
1319        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1320
1321        data_points = zeros((ab_points.shape), Float)
1322        #Shift datapoints according to new origins
1323        for k in range(len(ab_points)):
1324            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1325            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1326        #print "data_points",data_points     
1327        new_geospatial = ensure_geospatial(data_points,
1328                                             geo_reference=mesh_origin)
1329        new_points = new_geospatial.get_data_points(absolute=True)
1330        #print "new_points",new_points
1331        #print "ab_points",ab_points
1332           
1333        assert allclose(new_points, ab_points)
1334
1335        geo = Geo_reference(56,67,-56)
1336
1337        data_points = geo.change_points_geo_ref(ab_points)   
1338        new_geospatial = ensure_geospatial(data_points,
1339                                             geo_reference=geo)
1340        new_points = new_geospatial.get_data_points(absolute=True)
1341        #print "new_points",new_points
1342        #print "ab_points",ab_points
1343           
1344        assert allclose(new_points, ab_points)
1345
1346
1347        geo_reference = Geo_reference(56, 100, 200)
1348        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1349        points = geo_reference.change_points_geo_ref(ab_points)
1350        attributes = [2, 4]
1351        #print "geo in points", points
1352        G = Geospatial_data(points, attributes,
1353                            geo_reference=geo_reference)
1354         
1355        new_geospatial  = ensure_geospatial(G)
1356        new_points = new_geospatial.get_data_points(absolute=True)
1357        #print "new_points",new_points
1358        #print "ab_points",ab_points
1359           
1360        assert allclose(new_points, ab_points)
1361       
1362    def test_isinstance(self):
1363
1364        import os
1365       
1366        fileName = tempfile.mktemp(".xya")
1367        file = open(fileName,"w")
1368        file.write("  elevation   speed \n\
13691.0 0.0 10.0 0.0\n\
13700.0 1.0 0.0 10.0\n\
13711.0 0.0 10.4 40.0\n\
1372#geocrap\n\
137356\n\
137456.6\n\
13753\n")
1376        file.close()
1377
1378        results = Geospatial_data(fileName)
1379        assert allclose(results.get_data_points(absolute=False), \
1380                        [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1381        assert allclose(results.get_attributes(attribute_name='elevation'), \
1382                        [10.0, 0.0, 10.4])
1383        assert allclose(results.get_attributes(attribute_name='speed'), \
1384                        [0.0, 10.0, 40.0])
1385
1386        os.remove(fileName)
1387       
1388    def test_delimiter(self):
1389       
1390        try:
1391            G = Geospatial_data(delimiter=',')
1392#            results = Geospatial_data(file_name = fileName)
1393#            dict = import_points_file(fileName)
1394        except ValueError:
1395            pass
1396        else:
1397            msg = 'Instance with No fileName but has a delimiter\
1398                  did not raise error!'
1399            raise msg
1400
1401    def test_no_constructors(self):
1402       
1403        try:
1404            G = Geospatial_data()
1405#            results = Geospatial_data(file_name = fileName)
1406#            dict = import_points_file(fileName)
1407        except ValueError:
1408            pass
1409        else:
1410            msg = 'Instance must have a filename or data points'
1411            raise msg       
1412
1413    def test_check_geo_reference(self):
1414        """
1415        checks geo reference details are OK. eg can be called '#geo reference'
1416        if not throws a clear error message
1417        """
1418        import os
1419        fileName = tempfile.mktemp(".xya")
1420        file = open(fileName,"w")
1421        file.write("  elevation  \n\
14221.0 0.0 10.0\n\
14230.0 1.0 0.0\n\
14241.0 0.0 10.4\n\
1425#ge oreference\n\
142656\n\
14271.1\n\
14281.0\n")
1429
1430        file.close()
1431        results = Geospatial_data(fileName)
1432        assert allclose(results.get_geo_reference().get_xllcorner(), 1.1)
1433        assert allclose(results.get_geo_reference().get_yllcorner(), 1.0)
1434
1435        os.remove(fileName)
1436       
1437        fileName = tempfile.mktemp(".xya")
1438        file = open(fileName,"w")
1439        file.write("  elevation  \n\
14401.0 0.0 10.0\n\
14410.0 1.0 0.0\n\
14421.0 0.0 10.4\n")
1443
1444        file.close()
1445        results = Geospatial_data(fileName)
1446       
1447        os.remove(fileName)
1448       
1449    def test_check_geo_reference1(self):
1450        """
1451        checks geo reference details are OK. eg can be called '#geo reference'
1452        if not throws a clear error message
1453        """
1454        import os
1455        fileName = tempfile.mktemp(".xya")
1456        file = open(fileName,"w")
1457        file.write("  elevation  \n\
14581.0 0.0 10.0\n\
14590.0 1.0 0.0\n\
14601.0 0.0 10.4\n\
1461#geo t t\n\
146256\n\
14631.1\n"
1464)
1465        file.close()
1466
1467        try:
1468            results = Geospatial_data(fileName, delimiter = " ")
1469        except IOError:
1470            pass
1471        else:
1472            msg = 'Geo reference data format is incorrect'
1473            raise msg       
1474
1475
1476        os.remove(fileName)
1477
1478
1479       
1480    def test_lat_long(self):
1481        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1482        lon_gong = degminsec2decimal_degrees(150,55,0.)
1483       
1484        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1485        lon_2 = degminsec2decimal_degrees(150,00,0.)
1486       
1487        lats = [lat_gong, lat_2]
1488        longs = [lon_gong, lon_2]
1489        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1490
1491        points = gsd.get_data_points(absolute=True)
1492       
1493        assert allclose(points[0][0], 308728.009)
1494        assert allclose(points[0][1], 6180432.601)
1495        assert allclose(points[1][0],  222908.705)
1496        assert allclose(points[1][1], 6233785.284)
1497        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1498                        'Bad zone error!')
1499       
1500        try:
1501            results = Geospatial_data(latitudes=lats)
1502        except ValueError:
1503            pass
1504        else:
1505            self.failUnless(0 ==1,  'Error not thrown error!')
1506        try:
1507            results = Geospatial_data(latitudes=lats)
1508        except ValueError:
1509            pass
1510        else:
1511            self.failUnless(0 ==1,  'Error not thrown error!')
1512        try:
1513            results = Geospatial_data(longitudes=lats)
1514        except ValueError:
1515            pass
1516        else:
1517            self.failUnless(0 ==1, 'Error not thrown error!')
1518        try:
1519            results = Geospatial_data(latitudes=lats, longitudes=longs,
1520                                      geo_reference="p")
1521        except ValueError:
1522            pass
1523        else:
1524            self.failUnless(0 ==1,  'Error not thrown error!')
1525           
1526        try:
1527            results = Geospatial_data(latitudes=lats, longitudes=longs,
1528                                      data_points=12)
1529        except ValueError:
1530            pass
1531        else:
1532            self.failUnless(0 ==1,  'Error not thrown error!')
1533
1534    def test_lat_long2(self):
1535        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1536        lon_gong = degminsec2decimal_degrees(150,55,0.)
1537       
1538        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1539        lon_2 = degminsec2decimal_degrees(150,00,0.)
1540       
1541        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1542        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1543
1544        points = gsd.get_data_points(absolute=True)
1545       
1546        assert allclose(points[0][0], 308728.009)
1547        assert allclose(points[0][1], 6180432.601)
1548        assert allclose(points[1][0],  222908.705)
1549        assert allclose(points[1][1], 6233785.284)
1550        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1551                        'Bad zone error!')
1552
1553        try:
1554            results = Geospatial_data(points_are_lats_longs=True)
1555        except ValueError:
1556            pass
1557        else:
1558            self.failUnless(0 ==1,  'Error not thrown error!')
1559
1560    def test_len(self):
1561       
1562        points = [[1.0, 2.1], [3.0, 5.3]]
1563        G = Geospatial_data(points)
1564        self.failUnless(2 ==len(G),  'Len error!')
1565       
1566        points = [[1.0, 2.1]]
1567        G = Geospatial_data(points)
1568        self.failUnless(1 ==len(G),  'Len error!')
1569
1570        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1571        G = Geospatial_data(points)
1572        self.failUnless(4 ==len(G),  'Len error!')
1573         
1574if __name__ == "__main__":
1575
1576    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_ensure_geospatial')
1577    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1578    runner = unittest.TextTestRunner()
1579    runner.run(suite)
1580
1581   
Note: See TracBrowser for help on using the repository browser.