source: branches/numpy/anuga/geospatial_data/test_geospatial_data.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 15 years ago

Initial commit of numpy changes. Still a long way to go.

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