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

Last change on this file since 3702 was 3702, checked in by ole, 18 years ago

Added clipping method to Geospatial data.

File size: 48.6 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_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        G = Geospatial_data(points)
409
410        # First try the unit square   
411        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
412        assert allclose(G.clip(U).get_data_points(),
413                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
414
415        # Then a more complex polygon
416        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
417        G = Geospatial_data(points)
418        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
419       
420
421        assert allclose(G.clip(polygon).get_data_points(),
422                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
423
424
425       
426
427
428    def test_create_from_xya_file(self):
429        """Check that object can be created from a points file (.pts and .xya)
430        """
431
432        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
433        attributes = [2, 4, 5, 76]
434        '''
435        # Use old pointsdict format
436        pointsdict = {'pointlist': points,
437                      'attributelist': {'att1': attributes,
438                                        'att2': array(attributes) + 1}}
439        '''
440        att_dict = {'att1': attributes,
441                    'att2': array(attributes) +1}
442                   
443        # Create points as an xya file
444        FN = 'test_points.xya'
445        G1 = Geospatial_data(points, att_dict)
446        G1.export_points_file(FN)
447#        G1.export_points_file(ofile)
448
449        #Create object from file
450        G = Geospatial_data(file_name = FN)
451       
452        assert allclose(G.get_data_points(), points)
453        assert allclose(G.get_attributes('att1'), attributes)
454        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
455       
456        os.remove(FN)
457
458    def test_create_from_xya_file1(self):
459        """
460        Check that object can be created from an Absolute xya file
461        """
462
463        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
464        attributes = [2, 4, 5, 76]
465
466        att_dict = {'att1': attributes,
467                    'att2': array(attributes) +1}
468
469        geo_ref = Geo_reference(56, 10, 5)
470                   
471        # Create points as an xya file
472        FN = 'test_points.xya'
473        G1 = Geospatial_data(points, att_dict, geo_ref)
474
475        G1.export_points_file(FN, absolute=True)
476
477        #Create object from file
478        G = Geospatial_data(file_name = FN)
479       
480        assert allclose(G.get_data_points(absolute=True), 
481                        G1.get_data_points(absolute=True))
482        assert allclose(G.get_attributes('att1'), attributes)
483        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
484       
485        os.remove(FN)
486       
487    def test_loadxya(self):
488        """
489        comma delimited
490        """
491        fileName = tempfile.mktemp(".xya")
492        file = open(fileName,"w")
493        file.write("elevation  , speed \n\
4941.0, 0.0, 10.0, 0.0\n\
4950.0, 1.0, 0.0, 10.0\n\
4961.0, 0.0, 10.4, 40.0\n")
497        file.close()
498        results = Geospatial_data(fileName, delimiter=',')
499        os.remove(fileName)
500#        print 'data', results.get_data_points()
501        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
502        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
503        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
504
505    def test_loadxya2(self):
506        """
507        space delimited
508        """
509        import os
510       
511        fileName = tempfile.mktemp(".xya")
512        file = open(fileName,"w")
513        file.write("  elevation   speed \n\
5141.0 0.0 10.0 0.0\n\
5150.0 1.0 0.0 10.0\n\
5161.0 0.0 10.4 40.0\n")
517        file.close()
518
519        results = Geospatial_data(fileName, delimiter=' ')
520
521        os.remove(fileName)
522
523        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
524        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
525        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
526     
527    def test_loadxya3(self):
528        """
529        space delimited
530        """
531        import os
532       
533        fileName = tempfile.mktemp(".xya")
534        file = open(fileName,"w")
535        file.write("  elevation   speed \n\
5361.0 0.0 10.0 0.0\n\
5370.0 1.0 0.0 10.0\n\
5381.0 0.0 10.4 40.0\n\
539#geocrap\n\
54056\n\
54156.6\n\
5423\n")
543        file.close()
544
545        results = Geospatial_data(fileName, delimiter=' ')
546
547        os.remove(fileName)
548        assert allclose(results.get_data_points(absolute=False), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
549        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
550        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
551
552    def BADtest_loadxya4(self):
553        """
554        comma delimited
555        """
556        fileName = tempfile.mktemp(".xya")
557        file = open(fileName,"w")
558        file.write("elevation  , speed \n\
5591.0, 0.0, splat, 0.0\n\
5600.0, 1.0, 0.0, 10.0\n\
5611.0, 0.0, 10.4, 40.0\n")
562        file.close()
563        results = Geospatial_data(fileName, delimiter=',')
564        os.remove(fileName)
565#        print 'data', results.get_data_points()
566        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
567        assert allclose(results.get_attributes(attribute_name='elevation'), ["splat", 0.0, 10.4])
568        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
569       
570    def test_read_write_points_file_bad2(self):
571        att_dict = {}
572        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
573        att_dict['elevation'] = array([10.0, 0.0, 10.4])
574        att_dict['brightness'] = array([10.0, 0.0, 10.4])
575        geo_reference=Geo_reference(56,1.9,1.9)
576       
577        G = Geospatial_data(pointlist, att_dict, geo_reference)
578       
579        try:
580            G.export_points_file("_???/yeah.xya")
581           
582        except IOError:
583            pass
584        else:
585            msg = 'bad points file extension did not raise error!'
586            raise msg
587#            self.failUnless(0 == 1,
588#                        'bad points file extension did not raise error!')
589                   
590    def test_loadxy_bad(self):
591        import os
592       
593        fileName = tempfile.mktemp(".xya")
594        file = open(fileName,"w")
595        file.write("  elevation   \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        #print fileName
601        try:
602            results = Geospatial_data(fileName, delimiter=' ')
603        except IOError:
604            pass
605        else:
606            msg = 'bad xya file did not raise error!'
607            raise msg
608#            self.failUnless(0 == 1,
609#                        'bad xya file did not raise error!')
610        os.remove(fileName)
611       
612    def test_loadxy_bad2(self):
613        import os
614       
615        fileName = tempfile.mktemp(".xya")
616        file = open(fileName,"w")
617        file.write("elevation\n\
6181.0 0.0 10.0 \n\
6190.0 1.0\n\
6201.0 \n")
621        file.close()
622        #print fileName
623        try:
624            results = Geospatial_data(fileName, delimiter=' ')
625        except IOError:
626            pass
627        else:
628            msg = 'bad xya file did not raise error!'
629            raise msg
630        os.remove(fileName)
631   
632    def test_loadxy_bad3(self):
633        """
634        specifying wrong delimiter
635        """
636        import os
637       
638        fileName = tempfile.mktemp(".xya")
639        file = open(fileName,"w")
640        file.write("  elevation  , speed \n\
6411.0, 0.0, 10.0, 0.0\n\
6420.0, 1.0, 0.0, 10.0\n\
6431.0, 0.0, 10.4, 40.0\n")
644        file.close()
645        try:
646            results = Geospatial_data(fileName, delimiter=' ')
647        except IOError:
648            pass
649        else:
650            msg = 'bad xya file did not raise error!'
651            raise msg
652        os.remove(fileName)
653     
654    def test_loadxy_bad4(self):
655        """
656         specifying wrong delimiter
657        """
658        import os
659        fileName = tempfile.mktemp(".xya")
660        file = open(fileName,"w")
661        file.write("  elevation   speed \n\
6621.0 0.0 10.0 0.0\n\
6630.0 1.0 0.0 10.0\n\
6641.0 0.0 10.4 40.0\n\
665#geocrap\n\
66656\n\
66756.6\n\
6683\n"
669)
670        file.close()
671        try:
672            results = Geospatial_data(fileName, delimiter=',')
673        except IOError:
674            pass
675        else:
676            msg = 'bad xya file did not raise error!'
677            raise msg
678
679        os.remove(fileName)
680
681    def test_loadxy_bad5(self):
682        """
683        specifying wrong delimiter
684        """
685        import os
686       
687        fileName = tempfile.mktemp(".xya")
688        file = open(fileName,"w")
689        file.write("  elevation   speed \n\
6901.0 0.0 10.0 0.0\n\
6910.0 1.0 0.0 10.0\n\
6921.0 0.0 10.4 40.0\n\
693#geocrap\n\
694crap")
695        file.close()
696        try:
697#            dict = import_points_file(fileName,delimiter=' ')
698#            results = Geospatial_data()
699            results = Geospatial_data(fileName, delimiter=' ', verbose=True)
700#            results.import_points_file(fileName, delimiter=' ')
701        except IOError:
702            pass
703        else:
704            msg = 'bad xya file did not raise error!'
705            raise msg
706
707#            self.failUnless(0 ==1,
708#                        'bad xya file did not raise error!')   
709        os.remove(fileName)             
710
711    def test_loadxy_bad_no_file_xya(self):
712        import os
713       
714        fileName = tempfile.mktemp(".xya")
715        try:
716            results = Geospatial_data(fileName, delimiter=' ')
717        except IOError:
718            pass
719        else:
720            msg = 'imaginary file did not raise error!'
721            raise msg
722
723#        except IOError:
724#            pass
725#        else:
726#            self.failUnless(0 == 1,
727#                        'imaginary file did not raise error!')
728
729                       
730  ###################### .XYA ##############################
731       
732    def test_export_xya_file(self):
733#        dict = {}
734        att_dict = {}
735        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
736        att_dict['elevation'] = array([10.0, 0.0, 10.4])
737        att_dict['brightness'] = array([10.0, 0.0, 10.4])
738#        dict['attributelist'] = att_dict
739        geo_reference=Geo_reference(56,1.9,1.9)
740       
741       
742        fileName = tempfile.mktemp(".xya")
743        G = Geospatial_data(pointlist, att_dict, geo_reference)
744        G.export_points_file(fileName, False)
745
746#        dict2 = import_points_file(fileName)
747        results = Geospatial_data(file_name = fileName)
748        #print "fileName",fileName
749        os.remove(fileName)
750       
751        assert allclose(results.get_data_points(absolute=False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
752        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
753        answer = [10.0, 0.0, 10.4]
754        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
755        #print "dict2['geo_reference']",dict2['geo_reference']
756        self.failUnless(results.get_geo_reference() == geo_reference,
757                         'test_writepts failed. Test geo_reference')
758
759    def test_export_xya_file2(self):
760        """test absolute xya file
761        """
762        att_dict = {}
763        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
764        att_dict['elevation'] = array([10.0, 0.0, 10.4])
765        att_dict['brightness'] = array([10.0, 0.0, 10.4])
766       
767        fileName = tempfile.mktemp(".xya")
768        G = Geospatial_data(pointlist, att_dict)
769        G.export_points_file(fileName)
770        results = Geospatial_data(file_name = fileName)
771#        dict2 = import_points_file(fileName)
772        os.remove(fileName)
773       
774        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
775        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
776        answer = [10.0, 0.0, 10.4]
777        assert allclose(results.get_attributes('brightness'), answer)
778
779    def test_export_xya_file3(self):
780        """test absolute xya file with geo_ref
781        """
782        att_dict = {}
783        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
784        att_dict['elevation'] = array([10.0, 0.0, 10.4])
785        att_dict['brightness'] = array([10.0, 0.0, 10.4])
786        geo_reference=Geo_reference(56,1.9,1.9)
787       
788       
789        fileName = tempfile.mktemp(".xya")
790        G = Geospatial_data(pointlist, att_dict, geo_reference)
791       
792        G.export_points_file(fileName, absolute=True)
793       
794        results = Geospatial_data(file_name = fileName)
795        os.remove(fileName)
796
797        assert allclose(results.get_data_points(),
798                        [[2.9, 1.9],[1.9, 2.9],[2.9, 1.9]])
799        assert allclose(results.get_attributes(attribute_name='elevation'),
800                         [10.0, 0.0, 10.4])
801        answer = [10.0, 0.0, 10.4]
802        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
803        self.failUnless(results.get_geo_reference() == geo_reference,
804                         'test_writepts failed. Test geo_reference')                         
805                       
806                       
807                       
808    def test_new_export_pts_file(self):
809        att_dict = {}
810        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
811        att_dict['elevation'] = array([10.1, 0.0, 10.4])
812        att_dict['brightness'] = array([10.0, 1.0, 10.4])
813       
814        fileName = tempfile.mktemp(".pts")
815       
816        G = Geospatial_data(pointlist, att_dict)
817       
818        G.export_points_file(fileName)
819
820        results = Geospatial_data(file_name = fileName)
821
822        os.remove(fileName)
823       
824        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
825        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
826        answer = [10.0, 1.0, 10.4]
827        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
828
829    def test_new_export_absolute_pts_file(self):
830        att_dict = {}
831        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
832        att_dict['elevation'] = array([10.1, 0.0, 10.4])
833        att_dict['brightness'] = array([10.0, 1.0, 10.4])
834        geo_ref = Geo_reference(50, 25, 55)
835       
836        fileName = tempfile.mktemp(".pts")
837       
838        G = Geospatial_data(pointlist, att_dict, geo_ref)
839       
840        G.export_points_file(fileName, absolute=True)
841
842        results = Geospatial_data(file_name = fileName)
843
844        os.remove(fileName)
845       
846        assert allclose(results.get_data_points(), G.get_data_points(True))
847        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
848        answer = [10.0, 1.0, 10.4]
849        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
850
851    def test_loadpts(self):
852       
853        from Scientific.IO.NetCDF import NetCDFFile
854
855        fileName = tempfile.mktemp(".pts")
856        # NetCDF file definition
857        outfile = NetCDFFile(fileName, 'w')
858       
859        # dimension definitions
860        outfile.createDimension('number_of_points', 3)   
861        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
862   
863        # variable definitions
864        outfile.createVariable('points', Float, ('number_of_points',
865                                                 'number_of_dimensions'))
866        outfile.createVariable('elevation', Float, ('number_of_points',))
867   
868        # Get handles to the variables
869        points = outfile.variables['points']
870        elevation = outfile.variables['elevation']
871 
872        points[0, :] = [1.0,0.0]
873        elevation[0] = 10.0 
874        points[1, :] = [0.0,1.0]
875        elevation[1] = 0.0 
876        points[2, :] = [1.0,0.0]
877        elevation[2] = 10.4   
878
879        outfile.close()
880       
881        results = Geospatial_data(file_name = fileName)
882        os.remove(fileName)
883        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
884        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
885        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
886       
887    def test_writepts(self):
888        att_dict = {}
889        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
890        att_dict['elevation'] = array([10.0, 0.0, 10.4])
891        att_dict['brightness'] = array([10.0, 0.0, 10.4])
892        geo_reference=Geo_reference(56,1.9,1.9)
893       
894        fileName = tempfile.mktemp(".pts")
895       
896        G = Geospatial_data(pointlist, att_dict, geo_reference)
897       
898        G.export_points_file(fileName, False)
899       
900        results = Geospatial_data(file_name = fileName)
901
902        os.remove(fileName)
903
904        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
905        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
906        answer = [10.0, 0.0, 10.4]
907        assert allclose(results.get_attributes('brightness'), answer)
908
909       
910        self.failUnless(geo_reference == geo_reference,
911                         'test_writepts failed. Test geo_reference')
912       
913 ########################## BAD .PTS ##########################         
914
915    def test_load_bad_no_file_pts(self):
916        import os
917        import tempfile
918       
919        fileName = tempfile.mktemp(".pts")
920        #print fileName
921        try:
922            results = Geospatial_data(file_name = fileName)
923#            dict = import_points_file(fileName)
924        except IOError:
925            pass
926        else:
927            msg = 'imaginary file did not raise error!'
928            raise msg
929#            self.failUnless(0 == 1,
930#                        'imaginary file did not raise error!')
931
932
933    def test_create_from_pts_file(self):
934       
935        from Scientific.IO.NetCDF import NetCDFFile
936
937#        fileName = tempfile.mktemp(".pts")
938        FN = 'test_points.pts'
939        # NetCDF file definition
940        outfile = NetCDFFile(FN, 'w')
941       
942        # dimension definitions
943        outfile.createDimension('number_of_points', 3)   
944        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
945   
946        # variable definitions
947        outfile.createVariable('points', Float, ('number_of_points',
948                                                 'number_of_dimensions'))
949        outfile.createVariable('elevation', Float, ('number_of_points',))
950   
951        # Get handles to the variables
952        points = outfile.variables['points']
953        elevation = outfile.variables['elevation']
954 
955        points[0, :] = [1.0,0.0]
956        elevation[0] = 10.0 
957        points[1, :] = [0.0,1.0]
958        elevation[1] = 0.0 
959        points[2, :] = [1.0,0.0]
960        elevation[2] = 10.4   
961
962        outfile.close()
963
964        G = Geospatial_data(file_name = FN)
965
966        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
967        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
968
969        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
970        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
971        os.remove(FN)
972
973    def test_create_from_pts_file_with_geo(self):
974        """This test reveals if Geospatial data is correctly instantiated from a pts file.
975        """
976       
977        from Scientific.IO.NetCDF import NetCDFFile
978
979        FN = 'test_points.pts'
980        # NetCDF file definition
981        outfile = NetCDFFile(FN, 'w')
982
983        # Make up an arbitrary georef
984        xll = 0.1
985        yll = 20
986        geo_reference=Geo_reference(56, xll, yll)
987        geo_reference.write_NetCDF(outfile)
988
989        # dimension definitions
990        outfile.createDimension('number_of_points', 3)   
991        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
992   
993        # variable definitions
994        outfile.createVariable('points', Float, ('number_of_points',
995                                                 'number_of_dimensions'))
996        outfile.createVariable('elevation', Float, ('number_of_points',))
997   
998        # Get handles to the variables
999        points = outfile.variables['points']
1000        elevation = outfile.variables['elevation']
1001
1002        points[0, :] = [1.0,0.0]
1003        elevation[0] = 10.0 
1004        points[1, :] = [0.0,1.0]
1005        elevation[1] = 0.0 
1006        points[2, :] = [1.0,0.0]
1007        elevation[2] = 10.4   
1008
1009        outfile.close()
1010
1011        G = Geospatial_data(file_name = FN)
1012
1013        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
1014        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
1015
1016        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
1017                                              [0.0+xll, 1.0+yll],
1018                                              [1.0+xll, 0.0+yll]])
1019       
1020        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1021        os.remove(FN)
1022
1023       
1024    def test_add_(self):
1025        '''adds an xya and pts files, reads the files and adds them
1026           checking results are correct
1027        '''
1028
1029        # create files
1030        att_dict1 = {}
1031        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1032        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
1033        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
1034        geo_reference1 = Geo_reference(56, 2.0, 1.0)
1035       
1036        att_dict2 = {}
1037        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
1038        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
1039        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
1040        geo_reference2 = Geo_reference(56, 1.0, 2.0) 
1041
1042        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1043        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1044       
1045        fileName1 = tempfile.mktemp(".xya")
1046        fileName2 = tempfile.mktemp(".pts")
1047
1048        #makes files
1049        G1.export_points_file(fileName1)
1050        G2.export_points_file(fileName2)
1051       
1052        # add files
1053       
1054        G3 = Geospatial_data(file_name = fileName1)
1055        G4 = Geospatial_data(file_name = fileName2)
1056       
1057        G = G3 + G4
1058
1059       
1060        #read results
1061#        print'res', G.get_data_points()
1062#        print'res1', G.get_data_points(False)
1063        assert allclose(G.get_data_points(),
1064                        [[ 3.0, 1.0], [ 2.0, 2.0],
1065                         [ 3.0, 1.0], [ 3.0, 3.0],
1066                         [ 2.0, 4.0], [ 3.0, 3.0]])
1067                         
1068        assert allclose(G.get_attributes(attribute_name='elevation'),
1069                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
1070       
1071        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
1072        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
1073       
1074        self.failUnless(G.get_geo_reference() == geo_reference1,
1075                         'test_writepts failed. Test geo_reference')
1076                         
1077        os.remove(fileName1)
1078        os.remove(fileName2)
1079       
1080    def test_ensure_absolute(self):
1081        points = [[2.0, 0.0],[1.0, 1.0],
1082                         [2.0, 0.0],[2.0, 2.0],
1083                         [1.0, 3.0],[2.0, 2.0]]
1084        new_points = ensure_absolute(points)
1085       
1086        assert allclose(new_points, points)
1087
1088        points = array([[2.0, 0.0],[1.0, 1.0],
1089                         [2.0, 0.0],[2.0, 2.0],
1090                         [1.0, 3.0],[2.0, 2.0]])
1091        new_points = ensure_absolute(points)
1092       
1093        assert allclose(new_points, points)
1094       
1095        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1096                         [2.0, 0.0],[2.0, 2.0],
1097                         [1.0, 3.0],[2.0, 2.0]])
1098       
1099        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1100
1101        data_points = zeros((ab_points.shape), Float)
1102        #Shift datapoints according to new origins
1103        for k in range(len(ab_points)):
1104            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1105            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1106        #print "data_points",data_points     
1107        new_points = ensure_absolute(data_points,
1108                                             geo_reference=mesh_origin)
1109        #print "new_points",new_points
1110        #print "ab_points",ab_points
1111           
1112        assert allclose(new_points, ab_points)
1113
1114        geo = Geo_reference(56,67,-56)
1115
1116        data_points = geo.change_points_geo_ref(ab_points)   
1117        new_points = ensure_absolute(data_points,
1118                                             geo_reference=geo)
1119        #print "new_points",new_points
1120        #print "ab_points",ab_points
1121           
1122        assert allclose(new_points, ab_points)
1123
1124
1125        geo_reference = Geo_reference(56, 100, 200)
1126        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1127        points = geo_reference.change_points_geo_ref(ab_points)
1128        attributes = [2, 4]
1129        #print "geo in points", points
1130        G = Geospatial_data(points, attributes,
1131                            geo_reference=geo_reference)
1132         
1133        new_points = ensure_absolute(G)
1134        #print "new_points",new_points
1135        #print "ab_points",ab_points
1136           
1137        assert allclose(new_points, ab_points)
1138
1139       
1140        fileName = tempfile.mktemp(".xya")
1141        file = open(fileName,"w")
1142        file.write("  elevation   speed \n\
11431.0 0.0 10.0 0.0\n\
11440.0 1.0 0.0 10.0\n\
11451.0 0.0 10.4 40.0\n\
1146#geocrap\n\
114756\n\
114810\n\
114920\n")
1150        file.close()
1151       
1152        ab_points = ensure_absolute(fileName)
1153        actual =  [[11, 20.0],[10.0, 21.0],[11.0, 20.0]]
1154        assert allclose(ab_points, actual)
1155        os.remove(fileName)
1156
1157       
1158    def test_ensure_geospatial(self):
1159        points = [[2.0, 0.0],[1.0, 1.0],
1160                         [2.0, 0.0],[2.0, 2.0],
1161                         [1.0, 3.0],[2.0, 2.0]]
1162        new_points = ensure_geospatial(points)
1163       
1164        assert allclose(new_points.get_data_points(absolute = True), points)
1165
1166        points = array([[2.0, 0.0],[1.0, 1.0],
1167                         [2.0, 0.0],[2.0, 2.0],
1168                         [1.0, 3.0],[2.0, 2.0]])
1169        new_points = ensure_geospatial(points)
1170       
1171        assert allclose(new_points.get_data_points(absolute = True), points)
1172       
1173        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1174                         [2.0, 0.0],[2.0, 2.0],
1175                         [1.0, 3.0],[2.0, 2.0]])
1176       
1177        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1178
1179        data_points = zeros((ab_points.shape), Float)
1180        #Shift datapoints according to new origins
1181        for k in range(len(ab_points)):
1182            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1183            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1184        #print "data_points",data_points     
1185        new_geospatial = ensure_geospatial(data_points,
1186                                             geo_reference=mesh_origin)
1187        new_points = new_geospatial.get_data_points(absolute=True)
1188        #print "new_points",new_points
1189        #print "ab_points",ab_points
1190           
1191        assert allclose(new_points, ab_points)
1192
1193        geo = Geo_reference(56,67,-56)
1194
1195        data_points = geo.change_points_geo_ref(ab_points)   
1196        new_geospatial = ensure_geospatial(data_points,
1197                                             geo_reference=geo)
1198        new_points = new_geospatial.get_data_points(absolute=True)
1199        #print "new_points",new_points
1200        #print "ab_points",ab_points
1201           
1202        assert allclose(new_points, ab_points)
1203
1204
1205        geo_reference = Geo_reference(56, 100, 200)
1206        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1207        points = geo_reference.change_points_geo_ref(ab_points)
1208        attributes = [2, 4]
1209        #print "geo in points", points
1210        G = Geospatial_data(points, attributes,
1211                            geo_reference=geo_reference)
1212         
1213        new_geospatial  = ensure_geospatial(G)
1214        new_points = new_geospatial.get_data_points(absolute=True)
1215        #print "new_points",new_points
1216        #print "ab_points",ab_points
1217           
1218        assert allclose(new_points, ab_points)
1219       
1220    def test_isinstance(self):
1221
1222        import os
1223       
1224        fileName = tempfile.mktemp(".xya")
1225        file = open(fileName,"w")
1226        file.write("  elevation   speed \n\
12271.0 0.0 10.0 0.0\n\
12280.0 1.0 0.0 10.0\n\
12291.0 0.0 10.4 40.0\n\
1230#geocrap\n\
123156\n\
123256.6\n\
12333\n")
1234        file.close()
1235
1236        results = Geospatial_data(fileName)
1237        assert allclose(results.get_data_points(absolute=False), \
1238                        [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1239        assert allclose(results.get_attributes(attribute_name='elevation'), \
1240                        [10.0, 0.0, 10.4])
1241        assert allclose(results.get_attributes(attribute_name='speed'), \
1242                        [0.0, 10.0, 40.0])
1243
1244        os.remove(fileName)
1245       
1246    def test_delimiter(self):
1247       
1248        try:
1249            G = Geospatial_data(delimiter=',')
1250#            results = Geospatial_data(file_name = fileName)
1251#            dict = import_points_file(fileName)
1252        except ValueError:
1253            pass
1254        else:
1255            msg = 'Instance with No fileName but has a delimiter\
1256                  did not raise error!'
1257            raise msg
1258
1259    def test_no_constructors(self):
1260       
1261        try:
1262            G = Geospatial_data()
1263#            results = Geospatial_data(file_name = fileName)
1264#            dict = import_points_file(fileName)
1265        except ValueError:
1266            pass
1267        else:
1268            msg = 'Instance must have a filename or data points'
1269            raise msg       
1270
1271    def test_check_geo_reference(self):
1272        """
1273        checks geo reference details are OK. eg can be called '#geo reference'
1274        if not throws a clear error message
1275        """
1276        import os
1277        fileName = tempfile.mktemp(".xya")
1278        file = open(fileName,"w")
1279        file.write("  elevation  \n\
12801.0 0.0 10.0\n\
12810.0 1.0 0.0\n\
12821.0 0.0 10.4\n\
1283#ge oreference\n\
128456\n\
12851.1\n\
12861.0\n")
1287
1288        file.close()
1289        results = Geospatial_data(fileName)
1290        assert allclose(results.get_geo_reference().get_xllcorner(), 1.1)
1291        assert allclose(results.get_geo_reference().get_yllcorner(), 1.0)
1292
1293        os.remove(fileName)
1294       
1295        fileName = tempfile.mktemp(".xya")
1296        file = open(fileName,"w")
1297        file.write("  elevation  \n\
12981.0 0.0 10.0\n\
12990.0 1.0 0.0\n\
13001.0 0.0 10.4\n")
1301
1302        file.close()
1303        results = Geospatial_data(fileName)
1304       
1305        os.remove(fileName)
1306       
1307    def test_check_geo_reference1(self):
1308        """
1309        checks geo reference details are OK. eg can be called '#geo reference'
1310        if not throws a clear error message
1311        """
1312        import os
1313        fileName = tempfile.mktemp(".xya")
1314        file = open(fileName,"w")
1315        file.write("  elevation  \n\
13161.0 0.0 10.0\n\
13170.0 1.0 0.0\n\
13181.0 0.0 10.4\n\
1319#geo t t\n\
132056\n\
13211.1\n"
1322)
1323        file.close()
1324
1325        try:
1326            results = Geospatial_data(fileName, delimiter = " ")
1327        except IOError:
1328            pass
1329        else:
1330            msg = 'Geo reference data format is incorrect'
1331            raise msg       
1332
1333
1334        os.remove(fileName)
1335
1336
1337       
1338    def test_lat_long(self):
1339        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1340        lon_gong = degminsec2decimal_degrees(150,55,0.)
1341       
1342        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1343        lon_2 = degminsec2decimal_degrees(150,00,0.)
1344       
1345        lats = [lat_gong, lat_2]
1346        longs = [lon_gong, lon_2]
1347        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1348
1349        points = gsd.get_data_points(absolute=True)
1350       
1351        assert allclose(points[0][0], 308728.009)
1352        assert allclose(points[0][1], 6180432.601)
1353        assert allclose(points[1][0],  222908.705)
1354        assert allclose(points[1][1], 6233785.284)
1355        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1356                        'Bad zone error!')
1357       
1358        try:
1359            results = Geospatial_data(latitudes=lats)
1360        except ValueError:
1361            pass
1362        else:
1363            self.failUnless(0 ==1,  'Error not thrown error!')
1364        try:
1365            results = Geospatial_data(latitudes=lats)
1366        except ValueError:
1367            pass
1368        else:
1369            self.failUnless(0 ==1,  'Error not thrown error!')
1370        try:
1371            results = Geospatial_data(longitudes=lats)
1372        except ValueError:
1373            pass
1374        else:
1375            self.failUnless(0 ==1, 'Error not thrown error!')
1376        try:
1377            results = Geospatial_data(latitudes=lats, longitudes=longs,
1378                                      geo_reference="p")
1379        except ValueError:
1380            pass
1381        else:
1382            self.failUnless(0 ==1,  'Error not thrown error!')
1383           
1384        try:
1385            results = Geospatial_data(latitudes=lats, longitudes=longs,
1386                                      data_points=12)
1387        except ValueError:
1388            pass
1389        else:
1390            self.failUnless(0 ==1,  'Error not thrown error!')
1391
1392    def test_lat_long2(self):
1393        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1394        lon_gong = degminsec2decimal_degrees(150,55,0.)
1395       
1396        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1397        lon_2 = degminsec2decimal_degrees(150,00,0.)
1398       
1399        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1400        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1401
1402        points = gsd.get_data_points(absolute=True)
1403       
1404        assert allclose(points[0][0], 308728.009)
1405        assert allclose(points[0][1], 6180432.601)
1406        assert allclose(points[1][0],  222908.705)
1407        assert allclose(points[1][1], 6233785.284)
1408        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1409                        'Bad zone error!')
1410
1411        try:
1412            results = Geospatial_data(points_are_lats_longs=True)
1413        except ValueError:
1414            pass
1415        else:
1416            self.failUnless(0 ==1,  'Error not thrown error!')
1417
1418    def test_len(self):
1419       
1420        points = [[1.0, 2.1], [3.0, 5.3]]
1421        G = Geospatial_data(points)
1422        self.failUnless(2 ==len(G),  'Len error!')
1423       
1424        points = [[1.0, 2.1]]
1425        G = Geospatial_data(points)
1426        self.failUnless(1 ==len(G),  'Len error!')
1427
1428        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1429        G = Geospatial_data(points)
1430        self.failUnless(4 ==len(G),  'Len error!')
1431         
1432if __name__ == "__main__":
1433
1434    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_ensure_geospatial')
1435    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1436    runner = unittest.TextTestRunner()
1437    runner.run(suite)
1438
1439   
Note: See TracBrowser for help on using the repository browser.