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

Last change on this file since 5729 was 5729, checked in by ole, 16 years ago

Refactored get_flow_through_cross_section and added a runtime version to
shallow_water_domain. It still needs some work and testing.
File size: 71.5 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import os
6from Numeric import zeros, array, allclose, concatenate,sort
7from math import sqrt, pi
8import tempfile
9from sets import ImmutableSet
10
11from anuga.geospatial_data.geospatial_data import *
12from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
13from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
14from anuga.utilities.anuga_exceptions import ANUGAError
15from anuga.utilities.system_tools import get_host_name
16
17class Test_Geospatial_data(unittest.TestCase):
18    def setUp(self):
19        pass
20
21    def tearDown(self):
22        pass
23
24
25    def test_0(self):
26        #Basic points
27        from anuga.coordinate_transforms.geo_reference import Geo_reference
28       
29        points = [[1.0, 2.1], [3.0, 5.3]]
30        G = Geospatial_data(points)
31
32        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
33
34        # Check __repr__
35        # FIXME (Ole): Is this really machine independent?
36        rep = `G`
37        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
38
39        msg = 'Representation %s is not equal to %s' %(rep, ref)
40        assert rep == ref, msg
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        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
58        assert allclose(G.attributes.values()[0], [2, 4])
59       
60
61    def test_2(self):
62        from anuga.coordinate_transforms.geo_reference import Geo_reference
63        points = [[1.0, 2.1], [3.0, 5.3]]
64        attributes = [2, 4]
65        G = Geospatial_data(points, attributes,
66                            geo_reference=Geo_reference(56, 100, 200))
67
68        assert G.geo_reference.zone == 56
69        assert G.geo_reference.xllcorner == 100
70        assert G.geo_reference.yllcorner == 200
71
72
73    def test_get_attributes_1(self):
74        from anuga.coordinate_transforms.geo_reference import Geo_reference
75        points = [[1.0, 2.1], [3.0, 5.3]]
76        attributes = [2, 4]
77        G = Geospatial_data(points, attributes,
78                            geo_reference=Geo_reference(56, 100, 200))
79
80
81        P = G.get_data_points(absolute=False)
82        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
83
84        P = G.get_data_points(absolute=True)
85        assert allclose(P, [[101.0, 202.1], [103.0, 205.3]])       
86
87        V = G.get_attributes() #Simply get them
88        assert allclose(V, [2, 4])
89
90        V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name
91        assert allclose(V, [2, 4])
92
93    def test_get_attributes_2(self):
94        #Multiple attributes
95       
96       
97        from anuga.coordinate_transforms.geo_reference import Geo_reference
98        points = [[1.0, 2.1], [3.0, 5.3]]
99        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
100        G = Geospatial_data(points, attributes,
101                            geo_reference=Geo_reference(56, 100, 200),
102                            default_attribute_name='a1')
103
104
105        P = G.get_data_points(absolute=False)
106        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
107       
108        V = G.get_attributes() #Get default attribute
109        assert allclose(V, [2, 4])
110
111        V = G.get_attributes('a0') #Get by name
112        assert allclose(V, [0, 0])
113
114        V = G.get_attributes('a1') #Get by name
115        assert allclose(V, [2, 4])
116
117        V = G.get_attributes('a2') #Get by name
118        assert allclose(V, [79.4, -7])
119
120        try:
121            V = G.get_attributes('hdnoatedu') #Invalid
122        except AssertionError:
123            pass
124        else:
125            raise 'Should have raised exception' 
126
127    def test_get_data_points(self):
128        points_ab = [[12.5,34.7],[-4.5,-60.0]]
129        x_p = -10
130        y_p = -40
131        geo_ref = Geo_reference(56, x_p, y_p)
132        points_rel = geo_ref.change_points_geo_ref(points_ab)
133       
134        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
135
136        results = spatial.get_data_points(absolute=False)
137       
138        assert allclose(results, points_rel)
139       
140        x_p = -1770
141        y_p = 4.01
142        geo_ref = Geo_reference(56, x_p, y_p)
143        points_rel = geo_ref.change_points_geo_ref(points_ab)
144        results = spatial.get_data_points \
145                  ( geo_reference=geo_ref)
146       
147        assert allclose(results, points_rel)
148
149 
150    def test_get_data_points_lat_long(self):
151        # lat long [-30.],[130]
152        #Zone:   52   
153        #Easting:  596450.153  Northing: 6680793.777
154        # lat long [-32.],[131]
155        #Zone:   52   
156        #Easting:  688927.638  Northing: 6457816.509
157       
158        points_Lat_long = [[-30.,130], [-32,131]]
159       
160        spatial = Geospatial_data(latitudes=[-30, -32.],
161                                  longitudes=[130, 131])
162
163        results = spatial.get_data_points(as_lat_long=True)
164        #print "test_get_data_points_lat_long - results", results
165        #print "points_Lat_long",points_Lat_long
166        assert allclose(results, points_Lat_long)
167     
168    def test_get_data_points_lat_longII(self):
169        # x,y  North,east long,lat
170        boundary_polygon = [[ 250000, 7630000]]
171        zone = 50
172       
173        geo_reference = Geo_reference(zone=zone)
174        geo = Geospatial_data(boundary_polygon,geo_reference=geo_reference)
175        seg_lat_long = geo.get_data_points(as_lat_long=True)
176        lat_result = degminsec2decimal_degrees(-21,24,54)
177        long_result = degminsec2decimal_degrees(114,35,17.89)
178        #print "seg_lat_long", seg_lat_long [0][0]
179        #print "lat_result",lat_result
180        assert allclose(seg_lat_long[0][0], lat_result)#lat
181        assert allclose(seg_lat_long[0][1], long_result)#long
182
183
184    def test_get_data_points_lat_longIII(self):
185        # x,y  North,east long,lat
186        #for northern hemisphere
187        boundary_polygon = [[419944.8, 918642.4]]
188        zone = 47
189       
190        geo_reference = Geo_reference(zone=zone)
191        geo = Geospatial_data(boundary_polygon,
192                              geo_reference=geo_reference)
193                             
194        seg_lat_long = geo.get_data_points(as_lat_long=True,
195                                           isSouthHemisphere=False)
196                                           
197        lat_result = degminsec2decimal_degrees(8.31,0,0)
198        long_result = degminsec2decimal_degrees(98.273,0,0)
199        #print "seg_lat_long", seg_lat_long [0]
200        #print "lat_result",lat_result
201        assert allclose(seg_lat_long[0][0], lat_result)#lat
202        assert allclose(seg_lat_long[0][1], long_result)#long
203
204
205             
206    def test_set_geo_reference(self):
207        points_ab = [[12.5,34.7],[-4.5,-60.0]]
208        x_p = -10
209        y_p = -40
210        geo_ref = Geo_reference(56, x_p, y_p)
211        points_rel = geo_ref.change_points_geo_ref(points_ab)
212        spatial = Geospatial_data(points_rel)
213        # since the geo_ref wasn't set
214        assert not allclose( points_ab, spatial.get_data_points(absolute=True))
215       
216        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
217        assert allclose( points_ab, spatial.get_data_points(absolute=True))
218       
219        x_p = 10
220        y_p = 400
221        new_geo_ref = Geo_reference(56, x_p, y_p)
222        spatial.set_geo_reference(new_geo_ref)
223        assert allclose( points_ab, spatial.get_data_points(absolute=True))
224       
225       
226       
227    def test_conversions_to_points_dict(self):
228        #test conversions to points_dict
229       
230       
231        from anuga.coordinate_transforms.geo_reference import Geo_reference
232        points = [[1.0, 2.1], [3.0, 5.3]]
233        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
234        G = Geospatial_data(points, attributes,
235                            geo_reference=Geo_reference(56, 100, 200),
236                            default_attribute_name='a1')
237
238
239        points_dict = geospatial_data2points_dictionary(G)
240
241        assert points_dict.has_key('pointlist')
242        assert points_dict.has_key('attributelist')       
243        assert points_dict.has_key('geo_reference')
244
245        assert allclose( points_dict['pointlist'], points )
246
247        A = points_dict['attributelist']
248        assert A.has_key('a0')
249        assert A.has_key('a1')
250        assert A.has_key('a2')       
251
252        assert allclose( A['a0'], [0, 0] )
253        assert allclose( A['a1'], [2, 4] )       
254        assert allclose( A['a2'], [79.4, -7] )
255
256
257        geo = points_dict['geo_reference']
258        assert geo is G.geo_reference
259
260
261    def test_conversions_from_points_dict(self):
262        """test conversions from points_dict
263        """
264
265        from anuga.coordinate_transforms.geo_reference import Geo_reference
266       
267        points = [[1.0, 2.1], [3.0, 5.3]]
268        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
269
270        points_dict = {}
271        points_dict['pointlist'] = points
272        points_dict['attributelist'] = attributes
273        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
274       
275
276        G = points_dictionary2geospatial_data(points_dict)
277
278        P = G.get_data_points(absolute=False)
279        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
280       
281        #V = G.get_attribute_values() #Get default attribute
282        #assert allclose(V, [2, 4])
283
284        V = G.get_attributes('a0') #Get by name
285        assert allclose(V, [0, 0])
286
287        V = G.get_attributes('a1') #Get by name
288        assert allclose(V, [2, 4])
289
290        V = G.get_attributes('a2') #Get by name
291        assert allclose(V, [79.4, -7])
292
293    def test_add(self):
294        """ test the addition of two geospatical objects
295            no geo_reference see next test
296        """
297        points = [[1.0, 2.1], [3.0, 5.3]]
298        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
299        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
300        G1 = Geospatial_data(points, attributes)       
301        G2 = Geospatial_data(points, attributes1) 
302       
303#        g3 = geospatial_data2points_dictionary(G1)
304#        print 'g3=', g3
305       
306        G = G1 + G2
307
308        assert G.attributes.has_key('depth')
309        assert G.attributes.has_key('elevation')
310        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
311        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
312        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
313                                              [1.0, 2.1], [3.0, 5.3]])
314       
315    def test_addII(self):
316        """ test the addition of two geospatical objects
317            no geo_reference see next test
318        """
319        points = [[1.0, 2.1], [3.0, 5.3]]
320        attributes = {'depth':[2, 4]}
321        G1 = Geospatial_data(points, attributes) 
322       
323        points = [[5.0, 2.1], [3.0, 50.3]]
324        attributes = {'depth':[200, 400]}
325        G2 = Geospatial_data(points, attributes)
326       
327#        g3 = geospatial_data2points_dictionary(G1)
328#        print 'g3=', g3
329       
330        G = G1 + G2
331
332        assert G.attributes.has_key('depth') 
333        assert G.attributes.keys(), ['depth']
334        assert allclose(G.attributes['depth'], [2, 4, 200, 400])
335        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
336                                              [5.0, 2.1], [3.0, 50.3]])
337    def test_add_with_geo (self):
338        """
339        Difference in Geo_reference resolved
340        """
341        points1 = [[1.0, 2.1], [3.0, 5.3]]
342        points2 = [[5.0, 6.1], [6.0, 3.3]]
343        attributes1 = [2, 4]
344        attributes2 = [5, 76]
345        geo_ref1= Geo_reference(55, 1.0, 2.0)
346        geo_ref2 = Geo_reference(zone=55,
347                                 xllcorner=0.1,
348                                 yllcorner=3.0,
349                                 datum='wgs84',
350                                 projection='UTM',
351                                 units='m')
352                               
353        G1 = Geospatial_data(points1, attributes1, geo_ref1)
354        G2 = Geospatial_data(points2, attributes2, geo_ref2)
355
356        #Check that absolute values are as expected
357        P1 = G1.get_data_points(absolute=True)
358        assert allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
359
360        P2 = G2.get_data_points(absolute=True)
361        assert allclose(P2, [[5.1, 9.1], [6.1, 6.3]])       
362       
363        G = G1 + G2
364
365        # Check absoluteness
366        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
367        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
368
369        P = G.get_data_points(absolute=True)
370
371        #P_relative = G.get_data_points(absolute=False)
372        #
373        #assert allclose(P_relative, P - [0.1, 2.0])
374
375        assert allclose(P, concatenate( (P1,P2) ))
376        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
377                            [5.1, 9.1], [6.1, 6.3]])
378       
379
380
381       
382
383    def test_add_with_geo_absolute (self):
384        """
385        Difference in Geo_reference resolved
386        """
387        points1 = array([[2.0, 4.1], [4.0, 7.3]])
388        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
389        attributes1 = [2, 4]
390        attributes2 = [5, 76]
391        geo_ref1= Geo_reference(55, 1.0, 2.0)
392        geo_ref2 = Geo_reference(55, 2.0, 3.0)
393
394       
395                               
396        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
397                             attributes1, geo_ref1)
398       
399        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
400                             attributes2, geo_ref2)
401
402        #Check that absolute values are as expected
403        P1 = G1.get_data_points(absolute=True)
404        assert allclose(P1, points1)
405
406        P1 = G1.get_data_points(absolute=False)
407        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
408
409        P2 = G2.get_data_points(absolute=True)
410        assert allclose(P2, points2)
411
412        P2 = G2.get_data_points(absolute=False)
413        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
414       
415        G = G1 + G2
416       
417        #assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
418        #assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
419
420        P = G.get_data_points(absolute=True)
421
422        #P_relative = G.get_data_points(absolute=False)
423        #
424        #assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
425
426        assert allclose(P, concatenate( (points1,points2) ))
427
428
429    def test_add_with_None(self):
430        """ test that None can be added to a geospatical objects
431        """
432       
433        points1 = array([[2.0, 4.1], [4.0, 7.3]])
434        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
435
436        geo_ref1= Geo_reference(55, 1.0, 2.0)
437        geo_ref2 = Geo_reference(zone=55,
438                                 xllcorner=0.1,
439                                 yllcorner=3.0,
440                                 datum='wgs84',
441                                 projection='UTM',
442                                 units='m')
443       
444
445        attributes1 = {'depth':[2, 4.7], 'elevation':[6.1, 5]}
446        attributes2 = {'depth':[-2.3, 4], 'elevation':[2.5, 1]}
447
448
449        G1 = Geospatial_data(points1, attributes1, geo_ref1)
450        assert allclose(G1.get_geo_reference().get_xllcorner(), 1.0)
451        assert allclose(G1.get_geo_reference().get_yllcorner(), 2.0)
452        assert G1.attributes.has_key('depth')
453        assert G1.attributes.has_key('elevation')
454        assert allclose(G1.attributes['depth'], [2, 4.7])
455        assert allclose(G1.attributes['elevation'], [6.1, 5])       
456       
457        G2 = Geospatial_data(points2, attributes2, geo_ref2)
458        assert allclose(G2.get_geo_reference().get_xllcorner(), 0.1)
459        assert allclose(G2.get_geo_reference().get_yllcorner(), 3.0)
460        assert G2.attributes.has_key('depth')
461        assert G2.attributes.has_key('elevation')
462        assert allclose(G2.attributes['depth'], [-2.3, 4])
463        assert allclose(G2.attributes['elevation'], [2.5, 1])       
464
465        #Check that absolute values are as expected
466        P1 = G1.get_data_points(absolute=True)
467        assert allclose(P1, [[3.0, 6.1], [5.0, 9.3]])
468
469        P2 = G2.get_data_points(absolute=True)
470        assert allclose(P2, [[5.2, 12.1], [6.2, 9.3]])       
471
472        # Normal add
473        G = G1 + None
474
475        assert G.attributes.has_key('depth')
476        assert G.attributes.has_key('elevation')
477        assert allclose(G.attributes['depth'], [2, 4.7])
478        assert allclose(G.attributes['elevation'], [6.1, 5])       
479
480        # Points are now absolute.
481        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
482        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
483        P = G.get_data_points(absolute=True)       
484        assert allclose(P, [[3.0, 6.1], [5.0, 9.3]])
485
486
487        G = G2 + None
488        assert G.attributes.has_key('depth')
489        assert G.attributes.has_key('elevation')
490        assert allclose(G.attributes['depth'], [-2.3, 4])
491        assert allclose(G.attributes['elevation'], [2.5, 1])       
492
493        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
494        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
495        P = G.get_data_points(absolute=True)       
496        assert allclose(P, [[5.2, 12.1], [6.2, 9.3]])
497       
498
499
500        # Reverse add
501        G = None + G1
502
503        assert G.attributes.has_key('depth')
504        assert G.attributes.has_key('elevation')
505        assert allclose(G.attributes['depth'], [2, 4.7])
506        assert allclose(G.attributes['elevation'], [6.1, 5])       
507
508        # Points are now absolute.
509        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
510        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
511        P = G.get_data_points(absolute=True)       
512        assert allclose(P, [[3.0, 6.1], [5.0, 9.3]])       
513
514
515        G = None + G2
516        assert G.attributes.has_key('depth')
517        assert G.attributes.has_key('elevation')
518        assert allclose(G.attributes['depth'], [-2.3, 4])
519        assert allclose(G.attributes['elevation'], [2.5, 1])       
520
521        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
522        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
523        P = G.get_data_points(absolute=True)       
524        assert allclose(P, [[5.2, 12.1], [6.2, 9.3]])
525
526       
527
528       
529                           
530       
531    def test_clip0(self):
532        """test_clip0(self):
533       
534        Test that point sets can be clipped by a polygon
535        """
536       
537        from anuga.coordinate_transforms.geo_reference import Geo_reference
538       
539        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
540                  [0, 0], [2.4, 3.3]]
541        G = Geospatial_data(points)
542
543        # First try the unit square   
544        U = [[0,0], [1,0], [1,1], [0,1]] 
545        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
546
547        # Then a more complex polygon
548        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
549        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
550        G = Geospatial_data(points)
551
552        assert allclose(G.clip(polygon).get_data_points(),
553                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
554
555    def test_clip0_with_attributes(self):
556        """test_clip0_with_attributes(self):
557       
558        Test that point sets with attributes can be clipped by a polygon
559        """
560       
561        from anuga.coordinate_transforms.geo_reference import Geo_reference
562       
563        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
564                  [0, 0], [2.4, 3.3]]
565
566        attributes = [2, -4, 5, 76, -2, 0.1, 3]
567        att_dict = {'att1': attributes,
568                    'att2': array(attributes)+1}
569       
570        G = Geospatial_data(points, att_dict)
571
572        # First try the unit square   
573        U = [[0,0], [1,0], [1,1], [0,1]] 
574        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
575        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
576        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])               
577
578        # Then a more complex polygon
579        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
580        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
581
582        # This time just one attribute
583        attributes = [2, -4, 5, 76, -2, 0.1]
584        G = Geospatial_data(points, attributes)
585
586        assert allclose(G.clip(polygon).get_data_points(),
587                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
588        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
589       
590
591    def test_clip1(self):
592        """test_clip1(self):
593       
594        Test that point sets can be clipped by a polygon given as
595        another Geospatial dataset
596        """
597       
598        from anuga.coordinate_transforms.geo_reference import Geo_reference
599       
600        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
601                  [0, 0], [2.4, 3.3]]
602        attributes = [2, -4, 5, 76, -2, 0.1, 3]
603        att_dict = {'att1': attributes,
604                    'att2': array(attributes)+1}
605        G = Geospatial_data(points, att_dict)
606       
607        # First try the unit square   
608        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
609        assert allclose(G.clip(U).get_data_points(),
610                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
611
612        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
613        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                       
614       
615        # Then a more complex polygon
616        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
617        attributes = [2, -4, 5, 76, -2, 0.1]       
618        G = Geospatial_data(points, attributes)
619        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
620       
621
622        assert allclose(G.clip(polygon).get_data_points(),
623                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
624        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
625       
626
627    def test_clip0_outside(self):
628        """test_clip0_outside(self):
629       
630        Test that point sets can be clipped outside of a polygon
631        """
632       
633        from anuga.coordinate_transforms.geo_reference import Geo_reference
634       
635        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
636                  [0, 0], [2.4, 3.3]]
637        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
638        G = Geospatial_data(points, attributes)
639
640        # First try the unit square   
641        U = [[0,0], [1,0], [1,1], [0,1]]
642        assert allclose(G.clip_outside(U).get_data_points(),
643                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
644        #print G.clip_outside(U).get_attributes()
645        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])       
646       
647
648        # Then a more complex polygon
649        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
650        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
651        attributes = [2, -4, 5, 76, -2, 0.1]       
652        G = Geospatial_data(points, attributes)
653
654        assert allclose(G.clip_outside(polygon).get_data_points(),
655                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
656        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])               
657
658
659    def test_clip1_outside(self):
660        """test_clip1_outside(self):
661       
662        Test that point sets can be clipped outside of a polygon given as
663        another Geospatial dataset
664        """
665       
666        from anuga.coordinate_transforms.geo_reference import Geo_reference
667       
668        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
669                  [0, 0], [2.4, 3.3]]
670        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
671        G = Geospatial_data(points, attributes)       
672
673        # First try the unit square   
674        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
675        assert allclose(G.clip_outside(U).get_data_points(),
676                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
677        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])       
678
679        # Then a more complex polygon
680        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
681        attributes = [2, -4, 5, 76, -2, 0.1]       
682        G = Geospatial_data(points, attributes)
683
684        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
685       
686
687        assert allclose(G.clip_outside(polygon).get_data_points(),
688                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
689        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
690       
691
692
693    def test_clip1_inside_outside(self):
694        """test_clip1_inside_outside(self):
695       
696        Test that point sets can be clipped outside of a polygon given as
697        another Geospatial dataset
698        """
699       
700        from anuga.coordinate_transforms.geo_reference import Geo_reference
701       
702        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
703                  [0, 0], [2.4, 3.3]]
704        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
705        G = Geospatial_data(points, attributes)
706
707        # First try the unit square   
708        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
709        G1 = G.clip(U)
710        assert allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
711        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
712       
713        G2 = G.clip_outside(U)
714        assert allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
715                                              [3.0, 5.3], [2.4, 3.3]])
716        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])               
717
718       
719        # New ordering
720        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\
721                     [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
722
723        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
724       
725        assert allclose((G1+G2).get_data_points(), new_points)
726        assert allclose((G1+G2).get_attributes(), new_attributes)
727
728        G = G1+G2
729        FN = 'test_combine.pts'
730        G.export_points_file(FN)
731
732
733        # Read it back in
734        G3 = Geospatial_data(FN)
735
736
737        # Check result
738        assert allclose(G3.get_data_points(), new_points)       
739        assert allclose(G3.get_attributes(), new_attributes)       
740       
741        os.remove(FN)
742
743       
744    def test_load_csv(self):
745       
746        import os
747        import tempfile
748       
749        fileName = tempfile.mktemp(".csv")
750        file = open(fileName,"w")
751        file.write("x,y,elevation speed \n\
7521.0 0.0 10.0 0.0\n\
7530.0 1.0 0.0 10.0\n\
7541.0 0.0 10.4 40.0\n")
755        file.close()
756        #print fileName
757        results = Geospatial_data(fileName)
758        os.remove(fileName)
759#        print 'data', results.get_data_points()
760        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
761                                                    [1.0, 0.0]])
762        assert allclose(results.get_attributes(attribute_name='elevation'),
763                        [10.0, 0.0, 10.4])
764        assert allclose(results.get_attributes(attribute_name='speed'),
765                        [0.0, 10.0, 40.0])
766
767
768  ###################### .CSV ##############################
769
770    def test_load_csv_lat_long_bad_blocking(self):
771        """
772        test_load_csv_lat_long_bad_blocking(self):
773        Different zones in Geo references
774        """
775        fileName = tempfile.mktemp(".csv")
776        file = open(fileName,"w")
777        file.write("Lati,LONG,z \n\
778-25.0,180.0,452.688000\n\
779-34,150.0,459.126000\n")
780        file.close()
781       
782        results = Geospatial_data(fileName, max_read_lines=1,
783                                  load_file_now=False)
784       
785        #for i in results:
786        #    pass
787        try:
788            for i in results:
789                pass
790        except ANUGAError:
791            pass
792        else:
793            msg = 'Different zones in Geo references not caught.'
794            raise msg       
795       
796        os.remove(fileName)
797       
798    def test_load_csv(self):
799       
800        fileName = tempfile.mktemp(".txt")
801        file = open(fileName,"w")
802        file.write(" x,y, elevation ,  speed \n\
8031.0, 0.0, 10.0, 0.0\n\
8040.0, 1.0, 0.0, 10.0\n\
8051.0, 0.0 ,10.4, 40.0\n")
806        file.close()
807
808        results = Geospatial_data(fileName, max_read_lines=2)
809
810
811        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
812        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
813        assert allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
814
815        # Blocking
816        geo_list = []
817        for i in results:
818            geo_list.append(i)
819           
820        assert allclose(geo_list[0].get_data_points(),
821                        [[1.0, 0.0],[0.0, 1.0]])
822
823        assert allclose(geo_list[0].get_attributes(attribute_name='elevation'),
824                        [10.0, 0.0])
825        assert allclose(geo_list[1].get_data_points(),
826                        [[1.0, 0.0]])       
827        assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
828                        [10.4])
829           
830        os.remove(fileName)         
831       
832    def test_load_csv_bad(self):
833        """ test_load_csv_bad(self):
834        header column, body column missmatch
835        (Format error)
836        """
837        import os
838       
839        fileName = tempfile.mktemp(".txt")
840        file = open(fileName,"w")
841        file.write(" elevation ,  speed \n\
8421.0, 0.0, 10.0, 0.0\n\
8430.0, 1.0, 0.0, 10.0\n\
8441.0, 0.0 ,10.4, 40.0\n")
845        file.close()
846
847        results = Geospatial_data(fileName, max_read_lines=2,
848                                  load_file_now=False)
849
850        # Blocking
851        geo_list = []
852        #for i in results:
853        #    geo_list.append(i)
854        try:
855            for i in results:
856                geo_list.append(i)
857        except SyntaxError:
858            pass
859        else:
860            msg = 'bad file did not raise error!'
861            raise msg
862        os.remove(fileName)
863
864    def test_load_csv_badII(self):
865        """ test_load_csv_bad(self):
866        header column, body column missmatch
867        (Format error)
868        """
869        import os
870       
871        fileName = tempfile.mktemp(".txt")
872        file = open(fileName,"w")
873        file.write(" x,y,elevation ,  speed \n\
8741.0, 0.0, 10.0, 0.0\n\
8750.0, 1.0, 0.0, 10\n\
8761.0, 0.0 ,10.4 yeah\n")
877        file.close()
878
879        results = Geospatial_data(fileName, max_read_lines=2,
880                                  load_file_now=False)
881
882        # Blocking
883        geo_list = []
884        #for i in results:
885        #    geo_list.append(i)
886        try:
887            for i in results:
888                geo_list.append(i)
889        except SyntaxError:
890            pass
891        else:
892            msg = 'bad file did not raise error!'
893            raise msg
894        os.remove(fileName)
895
896    def test_load_csv_badIII(self):
897        """ test_load_csv_bad(self):
898        space delimited
899        """
900        import os
901       
902        fileName = tempfile.mktemp(".txt")
903        file = open(fileName,"w")
904        file.write(" x y elevation   speed \n\
9051. 0.0 10.0 0.0\n\
9060.0 1.0 0.0 10.0\n\
9071.0 0.0 10.4 40.0\n")
908        file.close()
909
910        try:
911            results = Geospatial_data(fileName, max_read_lines=2,
912                                      load_file_now=True)
913        except SyntaxError:
914            pass
915        else:
916            msg = 'bad file did not raise error!'
917            raise msg
918        os.remove(fileName)
919       
920    def test_load_csv_badIV(self):
921        """ test_load_csv_bad(self):
922        header column, body column missmatch
923        (Format error)
924        """
925        import os
926       
927        fileName = tempfile.mktemp(".txt")
928        file = open(fileName,"w")
929        file.write(" x,y,elevation ,  speed \n\
9301.0, 0.0, 10.0, wow\n\
9310.0, 1.0, 0.0, ha\n\
9321.0, 0.0 ,10.4, yeah\n")
933        file.close()
934
935        results = Geospatial_data(fileName, max_read_lines=2,
936                                  load_file_now=False)
937
938        # Blocking
939        geo_list = []
940        #for i in results:
941         #   geo_list.append(i)
942        try:
943            for i in results:
944                geo_list.append(i)
945        except SyntaxError:
946            pass
947        else:
948            msg = 'bad file did not raise error!'
949            raise msg
950        os.remove(fileName)
951
952    def test_load_pts_blocking(self):
953        #This is pts!
954       
955        import os
956       
957        fileName = tempfile.mktemp(".txt")
958        file = open(fileName,"w")
959        file.write(" x,y, elevation ,  speed \n\
9601.0, 0.0, 10.0, 0.0\n\
9610.0, 1.0, 0.0, 10.0\n\
9621.0, 0.0 ,10.4, 40.0\n")
963        file.close()
964
965        pts_file = tempfile.mktemp(".pts")
966       
967        convert = Geospatial_data(fileName)
968        convert.export_points_file(pts_file)
969        results = Geospatial_data(pts_file, max_read_lines=2)
970
971        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
972                                                    [1.0, 0.0]])
973        assert allclose(results.get_attributes(attribute_name='elevation'),
974                        [10.0, 0.0, 10.4])
975        assert allclose(results.get_attributes(attribute_name='speed'),
976                        [0.0, 10.0, 40.0])
977
978        # Blocking
979        geo_list = []
980        for i in results:
981            geo_list.append(i) 
982        assert allclose(geo_list[0].get_data_points(),
983                        [[1.0, 0.0],[0.0, 1.0]])
984        assert allclose(geo_list[0].get_attributes(attribute_name='elevation'),
985                        [10.0, 0.0])
986        assert allclose(geo_list[1].get_data_points(),
987                        [[1.0, 0.0]])       
988        assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
989                        [10.4])
990           
991        os.remove(fileName) 
992        os.remove(pts_file)               
993
994    def verbose_test_load_pts_blocking(self):
995       
996        import os
997       
998        fileName = tempfile.mktemp(".txt")
999        file = open(fileName,"w")
1000        file.write(" x,y, elevation ,  speed \n\
10011.0, 0.0, 10.0, 0.0\n\
10020.0, 1.0, 0.0, 10.0\n\
10031.0, 0.0, 10.0, 0.0\n\
10040.0, 1.0, 0.0, 10.0\n\
10051.0, 0.0, 10.0, 0.0\n\
10060.0, 1.0, 0.0, 10.0\n\
10071.0, 0.0, 10.0, 0.0\n\
10080.0, 1.0, 0.0, 10.0\n\
10091.0, 0.0, 10.0, 0.0\n\
10100.0, 1.0, 0.0, 10.0\n\
10111.0, 0.0, 10.0, 0.0\n\
10120.0, 1.0, 0.0, 10.0\n\
10131.0, 0.0, 10.0, 0.0\n\
10140.0, 1.0, 0.0, 10.0\n\
10151.0, 0.0, 10.0, 0.0\n\
10160.0, 1.0, 0.0, 10.0\n\
10171.0, 0.0, 10.0, 0.0\n\
10180.0, 1.0, 0.0, 10.0\n\
10191.0, 0.0, 10.0, 0.0\n\
10200.0, 1.0, 0.0, 10.0\n\
10211.0, 0.0, 10.0, 0.0\n\
10220.0, 1.0, 0.0, 10.0\n\
10231.0, 0.0, 10.0, 0.0\n\
10240.0, 1.0, 0.0, 10.0\n\
10251.0, 0.0, 10.0, 0.0\n\
10260.0, 1.0, 0.0, 10.0\n\
10271.0, 0.0, 10.0, 0.0\n\
10280.0, 1.0, 0.0, 10.0\n\
10291.0, 0.0, 10.0, 0.0\n\
10300.0, 1.0, 0.0, 10.0\n\
10311.0, 0.0, 10.0, 0.0\n\
10320.0, 1.0, 0.0, 10.0\n\
10331.0, 0.0, 10.0, 0.0\n\
10340.0, 1.0, 0.0, 10.0\n\
10351.0, 0.0, 10.0, 0.0\n\
10360.0, 1.0, 0.0, 10.0\n\
10371.0, 0.0, 10.0, 0.0\n\
10380.0, 1.0, 0.0, 10.0\n\
10391.0, 0.0, 10.0, 0.0\n\
10400.0, 1.0, 0.0, 10.0\n\
10411.0, 0.0, 10.0, 0.0\n\
10420.0, 1.0, 0.0, 10.0\n\
10431.0, 0.0, 10.0, 0.0\n\
10440.0, 1.0, 0.0, 10.0\n\
10451.0, 0.0, 10.0, 0.0\n\
10460.0, 1.0, 0.0, 10.0\n\
10471.0, 0.0 ,10.4, 40.0\n")
1048        file.close()
1049
1050        pts_file = tempfile.mktemp(".pts")
1051       
1052        convert = Geospatial_data(fileName)
1053        convert.export_points_file(pts_file)
1054        results = Geospatial_data(pts_file, max_read_lines=2, verbose=True)
1055
1056        # Blocking
1057        geo_list = []
1058        for i in results:
1059            geo_list.append(i) 
1060        assert allclose(geo_list[0].get_data_points(),
1061                        [[1.0, 0.0],[0.0, 1.0]])
1062        assert allclose(geo_list[0].get_attributes(attribute_name='elevation'),
1063                        [10.0, 0.0])
1064        assert allclose(geo_list[1].get_data_points(),
1065                        [[1.0, 0.0],[0.0, 1.0] ])       
1066        assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
1067                        [10.0, 0.0])
1068           
1069        os.remove(fileName) 
1070        os.remove(pts_file)
1071       
1072       
1073
1074    def test_new_export_pts_file(self):
1075        att_dict = {}
1076        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1077        att_dict['elevation'] = array([10.1, 0.0, 10.4])
1078        att_dict['brightness'] = array([10.0, 1.0, 10.4])
1079       
1080        fileName = tempfile.mktemp(".pts")
1081       
1082        G = Geospatial_data(pointlist, att_dict)
1083       
1084        G.export_points_file(fileName)
1085
1086        results = Geospatial_data(file_name = fileName)
1087
1088        os.remove(fileName)
1089       
1090        assert allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1091        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
1092        answer = [10.0, 1.0, 10.4]
1093        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
1094
1095    def test_new_export_absolute_pts_file(self):
1096        att_dict = {}
1097        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1098        att_dict['elevation'] = array([10.1, 0.0, 10.4])
1099        att_dict['brightness'] = array([10.0, 1.0, 10.4])
1100        geo_ref = Geo_reference(50, 25, 55)
1101       
1102        fileName = tempfile.mktemp(".pts")
1103       
1104        G = Geospatial_data(pointlist, att_dict, geo_ref)
1105       
1106        G.export_points_file(fileName, absolute=True)
1107
1108        results = Geospatial_data(file_name = fileName)
1109
1110        os.remove(fileName)
1111       
1112        assert allclose(results.get_data_points(), G.get_data_points(True))
1113        assert allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
1114        answer = [10.0, 1.0, 10.4]
1115        assert allclose(results.get_attributes(attribute_name='brightness'), answer)
1116
1117    def test_loadpts(self):
1118       
1119        from Scientific.IO.NetCDF import NetCDFFile
1120
1121        fileName = tempfile.mktemp(".pts")
1122        # NetCDF file definition
1123        outfile = NetCDFFile(fileName, 'w')
1124       
1125        # dimension definitions
1126        outfile.createDimension('number_of_points', 3)   
1127        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1128   
1129        # variable definitions
1130        outfile.createVariable('points', Float, ('number_of_points',
1131                                                 'number_of_dimensions'))
1132        outfile.createVariable('elevation', Float, ('number_of_points',))
1133   
1134        # Get handles to the variables
1135        points = outfile.variables['points']
1136        elevation = outfile.variables['elevation']
1137 
1138        points[0, :] = [1.0,0.0]
1139        elevation[0] = 10.0 
1140        points[1, :] = [0.0,1.0]
1141        elevation[1] = 0.0 
1142        points[2, :] = [1.0,0.0]
1143        elevation[2] = 10.4   
1144
1145        outfile.close()
1146       
1147        results = Geospatial_data(file_name = fileName)
1148        os.remove(fileName)
1149        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
1150        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1151        assert allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
1152       
1153    def test_writepts(self):
1154        #test_writepts: Test that storage of x,y,attributes works
1155       
1156        att_dict = {}
1157        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1158        att_dict['elevation'] = array([10.0, 0.0, 10.4])
1159        att_dict['brightness'] = array([10.0, 0.0, 10.4])
1160        geo_reference=Geo_reference(56,1.9,1.9)
1161
1162        # Test pts format
1163        fileName = tempfile.mktemp(".pts")
1164        G = Geospatial_data(pointlist, att_dict, geo_reference)
1165        G.export_points_file(fileName, False)
1166        results = Geospatial_data(file_name=fileName)
1167        os.remove(fileName)
1168
1169        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1170        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1171        answer = [10.0, 0.0, 10.4]
1172        assert allclose(results.get_attributes('brightness'), answer)
1173        self.failUnless(geo_reference == geo_reference,
1174                         'test_writepts failed. Test geo_reference')
1175
1176    def test_write_csv_attributes(self):
1177        #test_write : Test that storage of x,y,attributes works
1178       
1179        att_dict = {}
1180        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1181        att_dict['elevation'] = array([10.0, 0.0, 10.4])
1182        att_dict['brightness'] = array([10.0, 0.0, 10.4])
1183        geo_reference=Geo_reference(56,0,0)
1184        # Test txt format
1185        fileName = tempfile.mktemp(".txt")
1186        G = Geospatial_data(pointlist, att_dict, geo_reference)
1187        G.export_points_file(fileName)
1188        #print "fileName", fileName
1189        results = Geospatial_data(file_name=fileName)
1190        os.remove(fileName)
1191        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1192        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1193        answer = [10.0, 0.0, 10.4]
1194        assert allclose(results.get_attributes('brightness'), answer)
1195       
1196 
1197    def test_write_csv_attributes_lat_long(self):
1198        #test_write : Test that storage of x,y,attributes works
1199       
1200        att_dict = {}
1201        pointlist = array([[-21.5,114.5],[-21.6,114.5],[-21.7,114.5]])
1202        att_dict['elevation'] = array([10.0, 0.0, 10.4])
1203        att_dict['brightness'] = array([10.0, 0.0, 10.4])
1204        # Test txt format
1205        fileName = tempfile.mktemp(".txt")
1206        G = Geospatial_data(pointlist, att_dict, points_are_lats_longs=True)
1207        G.export_points_file(fileName, as_lat_long=True)
1208        #print "fileName", fileName
1209        results = Geospatial_data(file_name=fileName)
1210        os.remove(fileName)
1211        assert allclose(results.get_data_points(False, as_lat_long=True),
1212                        pointlist)
1213        assert allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
1214        answer = [10.0, 0.0, 10.4]
1215        assert allclose(results.get_attributes('brightness'), answer)
1216       
1217    def test_writepts_no_attributes(self):
1218
1219        #test_writepts_no_attributes: Test that storage of x,y alone works
1220       
1221        att_dict = {}
1222        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1223        geo_reference=Geo_reference(56,1.9,1.9)
1224
1225        # Test pts format
1226        fileName = tempfile.mktemp(".pts")
1227        G = Geospatial_data(pointlist, None, geo_reference)
1228        G.export_points_file(fileName, False)
1229        results = Geospatial_data(file_name=fileName)
1230        os.remove(fileName)
1231
1232        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1233        self.failUnless(geo_reference == geo_reference,
1234                         'test_writepts failed. Test geo_reference')
1235       
1236       
1237    def test_write_csv_no_attributes(self):
1238        #test_write txt _no_attributes: Test that storage of x,y alone works
1239       
1240        att_dict = {}
1241        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1242        geo_reference=Geo_reference(56,0,0)
1243        # Test format
1244        fileName = tempfile.mktemp(".txt")
1245        G = Geospatial_data(pointlist, None, geo_reference)
1246        G.export_points_file(fileName)
1247        results = Geospatial_data(file_name=fileName)
1248        os.remove(fileName)
1249        assert allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1250
1251         
1252       
1253 ########################## BAD .PTS ##########################         
1254
1255    def test_load_bad_no_file_pts(self):
1256        import os
1257        import tempfile
1258       
1259        fileName = tempfile.mktemp(".pts")
1260        #print fileName
1261        try:
1262            results = Geospatial_data(file_name = fileName)
1263#            dict = import_points_file(fileName)
1264        except IOError:
1265            pass
1266        else:
1267            msg = 'imaginary file did not raise error!'
1268            raise msg
1269#            self.failUnless(0 == 1,
1270#                        'imaginary file did not raise error!')
1271
1272
1273    def test_create_from_pts_file(self):
1274       
1275        from Scientific.IO.NetCDF import NetCDFFile
1276
1277#        fileName = tempfile.mktemp(".pts")
1278        FN = 'test_points.pts'
1279        # NetCDF file definition
1280        outfile = NetCDFFile(FN, 'w')
1281       
1282        # dimension definitions
1283        outfile.createDimension('number_of_points', 3)   
1284        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1285   
1286        # variable definitions
1287        outfile.createVariable('points', Float, ('number_of_points',
1288                                                 'number_of_dimensions'))
1289        outfile.createVariable('elevation', Float, ('number_of_points',))
1290   
1291        # Get handles to the variables
1292        points = outfile.variables['points']
1293        elevation = outfile.variables['elevation']
1294 
1295        points[0, :] = [1.0,0.0]
1296        elevation[0] = 10.0 
1297        points[1, :] = [0.0,1.0]
1298        elevation[1] = 0.0 
1299        points[2, :] = [1.0,0.0]
1300        elevation[2] = 10.4   
1301
1302        outfile.close()
1303
1304        G = Geospatial_data(file_name = FN)
1305
1306        assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
1307        assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
1308
1309        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1310        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1311        os.remove(FN)
1312
1313    def test_create_from_pts_file_with_geo(self):
1314        """This test reveals if Geospatial data is correctly instantiated from a pts file.
1315        """
1316       
1317        from Scientific.IO.NetCDF import NetCDFFile
1318
1319        FN = 'test_points.pts'
1320        # NetCDF file definition
1321        outfile = NetCDFFile(FN, 'w')
1322
1323        # Make up an arbitrary georef
1324        xll = 0.1
1325        yll = 20
1326        geo_reference=Geo_reference(56, xll, yll)
1327        geo_reference.write_NetCDF(outfile)
1328
1329        # dimension definitions
1330        outfile.createDimension('number_of_points', 3)   
1331        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1332   
1333        # variable definitions
1334        outfile.createVariable('points', Float, ('number_of_points',
1335                                                 'number_of_dimensions'))
1336        outfile.createVariable('elevation', Float, ('number_of_points',))
1337   
1338        # Get handles to the variables
1339        points = outfile.variables['points']
1340        elevation = outfile.variables['elevation']
1341
1342        points[0, :] = [1.0,0.0]
1343        elevation[0] = 10.0 
1344        points[1, :] = [0.0,1.0]
1345        elevation[1] = 0.0 
1346        points[2, :] = [1.0,0.0]
1347        elevation[2] = 10.4   
1348
1349        outfile.close()
1350
1351        G = Geospatial_data(file_name = FN)
1352
1353        assert allclose(G.get_geo_reference().get_xllcorner(), xll)
1354        assert allclose(G.get_geo_reference().get_yllcorner(), yll)
1355
1356        assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
1357                                              [0.0+xll, 1.0+yll],
1358                                              [1.0+xll, 0.0+yll]])
1359       
1360        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1361        os.remove(FN)
1362
1363       
1364    def test_add_(self):
1365        '''test_add_(self):
1366        adds an txt and pts files, reads the files and adds them
1367           checking results are correct
1368        '''
1369        # create files
1370        att_dict1 = {}
1371        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1372        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
1373        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
1374        geo_reference1 = Geo_reference(56, 2.0, 1.0)
1375       
1376        att_dict2 = {}
1377        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
1378        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
1379        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
1380        geo_reference2 = Geo_reference(56, 1.0, 2.0) 
1381
1382        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1383        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1384       
1385        fileName1 = tempfile.mktemp(".txt")
1386        fileName2 = tempfile.mktemp(".pts")
1387
1388        #makes files
1389        G1.export_points_file(fileName1)
1390        G2.export_points_file(fileName2)
1391       
1392        # add files
1393       
1394        G3 = Geospatial_data(file_name = fileName1)
1395        G4 = Geospatial_data(file_name = fileName2)
1396       
1397        G = G3 + G4
1398
1399       
1400        #read results
1401#        print'res', G.get_data_points()
1402#        print'res1', G.get_data_points(False)
1403        assert allclose(G.get_data_points(),
1404                        [[ 3.0, 1.0], [ 2.0, 2.0],
1405                         [ 3.0, 1.0], [ 3.0, 3.0],
1406                         [ 2.0, 4.0], [ 3.0, 3.0]])
1407                         
1408        assert allclose(G.get_attributes(attribute_name='elevation'),
1409                        [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
1410       
1411        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
1412        assert allclose(G.get_attributes(attribute_name='brightness'), answer)
1413       
1414        self.failUnless(G.get_geo_reference() == geo_reference1,
1415                         'test_writepts failed. Test geo_reference')
1416                         
1417        os.remove(fileName1)
1418        os.remove(fileName2)
1419       
1420    def test_ensure_absolute(self):
1421        points = [[2.0, 0.0],[1.0, 1.0],
1422                         [2.0, 0.0],[2.0, 2.0],
1423                         [1.0, 3.0],[2.0, 2.0]]
1424        new_points = ensure_absolute(points)
1425       
1426        assert allclose(new_points, points)
1427
1428        points = array([[2.0, 0.0],[1.0, 1.0],
1429                         [2.0, 0.0],[2.0, 2.0],
1430                         [1.0, 3.0],[2.0, 2.0]])
1431        new_points = ensure_absolute(points)
1432       
1433        assert allclose(new_points, points)
1434       
1435        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1436                         [2.0, 0.0],[2.0, 2.0],
1437                         [1.0, 3.0],[2.0, 2.0]])
1438       
1439        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1440
1441        data_points = zeros((ab_points.shape), Float)
1442        #Shift datapoints according to new origins
1443        for k in range(len(ab_points)):
1444            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1445            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1446        #print "data_points",data_points     
1447        new_points = ensure_absolute(data_points,
1448                                             geo_reference=mesh_origin)
1449        #print "new_points",new_points
1450        #print "ab_points",ab_points
1451           
1452        assert allclose(new_points, ab_points)
1453
1454        geo = Geo_reference(56,67,-56)
1455
1456        data_points = geo.change_points_geo_ref(ab_points)   
1457        new_points = ensure_absolute(data_points,
1458                                             geo_reference=geo)
1459        #print "new_points",new_points
1460        #print "ab_points",ab_points
1461           
1462        assert allclose(new_points, ab_points)
1463
1464
1465        geo_reference = Geo_reference(56, 100, 200)
1466        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1467        points = geo_reference.change_points_geo_ref(ab_points)
1468        attributes = [2, 4]
1469        #print "geo in points", points
1470        G = Geospatial_data(points, attributes,
1471                            geo_reference=geo_reference)
1472         
1473        new_points = ensure_absolute(G)
1474        #print "new_points",new_points
1475        #print "ab_points",ab_points
1476           
1477        assert allclose(new_points, ab_points)
1478
1479
1480       
1481    def test_ensure_geospatial(self):
1482        points = [[2.0, 0.0],[1.0, 1.0],
1483                         [2.0, 0.0],[2.0, 2.0],
1484                         [1.0, 3.0],[2.0, 2.0]]
1485        new_points = ensure_geospatial(points)
1486       
1487        assert allclose(new_points.get_data_points(absolute = True), points)
1488
1489        points = array([[2.0, 0.0],[1.0, 1.0],
1490                         [2.0, 0.0],[2.0, 2.0],
1491                         [1.0, 3.0],[2.0, 2.0]])
1492        new_points = ensure_geospatial(points)
1493       
1494        assert allclose(new_points.get_data_points(absolute = True), points)
1495       
1496        ab_points = array([[2.0, 0.0],[1.0, 1.0],
1497                         [2.0, 0.0],[2.0, 2.0],
1498                         [1.0, 3.0],[2.0, 2.0]])
1499       
1500        mesh_origin = (56, 290000, 618000) #zone, easting, northing
1501
1502        data_points = zeros((ab_points.shape), Float)
1503        #Shift datapoints according to new origins
1504        for k in range(len(ab_points)):
1505            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1506            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1507        #print "data_points",data_points     
1508        new_geospatial = ensure_geospatial(data_points,
1509                                             geo_reference=mesh_origin)
1510        new_points = new_geospatial.get_data_points(absolute=True)
1511        #print "new_points",new_points
1512        #print "ab_points",ab_points
1513           
1514        assert allclose(new_points, ab_points)
1515
1516        geo = Geo_reference(56,67,-56)
1517
1518        data_points = geo.change_points_geo_ref(ab_points)   
1519        new_geospatial = ensure_geospatial(data_points,
1520                                             geo_reference=geo)
1521        new_points = new_geospatial.get_data_points(absolute=True)
1522        #print "new_points",new_points
1523        #print "ab_points",ab_points
1524           
1525        assert allclose(new_points, ab_points)
1526
1527
1528        geo_reference = Geo_reference(56, 100, 200)
1529        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1530        points = geo_reference.change_points_geo_ref(ab_points)
1531        attributes = [2, 4]
1532        #print "geo in points", points
1533        G = Geospatial_data(points, attributes,
1534                            geo_reference=geo_reference)
1535         
1536        new_geospatial  = ensure_geospatial(G)
1537        new_points = new_geospatial.get_data_points(absolute=True)
1538        #print "new_points",new_points
1539        #print "ab_points",ab_points
1540           
1541        assert allclose(new_points, ab_points)
1542       
1543    def test_isinstance(self):
1544
1545        import os
1546       
1547        fileName = tempfile.mktemp(".csv")
1548        file = open(fileName,"w")
1549        file.write("x,y,  elevation ,  speed \n\
15501.0, 0.0, 10.0, 0.0\n\
15510.0, 1.0, 0.0, 10.0\n\
15521.0, 0.0, 10.4, 40.0\n")
1553        file.close()
1554
1555        results = Geospatial_data(fileName)
1556        assert allclose(results.get_data_points(absolute=True), \
1557                        [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
1558        assert allclose(results.get_attributes(attribute_name='elevation'), \
1559                        [10.0, 0.0, 10.4])
1560        assert allclose(results.get_attributes(attribute_name='speed'), \
1561                        [0.0, 10.0, 40.0])
1562
1563        os.remove(fileName)
1564       
1565
1566    def test_no_constructors(self):
1567       
1568        try:
1569            G = Geospatial_data()
1570#            results = Geospatial_data(file_name = fileName)
1571#            dict = import_points_file(fileName)
1572        except ValueError:
1573            pass
1574        else:
1575            msg = 'Instance must have a filename or data points'
1576            raise msg       
1577
1578    def test_load_csv_lat_long(self):
1579        """
1580        comma delimited
1581
1582        """
1583        fileName = tempfile.mktemp(".csv")
1584        file = open(fileName,"w")
1585        file.write("long,lat, elevation, yeah \n\
1586150.916666667,-34.50,452.688000, 10\n\
1587150.0,-34,459.126000, 10\n")
1588        file.close()
1589        results = Geospatial_data(fileName)
1590        os.remove(fileName)
1591        points = results.get_data_points()
1592       
1593        assert allclose(points[0][0], 308728.009)
1594        assert allclose(points[0][1], 6180432.601)
1595        assert allclose(points[1][0],  222908.705)
1596        assert allclose(points[1][1], 6233785.284)
1597       
1598     
1599    def test_load_csv_lat_longII(self):
1600        """
1601        comma delimited
1602
1603        """
1604        fileName = tempfile.mktemp(".csv")
1605        file = open(fileName,"w")
1606        file.write("Lati,LONG,z \n\
1607-34.50,150.916666667,452.688000\n\
1608-34,150.0,459.126000\n")
1609        file.close()
1610        results = Geospatial_data(fileName)
1611        os.remove(fileName)
1612        points = results.get_data_points()
1613       
1614        assert allclose(points[0][0], 308728.009)
1615        assert allclose(points[0][1], 6180432.601)
1616        assert allclose(points[1][0],  222908.705)
1617        assert allclose(points[1][1], 6233785.284)
1618
1619         
1620    def test_load_csv_lat_long_bad(self):
1621        """
1622        comma delimited
1623
1624        """
1625        fileName = tempfile.mktemp(".csv")
1626        file = open(fileName,"w")
1627        file.write("Lati,LONG,z \n\
1628-25.0,180.0,452.688000\n\
1629-34,150.0,459.126000\n")
1630        file.close()
1631        try:
1632            results = Geospatial_data(fileName)
1633        except ANUGAError:
1634            pass
1635        else:
1636            msg = 'Different zones in Geo references not caught.'
1637            raise msg       
1638       
1639        os.remove(fileName)
1640       
1641    def test_lat_long(self):
1642        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1643        lon_gong = degminsec2decimal_degrees(150,55,0.)
1644       
1645        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1646        lon_2 = degminsec2decimal_degrees(150,00,0.)
1647       
1648        lats = [lat_gong, lat_2]
1649        longs = [lon_gong, lon_2]
1650        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1651
1652        points = gsd.get_data_points(absolute=True)
1653       
1654        assert allclose(points[0][0], 308728.009)
1655        assert allclose(points[0][1], 6180432.601)
1656        assert allclose(points[1][0],  222908.705)
1657        assert allclose(points[1][1], 6233785.284)
1658        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1659                        'Bad zone error!')
1660       
1661        try:
1662            results = Geospatial_data(latitudes=lats)
1663        except ValueError:
1664            pass
1665        else:
1666            self.failUnless(0 ==1,  'Error not thrown error!')
1667        try:
1668            results = Geospatial_data(latitudes=lats)
1669        except ValueError:
1670            pass
1671        else:
1672            self.failUnless(0 ==1,  'Error not thrown error!')
1673        try:
1674            results = Geospatial_data(longitudes=lats)
1675        except ValueError:
1676            pass
1677        else:
1678            self.failUnless(0 ==1, 'Error not thrown error!')
1679        try:
1680            results = Geospatial_data(latitudes=lats, longitudes=longs,
1681                                      geo_reference="p")
1682        except ValueError:
1683            pass
1684        else:
1685            self.failUnless(0 ==1,  'Error not thrown error!')
1686           
1687        try:
1688            results = Geospatial_data(latitudes=lats, longitudes=longs,
1689                                      data_points=12)
1690        except ValueError:
1691            pass
1692        else:
1693            self.failUnless(0 ==1,  'Error not thrown error!')
1694
1695    def test_lat_long2(self):
1696        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1697        lon_gong = degminsec2decimal_degrees(150,55,0.)
1698       
1699        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1700        lon_2 = degminsec2decimal_degrees(150,00,0.)
1701       
1702        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1703        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1704
1705        points = gsd.get_data_points(absolute=True)
1706       
1707        assert allclose(points[0][0], 308728.009)
1708        assert allclose(points[0][1], 6180432.601)
1709        assert allclose(points[1][0],  222908.705)
1710        assert allclose(points[1][1], 6233785.284)
1711        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1712                        'Bad zone error!')
1713
1714        try:
1715            results = Geospatial_data(points_are_lats_longs=True)
1716        except ValueError:
1717            pass
1718        else:
1719            self.failUnless(0 ==1,  'Error not thrown error!')
1720
1721
1722    def test_write_urs_file(self):
1723        lat_gong = degminsec2decimal_degrees(-34,00,0)
1724        lon_gong = degminsec2decimal_degrees(150,30,0.)
1725       
1726        lat_2 = degminsec2decimal_degrees(-34,00,1)
1727        lon_2 = degminsec2decimal_degrees(150,00,0.)
1728        p1 = (lat_gong, lon_gong)
1729        p2 = (lat_2, lon_2)
1730        points = ImmutableSet([p1, p2, p1])
1731        gsd = Geospatial_data(data_points=list(points),
1732                              points_are_lats_longs=True)
1733       
1734        fn = tempfile.mktemp(".urs")
1735        gsd.export_points_file(fn)
1736        #print "fn", fn
1737        handle = open(fn)
1738        lines = handle.readlines()
1739        assert lines[0],'2'
1740        assert lines[1],'-34.0002778 150.0 0'
1741        assert lines[2],'-34.0 150.5 1'
1742        handle.close()
1743        os.remove(fn)
1744       
1745    def test_lat_long_set(self):
1746        lat_gong = degminsec2decimal_degrees(-34,30,0.)
1747        lon_gong = degminsec2decimal_degrees(150,55,0.)
1748       
1749        lat_2 = degminsec2decimal_degrees(-34,00,0.)
1750        lon_2 = degminsec2decimal_degrees(150,00,0.)
1751        p1 = (lat_gong, lon_gong)
1752        p2 = (lat_2, lon_2)
1753        points = ImmutableSet([p1, p2, p1])
1754        gsd = Geospatial_data(data_points=list(points),
1755                              points_are_lats_longs=True)
1756
1757        points = gsd.get_data_points(absolute=True)
1758        #print "points[0][0]", points[0][0]
1759        #Note the order is unknown, due to using sets
1760        # and it changes from windows to linux
1761        try:
1762            assert allclose(points[1][0], 308728.009)
1763            assert allclose(points[1][1], 6180432.601)
1764            assert allclose(points[0][0],  222908.705)
1765            assert allclose(points[0][1], 6233785.284)
1766        except AssertionError:
1767            assert allclose(points[0][0], 308728.009)
1768            assert allclose(points[0][1], 6180432.601)
1769            assert allclose(points[1][0],  222908.705)
1770            assert allclose(points[1][1], 6233785.284)
1771           
1772        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1773                        'Bad zone error!')
1774        points = gsd.get_data_points(as_lat_long=True)
1775        #print "test_lat_long_set points", points
1776        try:
1777            assert allclose(points[0][0], -34)
1778            assert allclose(points[0][1], 150)
1779        except AssertionError:
1780            assert allclose(points[1][0], -34)
1781            assert allclose(points[1][1], 150)
1782
1783    def test_len(self):
1784       
1785        points = [[1.0, 2.1], [3.0, 5.3]]
1786        G = Geospatial_data(points)
1787        self.failUnless(2 ==len(G),  'Len error!')
1788       
1789        points = [[1.0, 2.1]]
1790        G = Geospatial_data(points)
1791        self.failUnless(1 ==len(G),  'Len error!')
1792
1793        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1794        G = Geospatial_data(points)
1795        self.failUnless(4 ==len(G),  'Len error!')
1796       
1797    def test_split(self):
1798        """test if the results from spilt are disjoin sets"""
1799       
1800        #below is a work around until the randint works on cyclones compute nodes
1801        if get_host_name()[8:9]!='0':
1802               
1803           
1804            points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
1805                      [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
1806                      [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0],
1807                      [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
1808                      [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
1809            attributes = {'depth':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
1810                          14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25],
1811                          'speed':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
1812                          14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]}
1813            G = Geospatial_data(points, attributes)
1814   
1815            factor = 0.21
1816   
1817            #will return G1 with 10% of points and G2 with 90%
1818            G1, G2  = G.split(factor,100) 
1819           
1820            assert allclose(len(G), len(G1)+len(G2))
1821            assert allclose(round(len(G)*factor), len(G1))
1822   
1823            P = G1.get_data_points(absolute=False)
1824            assert allclose(P, [[5.0,4.0],[4.0,3.0],[4.0,2.0],[3.0,1.0],[2.0,3.0]])
1825   
1826            A = G1.get_attributes()
1827            assert allclose(A,[24, 18, 17, 11, 8])
1828       
1829    def test_split1(self):
1830        """test if the results from spilt are disjoin sets"""
1831        #below is a work around until the randint works on cyclones compute nodes
1832        if get_host_name()[8:9]!='0':
1833
1834            from RandomArray import randint,seed
1835            seed(100,100)
1836            a_points = randint(0,999999,(10,2))
1837            points = a_points.tolist()
1838    #        print points
1839   
1840            G = Geospatial_data(points)
1841   
1842            factor = 0.1
1843   
1844            #will return G1 with 10% of points and G2 with 90%
1845            G1, G2  = G.split(factor,100) 
1846           
1847    #        print 'G1',G1
1848            assert allclose(len(G), len(G1)+len(G2))
1849            assert allclose(round(len(G)*factor), len(G1))
1850   
1851            P = G1.get_data_points(absolute=False)
1852            assert allclose(P, [[982420.,28233.]])
1853
1854 
1855    def test_find_optimal_smoothing_parameter(self):
1856        """
1857        Creates a elevation file represting hill (sort of) and runs
1858        find_optimal_smoothing_parameter for 3 different alphas,
1859       
1860        NOTE the random number seed is provided to control the results
1861        """
1862        from cmath import cos
1863
1864        #below is a work around until the randint works on cyclones compute nodes
1865        if get_host_name()[8:9]!='0':
1866
1867            filename = tempfile.mktemp(".csv")
1868            file = open(filename,"w")
1869            file.write("x,y,elevation \n")
1870   
1871            for i in range(-5,6):
1872                for j in range(-5,6):
1873                    #this equation made surface like a circle ripple
1874                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1875    #                print 'x,y,f',i,j,z
1876                    file.write("%s, %s, %s\n" %(i, j, z))
1877                   
1878            file.close()
1879     
1880            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
1881                                                 alpha_list=[0.0001, 0.01, 1],
1882                                                 mesh_file=None,
1883                                                 mesh_resolution=3,
1884                                                 north_boundary=5,
1885                                                 south_boundary=-5,
1886                                                 east_boundary=5,
1887                                                 west_boundary=-5,
1888                                                 plot_name=None,
1889                                                 seed_num=100000,
1890                                                 verbose=False)
1891   
1892            os.remove(filename)
1893           
1894            # print value, alpha
1895            assert (alpha==0.01)
1896
1897    def test_find_optimal_smoothing_parameter1(self):
1898        """
1899        Creates a elevation file represting hill (sort of) and
1900        Then creates a mesh file and passes the mesh file and the elevation
1901        file to find_optimal_smoothing_parameter for 3 different alphas,
1902       
1903        NOTE the random number seed is provided to control the results
1904        """
1905        #below is a work around until the randint works on cyclones compute nodes
1906        if get_host_name()[8:9]!='0':
1907
1908            from cmath import cos
1909            from anuga.pmesh.mesh_interface import create_mesh_from_regions
1910           
1911            filename = tempfile.mktemp(".csv")
1912            file = open(filename,"w")
1913            file.write("x,y,elevation \n")
1914   
1915            for i in range(-5,6):
1916                for j in range(-5,6):
1917                    #this equation made surface like a circle ripple
1918                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1919    #                print 'x,y,f',i,j,z
1920                    file.write("%s, %s, %s\n" %(i, j, z))
1921                   
1922            file.close()
1923            poly=[[5,5],[5,-5],[-5,-5],[-5,5]]
1924            internal_poly=[[[[1,1],[1,-1],[-1,-1],[-1,1]],.5]]
1925            mesh_filename= tempfile.mktemp(".msh")
1926           
1927            create_mesh_from_regions(poly,
1928                                 boundary_tags={'back': [2],
1929                                                'side': [1,3],
1930                                                'ocean': [0]},
1931                             maximum_triangle_area=3,
1932                             interior_regions=internal_poly,
1933                             filename=mesh_filename,
1934                             use_cache=False,
1935                             verbose=False)
1936     
1937            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
1938                                                 alpha_list=[0.0001, 0.01, 1],
1939                                                 mesh_file=mesh_filename,
1940                                                 plot_name=None,
1941                                                 seed_num=174,
1942                                                 verbose=False)
1943   
1944            os.remove(filename)
1945            os.remove(mesh_filename)
1946           
1947    #        print value, alpha
1948            assert (alpha==0.01)
1949
1950    def test_find_optimal_smoothing_parameter2(self):
1951        """
1952        Tests requirement that mesh file must exist or IOError is thrown
1953       
1954        NOTE the random number seed is provided to control the results
1955        """
1956        from cmath import cos
1957        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1958       
1959        filename = tempfile.mktemp(".csv")
1960        mesh_filename= tempfile.mktemp(".msh")
1961       
1962        try:
1963            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
1964                                             alpha_list=[0.0001, 0.01, 1],
1965                                             mesh_file=mesh_filename,
1966                                             plot_name=None,
1967                                             seed_num=174,
1968                                             verbose=False)
1969        except IOError:
1970            pass
1971        else:
1972            self.failUnless(0 ==1,  'Error not thrown error!')
1973       
1974         
1975if __name__ == "__main__":
1976
1977    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')
1978    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
1979    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
1980    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1981    runner = unittest.TextTestRunner() #verbosity=2)
1982    runner.run(suite)
1983
1984   
Note: See TracBrowser for help on using the repository browser.