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

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

numpy changes.

File size: 67.7 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=\n%s\nbut G.get_data_points(absolute=True)=\n%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=\n%s\nshould be close to\n'
345               '[[2.0, 4.1], [4.0, 7.3],\n'
346               ' [5.1, 9.1], [6.1, 6.3]]'
347               % str(P))
348        assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3],
349                                [5.1, 9.1], [6.1, 6.3]]), msg
350
351    def test_add_with_geo_absolute (self):
352        '''Difference in Geo_reference resolved'''
353
354        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
355        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
356        attributes1 = [2, 4]
357        attributes2 = [5, 76]
358        geo_ref1= Geo_reference(55, 1.0, 2.0)
359        geo_ref2 = Geo_reference(55, 2.0, 3.0)
360
361        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(),
362                                        geo_ref1.get_yllcorner()],
363                             attributes1, geo_ref1)
364
365        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(),
366                                        geo_ref2.get_yllcorner()],
367                             attributes2, geo_ref2)
368
369        #Check that absolute values are as expected
370        P1 = G1.get_data_points(absolute=True)
371        assert num.allclose(P1, points1)
372
373        P1 = G1.get_data_points(absolute=False)
374        msg = ('P1=\n%s\nbut points1 - <...>=\n%s'
375               % (str(P1),
376                  str(points1 - [geo_ref1.get_xllcorner(),
377                                 geo_ref1.get_yllcorner()])))
378        assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(),
379                                           geo_ref1.get_yllcorner()]), msg
380
381        P2 = G2.get_data_points(absolute=True)
382        assert num.allclose(P2, points2)
383
384        P2 = G2.get_data_points(absolute=False)
385        assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(),
386                                           geo_ref2.get_yllcorner()])
387
388        G = G1 + G2
389        P = G.get_data_points(absolute=True)
390
391        assert num.allclose(P, num.concatenate((points1, points2)))
392
393    def test_add_with_None(self):
394        '''test that None can be added to a geospatical objects'''
395
396        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
397        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
398
399        geo_ref1= Geo_reference(55, 1.0, 2.0)
400        geo_ref2 = Geo_reference(zone=55,
401                                 xllcorner=0.1,
402                                 yllcorner=3.0,
403                                 datum='wgs84',
404                                 projection='UTM',
405                                 units='m')
406
407        attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]}
408        attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]}
409
410        G1 = Geospatial_data(points1, attributes1, geo_ref1)
411        assert num.allclose(G1.get_geo_reference().get_xllcorner(), 1.0)
412        assert num.allclose(G1.get_geo_reference().get_yllcorner(), 2.0)
413        assert G1.attributes.has_key('depth')
414        assert G1.attributes.has_key('elevation')
415        assert num.allclose(G1.attributes['depth'], [2, 4.7])
416        assert num.allclose(G1.attributes['elevation'], [6.1, 5])
417
418        G2 = Geospatial_data(points2, attributes2, geo_ref2)
419        assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1)
420        assert num.allclose(G2.get_geo_reference().get_yllcorner(), 3.0)
421        assert G2.attributes.has_key('depth')
422        assert G2.attributes.has_key('elevation')
423        assert num.allclose(G2.attributes['depth'], [-2.3, 4])
424        assert num.allclose(G2.attributes['elevation'], [2.5, 1])
425
426        #Check that absolute values are as expected
427        P1 = G1.get_data_points(absolute=True)
428        assert num.allclose(P1, [[3.0, 6.1], [5.0, 9.3]])
429
430        P2 = G2.get_data_points(absolute=True)
431        assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])
432
433        # Normal add
434        G = G1 + None
435        assert G.attributes.has_key('depth')
436        assert G.attributes.has_key('elevation')
437        assert num.allclose(G.attributes['depth'], [2, 4.7])
438        assert num.allclose(G.attributes['elevation'], [6.1, 5])
439
440        # Points are now absolute.
441        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
442        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
443        P = G.get_data_points(absolute=True)
444        msg = 'P=\n%s' % str(P)
445        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg
446
447        G = G2 + None
448        assert G.attributes.has_key('depth')
449        assert G.attributes.has_key('elevation')
450        assert num.allclose(G.attributes['depth'], [-2.3, 4])
451        assert num.allclose(G.attributes['elevation'], [2.5, 1])
452        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
453        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
454
455        P = G.get_data_points(absolute=True)
456        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
457
458        # Reverse add
459        G = None + G1
460        assert G.attributes.has_key('depth')
461        assert G.attributes.has_key('elevation')
462        assert num.allclose(G.attributes['depth'], [2, 4.7])
463        assert num.allclose(G.attributes['elevation'], [6.1, 5])
464
465        # Points are now absolute.
466        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
467        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
468
469        P = G.get_data_points(absolute=True)
470        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
471
472        G = None + G2
473        assert G.attributes.has_key('depth')
474        assert G.attributes.has_key('elevation')
475        assert num.allclose(G.attributes['depth'], [-2.3, 4])
476        assert num.allclose(G.attributes['elevation'], [2.5, 1])
477
478        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
479        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
480
481        P = G.get_data_points(absolute=True)
482        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
483
484    def test_clip0(self):
485        '''test_clip0(self):
486
487        Test that point sets can be clipped by a polygon
488        '''
489
490        from anuga.coordinate_transforms.geo_reference import Geo_reference
491
492        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
493                  [0, 0], [2.4, 3.3]]
494        G = Geospatial_data(points)
495
496        # First try the unit square
497        U = [[0,0], [1,0], [1,1], [0,1]]
498        assert num.allclose(G.clip(U).get_data_points(),
499                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
500
501        # Then a more complex polygon
502        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
503        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
504                  [0.5, 1.5], [0.5, -0.5]]
505        G = Geospatial_data(points)
506
507        assert num.allclose(G.clip(polygon).get_data_points(),
508                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
509
510    def test_clip0_with_attributes(self):
511        '''test_clip0_with_attributes(self):
512
513        Test that point sets with attributes can be clipped by a polygon
514        '''
515
516        from anuga.coordinate_transforms.geo_reference import Geo_reference
517
518        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
519                  [0, 0], [2.4, 3.3]]
520
521        attributes = [2, -4, 5, 76, -2, 0.1, 3]
522        att_dict = {'att1': attributes,
523                    'att2': num.array(attributes)+1}
524
525        G = Geospatial_data(points, att_dict)
526
527        # First try the unit square
528        U = [[0,0], [1,0], [1,1], [0,1]]
529        assert num.allclose(G.clip(U).get_data_points(),
530                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
531        assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
532        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
533
534        # Then a more complex polygon
535        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
536        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
537                  [0.5, 1.5], [0.5, -0.5]]
538
539        # This time just one attribute
540        attributes = [2, -4, 5, 76, -2, 0.1]
541        G = Geospatial_data(points, attributes)
542        assert num.allclose(G.clip(polygon).get_data_points(),
543                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
544        assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
545
546    def test_clip1(self):
547        '''test_clip1(self):
548
549        Test that point sets can be clipped by a polygon given as
550        another Geospatial dataset
551        '''
552
553        from anuga.coordinate_transforms.geo_reference import Geo_reference
554
555        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
556                  [0, 0], [2.4, 3.3]]
557        attributes = [2, -4, 5, 76, -2, 0.1, 3]
558        att_dict = {'att1': attributes,
559                    'att2': num.array(attributes)+1}
560        G = Geospatial_data(points, att_dict)
561
562        # First try the unit square
563        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
564        assert num.allclose(G.clip(U).get_data_points(),
565                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
566        assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
567        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
568
569        # Then a more complex polygon
570        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
571                  [0.5, 1.5], [0.5, -0.5]]
572        attributes = [2, -4, 5, 76, -2, 0.1]
573        G = Geospatial_data(points, attributes)
574        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1],
575                                   [2,1], [0,1]])
576
577        assert num.allclose(G.clip(polygon).get_data_points(),
578                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
579        assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
580
581    def test_clip0_outside(self):
582        '''test_clip0_outside(self):
583
584        Test that point sets can be clipped outside of a polygon
585        '''
586
587        from anuga.coordinate_transforms.geo_reference import Geo_reference
588
589        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
590                  [0, 0], [2.4, 3.3]]
591        attributes = [2, -4, 5, 76, -2, 0.1, 3]
592        G = Geospatial_data(points, attributes)
593
594        # First try the unit square
595        U = [[0,0], [1,0], [1,1], [0,1]]
596        assert num.allclose(G.clip_outside(U).get_data_points(),
597                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
598        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
599
600        # Then a more complex polygon
601        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
602        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
603                  [0.5, 1.5], [0.5, -0.5]]
604        attributes = [2, -4, 5, 76, -2, 0.1]
605        G = Geospatial_data(points, attributes)
606        assert num.allclose(G.clip_outside(polygon).get_data_points(),
607                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
608        assert num.allclose(G.clip_outside(polygon).get_attributes(),
609                            [2, -2, 0.1])
610
611    def test_clip1_outside(self):
612        '''test_clip1_outside(self):
613
614        Test that point sets can be clipped outside of a polygon given as
615        another Geospatial dataset
616        '''
617
618        from anuga.coordinate_transforms.geo_reference import Geo_reference
619
620        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
621                  [0, 0], [2.4, 3.3]]
622        attributes = [2, -4, 5, 76, -2, 0.1, 3]
623        G = Geospatial_data(points, attributes)
624
625        # First try the unit square
626        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
627        assert num.allclose(G.clip_outside(U).get_data_points(),
628                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
629        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
630
631        # Then a more complex polygon
632        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
633                  [0.5, 1.5], [0.5, -0.5]]
634        attributes = [2, -4, 5, 76, -2, 0.1]
635        G = Geospatial_data(points, attributes)
636        polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1],
637                                   [2, 1], [0, 1]])
638        assert num.allclose(G.clip_outside(polygon).get_data_points(),
639                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
640        assert num.allclose(G.clip_outside(polygon).get_attributes(),
641                            [2, -2, 0.1])
642
643    def test_clip1_inside_outside(self):
644        '''test_clip1_inside_outside(self):
645
646        Test that point sets can be clipped outside of a polygon given as
647        another Geospatial dataset
648        '''
649
650        from anuga.coordinate_transforms.geo_reference import Geo_reference
651
652        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
653                  [0, 0], [2.4, 3.3]]
654        attributes = [2, -4, 5, 76, -2, 0.1, 3]
655        G = Geospatial_data(points, attributes)
656
657        # First try the unit square
658        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
659        G1 = G.clip(U)
660        assert num.allclose(G1.get_data_points(),
661                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
662        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
663        G2 = G.clip_outside(U)
664        assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
665                                                  [3.0, 5.3], [2.4, 3.3]])
666        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
667
668        # New ordering
669        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0], [-1, 4],
670                      [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
671        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]
672
673        assert num.allclose((G1+G2).get_data_points(), new_points)
674        assert num.allclose((G1+G2).get_attributes(), new_attributes)
675
676        G = G1+G2
677        FN = 'test_combine.pts'
678        G.export_points_file(FN)
679
680        # Read it back in
681        G3 = Geospatial_data(FN)
682
683        # Check result
684        assert num.allclose(G3.get_data_points(), new_points)
685        assert num.allclose(G3.get_attributes(), new_attributes)
686
687        os.remove(FN)
688
689    def test_load_csv(self):
690        fileName = tempfile.mktemp('.csv')
691        file = open(fileName,'w')
692        file.write('x,y,elevation speed \n\
6931.0 0.0 10.0 0.0\n\
6940.0 1.0 0.0 10.0\n\
6951.0 0.0 10.4 40.0\n')
696        file.close()
697
698        results = Geospatial_data(fileName)
699        os.remove(fileName)
700        assert num.allclose(results.get_data_points(),
701                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
702        assert num.allclose(results.get_attributes(attribute_name='elevation'),
703                            [10.0, 0.0, 10.4])
704        assert num.allclose(results.get_attributes(attribute_name='speed'),
705                            [0.0, 10.0, 40.0])
706
707################################################################################
708# Test CSV files
709################################################################################
710
711    def test_load_csv_lat_long_bad_blocking(self):
712        '''test_load_csv_lat_long_bad_blocking(self):
713        Different zones in Geo references
714        '''
715
716        fileName = tempfile.mktemp('.csv')
717        file = open(fileName, 'w')
718        file.write('Lati,LONG,z \n\
719-25.0,180.0,452.688000\n\
720-34,150.0,459.126000\n')
721        file.close()
722
723        results = Geospatial_data(fileName, max_read_lines=1,
724                                  load_file_now=False)
725
726        try:
727            for i in results:
728                pass
729        except ANUGAError:
730            pass
731        else:
732            msg = 'Different zones in Geo references not caught.'
733            raise Exception, msg
734
735        os.remove(fileName)
736
737    def test_load_csv(self):
738        fileName = tempfile.mktemp('.txt')
739        file = open(fileName, 'w')
740        file.write(' x,y, elevation ,  speed \n\
7411.0, 0.0, 10.0, 0.0\n\
7420.0, 1.0, 0.0, 10.0\n\
7431.0, 0.0 ,10.4, 40.0\n')
744        file.close()
745
746        results = Geospatial_data(fileName, max_read_lines=2)
747
748        assert num.allclose(results.get_data_points(),
749                            [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
750        assert num.allclose(results.get_attributes(attribute_name='elevation'),
751                            [10.0, 0.0, 10.4])
752        assert num.allclose(results.get_attributes(attribute_name='speed'),
753                            [0.0, 10.0, 40.0])
754
755        # Blocking
756        geo_list = []
757        for i in results:
758            geo_list.append(i)
759
760        assert num.allclose(geo_list[0].get_data_points(),
761                            [[1.0, 0.0],[0.0, 1.0]])
762        assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
763                            [10.0, 0.0])
764        assert num.allclose(geo_list[1].get_data_points(),
765                            [[1.0, 0.0]])
766        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
767                            [10.4])
768
769        os.remove(fileName)
770
771    def test_load_csv_bad(self):
772        '''test_load_csv_bad(self):
773        header column, body column missmatch
774        (Format error)
775        '''
776
777        fileName = tempfile.mktemp('.txt')
778        file = open(fileName, 'w')
779        file.write(' elevation ,  speed \n\
7801.0, 0.0, 10.0, 0.0\n\
7810.0, 1.0, 0.0, 10.0\n\
7821.0, 0.0 ,10.4, 40.0\n')
783        file.close()
784
785        results = Geospatial_data(fileName, max_read_lines=2,
786                                  load_file_now=False)
787
788        # Blocking
789        geo_list = []
790        try:
791            for i in results:
792                geo_list.append(i)
793        except SyntaxError:
794            pass
795        else:
796            msg = 'Bad file did not raise error!'
797            raise Exception, msg
798        os.remove(fileName)
799
800    def test_load_csv_badII(self):
801        '''test_load_csv_bad(self):
802        header column, body column missmatch
803        (Format error)
804        '''
805
806        fileName = tempfile.mktemp('.txt')
807        file = open(fileName, 'w')
808        file.write(' x,y,elevation ,  speed \n\
8091.0, 0.0, 10.0, 0.0\n\
8100.0, 1.0, 0.0, 10\n\
8111.0, 0.0 ,10.4 yeah\n')
812        file.close()
813
814        results = Geospatial_data(fileName, max_read_lines=2,
815                                  load_file_now=False)
816
817        # Blocking
818        geo_list = []
819        try:
820            for i in results:
821                geo_list.append(i)
822        except SyntaxError:
823            pass
824        else:
825            msg = 'Bad file did not raise error!'
826            raise Exception, msg
827        os.remove(fileName)
828
829    def test_load_csv_badIII(self):
830        '''test_load_csv_bad(self):
831        space delimited
832        '''
833
834        fileName = tempfile.mktemp('.txt')
835        file = open(fileName, 'w')
836        file.write(' x y elevation   speed \n\
8371. 0.0 10.0 0.0\n\
8380.0 1.0 0.0 10.0\n\
8391.0 0.0 10.4 40.0\n')
840        file.close()
841
842        try:
843            results = Geospatial_data(fileName, max_read_lines=2,
844                                      load_file_now=True)
845        except SyntaxError:
846            pass
847        else:
848            msg = 'Bad file did not raise error!'
849            raise Exception, msg
850        os.remove(fileName)
851
852    def test_load_csv_badIV(self):
853        ''' test_load_csv_bad(self):
854        header column, body column missmatch
855        (Format error)
856        '''
857
858        fileName = tempfile.mktemp('.txt')
859        file = open(fileName, 'w')
860        file.write(' x,y,elevation ,  speed \n\
8611.0, 0.0, 10.0, wow\n\
8620.0, 1.0, 0.0, ha\n\
8631.0, 0.0 ,10.4, yeah\n')
864        file.close()
865
866        results = Geospatial_data(fileName, max_read_lines=2,
867                                  load_file_now=False)
868
869        # Blocking
870        geo_list = []
871        try:
872            for i in results:
873                geo_list.append(i)
874        except SyntaxError:
875            pass
876        else:
877            msg = 'Bad file did not raise error!'
878            raise Exception, msg
879        os.remove(fileName)
880
881    def test_load_pts_blocking(self):
882        #This is pts!
883        fileName = tempfile.mktemp('.txt')
884        file = open(fileName, 'w')
885        file.write(' x,y, elevation ,  speed \n\
8861.0, 0.0, 10.0, 0.0\n\
8870.0, 1.0, 0.0, 10.0\n\
8881.0, 0.0 ,10.4, 40.0\n')
889        file.close()
890
891        pts_file = tempfile.mktemp('.pts')
892
893        convert = Geospatial_data(fileName)
894        convert.export_points_file(pts_file)
895        results = Geospatial_data(pts_file, max_read_lines=2)
896
897        assert num.allclose(results.get_data_points(),
898                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
899        assert num.allclose(results.get_attributes(attribute_name='elevation'),
900                            [10.0, 0.0, 10.4])
901        assert num.allclose(results.get_attributes(attribute_name='speed'),
902                            [0.0, 10.0, 40.0])
903
904        # Blocking
905        geo_list = []
906        for i in results:
907            geo_list.append(i)
908        assert num.allclose(geo_list[0].get_data_points(),
909                            [[1.0, 0.0],[0.0, 1.0]])
910        assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
911                            [10.0, 0.0])
912        assert num.allclose(geo_list[1].get_data_points(),
913                            [[1.0, 0.0]])
914        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
915                            [10.4])
916
917        os.remove(fileName)
918        os.remove(pts_file)
919
920    def verbose_test_load_pts_blocking(self):
921        fileName = tempfile.mktemp('.txt')
922        file = open(fileName, 'w')
923        file.write(' x,y, elevation ,  speed \n\
9241.0, 0.0, 10.0, 0.0\n\
9250.0, 1.0, 0.0, 10.0\n\
9261.0, 0.0, 10.0, 0.0\n\
9270.0, 1.0, 0.0, 10.0\n\
9281.0, 0.0, 10.0, 0.0\n\
9290.0, 1.0, 0.0, 10.0\n\
9301.0, 0.0, 10.0, 0.0\n\
9310.0, 1.0, 0.0, 10.0\n\
9321.0, 0.0, 10.0, 0.0\n\
9330.0, 1.0, 0.0, 10.0\n\
9341.0, 0.0, 10.0, 0.0\n\
9350.0, 1.0, 0.0, 10.0\n\
9361.0, 0.0, 10.0, 0.0\n\
9370.0, 1.0, 0.0, 10.0\n\
9381.0, 0.0, 10.0, 0.0\n\
9390.0, 1.0, 0.0, 10.0\n\
9401.0, 0.0, 10.0, 0.0\n\
9410.0, 1.0, 0.0, 10.0\n\
9421.0, 0.0, 10.0, 0.0\n\
9430.0, 1.0, 0.0, 10.0\n\
9441.0, 0.0, 10.0, 0.0\n\
9450.0, 1.0, 0.0, 10.0\n\
9461.0, 0.0, 10.0, 0.0\n\
9470.0, 1.0, 0.0, 10.0\n\
9481.0, 0.0, 10.0, 0.0\n\
9490.0, 1.0, 0.0, 10.0\n\
9501.0, 0.0, 10.0, 0.0\n\
9510.0, 1.0, 0.0, 10.0\n\
9521.0, 0.0, 10.0, 0.0\n\
9530.0, 1.0, 0.0, 10.0\n\
9541.0, 0.0, 10.0, 0.0\n\
9550.0, 1.0, 0.0, 10.0\n\
9561.0, 0.0, 10.0, 0.0\n\
9570.0, 1.0, 0.0, 10.0\n\
9581.0, 0.0, 10.0, 0.0\n\
9590.0, 1.0, 0.0, 10.0\n\
9601.0, 0.0, 10.0, 0.0\n\
9610.0, 1.0, 0.0, 10.0\n\
9621.0, 0.0, 10.0, 0.0\n\
9630.0, 1.0, 0.0, 10.0\n\
9641.0, 0.0, 10.0, 0.0\n\
9650.0, 1.0, 0.0, 10.0\n\
9661.0, 0.0, 10.0, 0.0\n\
9670.0, 1.0, 0.0, 10.0\n\
9681.0, 0.0, 10.0, 0.0\n\
9690.0, 1.0, 0.0, 10.0\n\
9701.0, 0.0 ,10.4, 40.0\n')
971        file.close()
972
973        pts_file = tempfile.mktemp('.pts')
974        convert = Geospatial_data(fileName)
975        convert.export_points_file(pts_file)
976        results = Geospatial_data(pts_file, max_read_lines=2, verbose=True)
977
978        # Blocking
979        geo_list = []
980        for i in results:
981            geo_list.append(i)
982        assert num.allclose(geo_list[0].get_data_points(),
983                            [[1.0, 0.0], [0.0, 1.0]])
984        assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
985                            [10.0, 0.0])
986        assert num.allclose(geo_list[1].get_data_points(),
987                            [[1.0, 0.0],[0.0, 1.0] ])
988        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
989                            [10.0, 0.0])
990
991        os.remove(fileName)
992        os.remove(pts_file)
993
994    def test_new_export_pts_file(self):
995        att_dict = {}
996        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
997        att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
998        att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
999
1000        fileName = tempfile.mktemp('.pts')
1001        G = Geospatial_data(pointlist, att_dict)
1002        G.export_points_file(fileName)
1003        results = Geospatial_data(file_name = fileName)
1004        os.remove(fileName)
1005
1006        assert num.allclose(results.get_data_points(),
1007                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1008        assert num.allclose(results.get_attributes(attribute_name='elevation'),
1009                            [10.1, 0.0, 10.4])
1010        answer = [10.0, 1.0, 10.4]
1011        assert num.allclose(results.get_attributes(attribute_name='brightness'),
1012                            answer)
1013
1014    def test_new_export_absolute_pts_file(self):
1015        att_dict = {}
1016        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1017        att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
1018        att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
1019        geo_ref = Geo_reference(50, 25, 55)
1020
1021        fileName = tempfile.mktemp('.pts')
1022        G = Geospatial_data(pointlist, att_dict, geo_ref)
1023        G.export_points_file(fileName, absolute=True)
1024        results = Geospatial_data(file_name = fileName)
1025        os.remove(fileName)
1026
1027        msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s'
1028               % (str(results.get_data_points()), str(G.get_data_points(True))))
1029        assert num.allclose(results.get_data_points(),
1030                            G.get_data_points(True)), msg
1031        msg = ("results.get_attributes(attribute_name='elevation')=%s"
1032               % str(results.get_attributes(attribute_name='elevation')))
1033        assert num.allclose(results.get_attributes(attribute_name='elevation'),
1034                            [10.1, 0.0, 10.4]), msg
1035        answer = [10.0, 1.0, 10.4]
1036        msg = ("results.get_attributes(attribute_name='brightness')=%s, "
1037               'answer=%s' %
1038               (str(results.get_attributes(attribute_name='brightness')),
1039                str(answer)))
1040        assert num.allclose(results.get_attributes(attribute_name='brightness'),
1041                            answer), msg
1042
1043    def test_loadpts(self):
1044        from Scientific.IO.NetCDF import NetCDFFile
1045
1046        fileName = tempfile.mktemp('.pts')
1047        # NetCDF file definition
1048        outfile = NetCDFFile(fileName, netcdf_mode_w)
1049
1050        # dimension definitions
1051        outfile.createDimension('number_of_points', 3)
1052        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1053
1054        # variable definitions
1055        outfile.createVariable('points', netcdf_float, ('number_of_points',
1056                                                        'number_of_dimensions'))
1057        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1058
1059        # Get handles to the variables
1060        points = outfile.variables['points']
1061        elevation = outfile.variables['elevation']
1062
1063        points[0, :] = [1.0,0.0]
1064        elevation[0] = 10.0
1065        points[1, :] = [0.0,1.0]
1066        elevation[1] = 0.0
1067        points[2, :] = [1.0,0.0]
1068        elevation[2] = 10.4
1069
1070        outfile.close()
1071
1072        results = Geospatial_data(file_name = fileName)
1073        os.remove(fileName)
1074
1075        answer = [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
1076        assert num.allclose(results.get_data_points(),
1077                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1078        assert num.allclose(results.get_attributes(attribute_name='elevation'),
1079                            [10.0, 0.0, 10.4])
1080
1081    def test_writepts(self):
1082        '''Test that storage of x,y,attributes works'''
1083
1084        att_dict = {}
1085        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1086        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1087        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
1088        geo_reference=Geo_reference(56, 1.9, 1.9)
1089
1090        # Test pts format
1091        fileName = tempfile.mktemp('.pts')
1092        G = Geospatial_data(pointlist, att_dict, geo_reference)
1093        G.export_points_file(fileName, False)
1094        results = Geospatial_data(file_name=fileName)
1095        os.remove(fileName)
1096
1097        assert num.allclose(results.get_data_points(False),
1098                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1099        assert num.allclose(results.get_attributes('elevation'),
1100                            [10.0, 0.0, 10.4])
1101        answer = [10.0, 0.0, 10.4]
1102        assert num.allclose(results.get_attributes('brightness'), answer)
1103        self.failUnless(geo_reference == geo_reference,
1104                        'test_writepts failed. Test geo_reference')
1105
1106    def test_write_csv_attributes(self):
1107        '''Test that storage of x,y,attributes works'''
1108
1109        att_dict = {}
1110        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1111        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1112        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
1113        geo_reference=Geo_reference(56, 0, 0)
1114
1115        # Test txt format
1116        fileName = tempfile.mktemp('.txt')
1117        G = Geospatial_data(pointlist, att_dict, geo_reference)
1118        G.export_points_file(fileName)
1119        results = Geospatial_data(file_name=fileName)
1120        os.remove(fileName)
1121
1122        assert num.allclose(results.get_data_points(False),
1123                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1124        assert num.allclose(results.get_attributes('elevation'),
1125                            [10.0, 0.0, 10.4])
1126        answer = [10.0, 0.0, 10.4]
1127        assert num.allclose(results.get_attributes('brightness'), answer)
1128
1129    def test_write_csv_attributes_lat_long(self):
1130        '''Test that storage of x,y,attributes works'''
1131
1132        att_dict = {}
1133        pointlist = num.array([[-21.5, 114.5], [-21.6, 114.5], [-21.7, 114.5]])
1134        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1135        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
1136
1137        # Test txt format
1138        fileName = tempfile.mktemp('.txt')
1139        G = Geospatial_data(pointlist, att_dict, points_are_lats_longs=True)
1140        G.export_points_file(fileName, as_lat_long=True)
1141        results = Geospatial_data(file_name=fileName)
1142        os.remove(fileName)
1143
1144        assert num.allclose(results.get_data_points(False, as_lat_long=True),
1145                            pointlist)
1146        assert num.allclose(results.get_attributes('elevation'),
1147                            [10.0, 0.0, 10.4])
1148        answer = [10.0, 0.0, 10.4]
1149        assert num.allclose(results.get_attributes('brightness'), answer)
1150
1151    def test_writepts_no_attributes(self):
1152        '''Test that storage of x,y alone works'''
1153
1154        att_dict = {}
1155        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1156        geo_reference=Geo_reference(56, 1.9, 1.9)
1157
1158        # Test pts format
1159        fileName = tempfile.mktemp('.pts')
1160        G = Geospatial_data(pointlist, None, geo_reference)
1161        G.export_points_file(fileName, False)
1162        results = Geospatial_data(file_name=fileName)
1163        os.remove(fileName)
1164
1165        assert num.allclose(results.get_data_points(False),
1166                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1167        self.failUnless(geo_reference == geo_reference,
1168                        'test_writepts failed. Test geo_reference')
1169
1170    def test_write_csv_no_attributes(self):
1171        '''Test that storage of x,y alone works'''
1172
1173        att_dict = {}
1174        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1175        geo_reference=Geo_reference(56,0,0)
1176
1177        # Test format
1178        fileName = tempfile.mktemp('.txt')
1179        G = Geospatial_data(pointlist, None, geo_reference)
1180        G.export_points_file(fileName)
1181        results = Geospatial_data(file_name=fileName)
1182        os.remove(fileName)
1183
1184        assert num.allclose(results.get_data_points(False),
1185                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1186
1187################################################################################
1188# Check bad PTS files.
1189################################################################################
1190
1191    def test_load_bad_no_file_pts(self):
1192        fileName = tempfile.mktemp('.pts')
1193        try:
1194            results = Geospatial_data(file_name=fileName)
1195        except IOError:
1196            pass
1197        else:
1198            msg = 'imaginary file did not raise error!'
1199            raise Exception, msg
1200
1201    def test_create_from_pts_file(self):
1202        from Scientific.IO.NetCDF import NetCDFFile
1203
1204        # NetCDF file definition
1205        FN = 'test_points.pts'
1206        outfile = NetCDFFile(FN, netcdf_mode_w)
1207
1208        # dimension definitions
1209        outfile.createDimension('number_of_points', 3)
1210        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1211
1212        # variable definitions
1213        outfile.createVariable('points', netcdf_float, ('number_of_points',
1214                                                        'number_of_dimensions'))
1215        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1216
1217        # Get handles to the variables
1218        points = outfile.variables['points']
1219        elevation = outfile.variables['elevation']
1220
1221        points[0, :] = [1.0,0.0]
1222        elevation[0] = 10.0
1223        points[1, :] = [0.0,1.0]
1224        elevation[1] = 0.0
1225        points[2, :] = [1.0,0.0]
1226        elevation[2] = 10.4
1227
1228        outfile.close()
1229
1230        G = Geospatial_data(file_name = FN)
1231
1232        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
1233        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
1234        assert num.allclose(G.get_data_points(),
1235                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1236        assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1237        os.remove(FN)
1238
1239    def test_create_from_pts_file_with_geo(self):
1240        '''Test if Geospatial data is correctly instantiated from a pts file.'''
1241
1242        from Scientific.IO.NetCDF import NetCDFFile
1243
1244        # NetCDF file definition
1245        FN = 'test_points.pts'
1246        outfile = NetCDFFile(FN, netcdf_mode_w)
1247
1248        # Make up an arbitrary georef
1249        xll = 0.1
1250        yll = 20
1251        geo_reference=Geo_reference(56, xll, yll)
1252        geo_reference.write_NetCDF(outfile)
1253
1254        # dimension definitions
1255        outfile.createDimension('number_of_points', 3)
1256        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1257
1258        # variable definitions
1259        outfile.createVariable('points', netcdf_float, ('number_of_points',
1260                                                        'number_of_dimensions'))
1261        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1262
1263        # Get handles to the variables
1264        points = outfile.variables['points']
1265        elevation = outfile.variables['elevation']
1266
1267        points[0, :] = [1.0,0.0]
1268        elevation[0] = 10.0
1269        points[1, :] = [0.0,1.0]
1270        elevation[1] = 0.0
1271        points[2, :] = [1.0,0.0]
1272        elevation[2] = 10.4
1273
1274        outfile.close()
1275
1276        G = Geospatial_data(file_name = FN)
1277
1278        assert num.allclose(G.get_geo_reference().get_xllcorner(), xll)
1279        assert num.allclose(G.get_geo_reference().get_yllcorner(), yll)
1280        assert num.allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
1281                                                  [0.0+xll, 1.0+yll],
1282                                                  [1.0+xll, 0.0+yll]])
1283        assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4])
1284
1285        os.remove(FN)
1286
1287    def test_add_(self):
1288        '''test_add_(self):
1289        adds an txt and pts files, reads the files and adds them
1290        checking results are correct
1291        '''
1292
1293        # create files
1294        att_dict1 = {}
1295        pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1296        att_dict1['elevation'] = num.array([-10.0, 0.0, 10.4])
1297        att_dict1['brightness'] = num.array([10.0, 0.0, 10.4])
1298        geo_reference1 = Geo_reference(56, 2.0, 1.0)
1299
1300        att_dict2 = {}
1301        pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0], [2.0, 1.0]])
1302        att_dict2['elevation'] = num.array([1.0, 15.0, 1.4])
1303        att_dict2['brightness'] = num.array([14.0, 1.0, -12.4])
1304        geo_reference2 = Geo_reference(56, 1.0, 2.0)
1305
1306        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1307        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1308
1309        fileName1 = tempfile.mktemp('.txt')
1310        fileName2 = tempfile.mktemp('.pts')
1311
1312        # makes files
1313        G1.export_points_file(fileName1)
1314        G2.export_points_file(fileName2)
1315
1316        # add files
1317        G3 = Geospatial_data(file_name=fileName1)
1318        G4 = Geospatial_data(file_name=fileName2)
1319        G = G3 + G4
1320
1321        #read results
1322        assert num.allclose(G.get_data_points(),
1323                            [[ 3.0, 1.0], [ 2.0, 2.0],
1324                             [ 3.0, 1.0], [ 3.0, 3.0],
1325                             [ 2.0, 4.0], [ 3.0, 3.0]])
1326        assert num.allclose(G.get_attributes(attribute_name='elevation'),
1327                            [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
1328        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
1329        assert num.allclose(G.get_attributes(attribute_name='brightness'),
1330                            answer)
1331        self.failUnless(G.get_geo_reference() == geo_reference1,
1332                        'test_writepts failed. Test geo_reference')
1333
1334        os.remove(fileName1)
1335        os.remove(fileName2)
1336
1337    def test_ensure_absolute(self):
1338        points = [[2.0, 0.0], [1.0, 1.0],
1339                  [2.0, 0.0], [2.0, 2.0],
1340                  [1.0, 3.0], [2.0, 2.0]]
1341        new_points = ensure_absolute(points)
1342
1343        assert num.allclose(new_points, points)
1344
1345        points = num.array([[2.0, 0.0], [1.0, 1.0],
1346                            [2.0, 0.0], [2.0, 2.0],
1347                            [1.0, 3.0], [2.0, 2.0]])
1348        new_points = ensure_absolute(points)
1349
1350        assert num.allclose(new_points, points)
1351
1352        ab_points = num.array([[2.0, 0.0], [1.0, 1.0],
1353                               [2.0, 0.0], [2.0, 2.0],
1354                               [1.0, 3.0], [2.0, 2.0]])
1355
1356        mesh_origin = (56, 290000, 618000) # zone, easting, northing
1357        data_points = num.zeros((ab_points.shape), num.float)
1358
1359        #Shift datapoints according to new origins
1360        for k in range(len(ab_points)):
1361            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1362            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1363        new_points = ensure_absolute(data_points, geo_reference=mesh_origin)
1364
1365        assert num.allclose(new_points, ab_points)
1366
1367        geo = Geo_reference(56,67,-56)
1368        data_points = geo.change_points_geo_ref(ab_points)
1369        new_points = ensure_absolute(data_points, geo_reference=geo)
1370
1371        assert num.allclose(new_points, ab_points)
1372
1373        geo_reference = Geo_reference(56, 100, 200)
1374        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1375        points = geo_reference.change_points_geo_ref(ab_points)
1376        attributes = [2, 4]
1377        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
1378        new_points = ensure_absolute(G)
1379
1380        assert num.allclose(new_points, ab_points)
1381
1382    def test_ensure_geospatial(self):
1383        points = [[2.0, 0.0], [1.0, 1.0],
1384                  [2.0, 0.0], [2.0, 2.0],
1385                  [1.0, 3.0], [2.0, 2.0]]
1386        new_points = ensure_geospatial(points)
1387        assert num.allclose(new_points.get_data_points(absolute=True), points)
1388
1389        points = num.array([[2.0, 0.0], [1.0, 1.0],
1390                            [2.0, 0.0], [2.0, 2.0],
1391                            [1.0, 3.0], [2.0, 2.0]])
1392        new_points = ensure_geospatial(points)
1393        assert num.allclose(new_points.get_data_points(absolute=True), points)
1394
1395        ab_points = num.array([[2.0, 0.0],[1.0, 1.0],
1396                               [2.0, 0.0],[2.0, 2.0],
1397                               [1.0, 3.0],[2.0, 2.0]])
1398        mesh_origin = (56, 290000, 618000)      # zone, easting, northing
1399        data_points = num.zeros((ab_points.shape), num.float)
1400
1401        #Shift datapoints according to new origins
1402        for k in range(len(ab_points)):
1403            data_points[k][0] = ab_points[k][0] - mesh_origin[1]
1404            data_points[k][1] = ab_points[k][1] - mesh_origin[2]
1405        new_geospatial = ensure_geospatial(data_points,
1406                                           geo_reference=mesh_origin)
1407        new_points = new_geospatial.get_data_points(absolute=True)
1408        assert num.allclose(new_points, ab_points)
1409
1410        geo = Geo_reference(56, 67, -56)
1411        data_points = geo.change_points_geo_ref(ab_points)
1412        new_geospatial = ensure_geospatial(data_points, geo_reference=geo)
1413        new_points = new_geospatial.get_data_points(absolute=True)
1414        assert num.allclose(new_points, ab_points)
1415
1416        geo_reference = Geo_reference(56, 100, 200)
1417        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1418        points = geo_reference.change_points_geo_ref(ab_points)
1419        attributes = [2, 4]
1420        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
1421        new_geospatial = ensure_geospatial(G)
1422        new_points = new_geospatial.get_data_points(absolute=True)
1423        assert num.allclose(new_points, ab_points)
1424
1425    def test_isinstance(self):
1426        fileName = tempfile.mktemp('.csv')
1427        file = open(fileName, 'w')
1428        file.write('x,y,  elevation ,  speed \n\
14291.0, 0.0, 10.0, 0.0\n\
14300.0, 1.0, 0.0, 10.0\n\
14311.0, 0.0, 10.4, 40.0\n')
1432        file.close()
1433
1434        results = Geospatial_data(fileName)
1435        assert num.allclose(results.get_data_points(absolute=True),
1436                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1437        assert num.allclose(results.get_attributes(attribute_name='elevation'),
1438                            [10.0, 0.0, 10.4])
1439        assert num.allclose(results.get_attributes(attribute_name='speed'),
1440                            [0.0, 10.0, 40.0])
1441
1442        os.remove(fileName)
1443
1444    def test_no_constructors(self):
1445        try:
1446            G = Geospatial_data()
1447        except ValueError:
1448            pass
1449        else:
1450            msg = 'Instance must have a filename or data points'
1451            raise Exception, msg
1452
1453    def test_load_csv_lat_long(self):
1454        '''comma delimited'''
1455
1456        fileName = tempfile.mktemp('.csv')
1457        file = open(fileName, 'w')
1458        file.write('long,lat, elevation, yeah \n\
1459150.916666667,-34.50,452.688000, 10\n\
1460150.0,-34,459.126000, 10\n')
1461        file.close()
1462
1463        results = Geospatial_data(fileName)
1464        os.remove(fileName)
1465        points = results.get_data_points()
1466
1467        assert num.allclose(points[0][0], 308728.009)
1468        assert num.allclose(points[0][1], 6180432.601)
1469        assert num.allclose(points[1][0],  222908.705)
1470        assert num.allclose(points[1][1], 6233785.284)
1471
1472    def test_load_csv_lat_longII(self):
1473        '''comma delimited'''
1474
1475        fileName = tempfile.mktemp('.csv')
1476        file = open(fileName, 'w')
1477        file.write('Lati,LONG,z \n\
1478-34.50,150.916666667,452.688000\n\
1479-34,150.0,459.126000\n')
1480        file.close()
1481
1482        results = Geospatial_data(fileName)
1483        os.remove(fileName)
1484        points = results.get_data_points()
1485
1486        assert num.allclose(points[0][0], 308728.009)
1487        assert num.allclose(points[0][1], 6180432.601)
1488        assert num.allclose(points[1][0], 222908.705)
1489        assert num.allclose(points[1][1], 6233785.284)
1490
1491    def test_load_csv_lat_long_bad(self):
1492        '''comma delimited'''
1493
1494        fileName = tempfile.mktemp('.csv')
1495        file = open(fileName, 'w')
1496        file.write('Lati,LONG,z \n\
1497-25.0,180.0,452.688000\n\
1498-34,150.0,459.126000\n')
1499        file.close()
1500
1501        try:
1502            results = Geospatial_data(fileName)
1503        except ANUGAError:
1504            pass
1505        else:
1506            msg = 'Different zones in Geo references not caught.'
1507            raise Exception, msg
1508
1509        os.remove(fileName)
1510
1511    def test_lat_long(self):
1512        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1513        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1514
1515        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1516        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1517
1518        lats = [lat_gong, lat_2]
1519        longs = [lon_gong, lon_2]
1520        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1521
1522        points = gsd.get_data_points(absolute=True)
1523
1524        assert num.allclose(points[0][0], 308728.009)
1525        assert num.allclose(points[0][1], 6180432.601)
1526        assert num.allclose(points[1][0], 222908.705)
1527        assert num.allclose(points[1][1], 6233785.284)
1528        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1529                        'Bad zone error!')
1530
1531        try:
1532            results = Geospatial_data(latitudes=lats)
1533        except ValueError:
1534            pass
1535        else:
1536            self.fail('Error not thrown error!')
1537        try:
1538            results = Geospatial_data(latitudes=lats)
1539        except ValueError:
1540            pass
1541        else:
1542            self.fail('Error not thrown error!')
1543        try:
1544            results = Geospatial_data(longitudes=lats)
1545        except ValueError:
1546            pass
1547        else:
1548            self.fail('Error not thrown error!')
1549        try:
1550            results = Geospatial_data(latitudes=lats, longitudes=longs,
1551                                      geo_reference="p")
1552        except ValueError:
1553            pass
1554        else:
1555            self.fail('Error not thrown error!')
1556
1557        try:
1558            results = Geospatial_data(latitudes=lats, longitudes=longs,
1559                                      data_points=12)
1560        except ValueError:
1561            pass
1562        else:
1563            self.fail('Error not thrown error!')
1564
1565    def test_lat_long2(self):
1566        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1567        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1568
1569        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1570        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1571
1572        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1573        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1574
1575        points = gsd.get_data_points(absolute=True)
1576
1577        assert num.allclose(points[0][0], 308728.009)
1578        assert num.allclose(points[0][1], 6180432.601)
1579        assert num.allclose(points[1][0], 222908.705)
1580        assert num.allclose(points[1][1], 6233785.284)
1581        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1582                        'Bad zone error!')
1583
1584        try:
1585            results = Geospatial_data(points_are_lats_longs=True)
1586        except ValueError:
1587            pass
1588        else:
1589            self.fail('Error not thrown error!')
1590
1591    def test_write_urs_file(self):
1592        lat_gong = degminsec2decimal_degrees(-34, 00, 0)
1593        lon_gong = degminsec2decimal_degrees(150, 30, 0.)
1594
1595        lat_2 = degminsec2decimal_degrees(-34, 00, 1)
1596        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1597        p1 = (lat_gong, lon_gong)
1598        p2 = (lat_2, lon_2)
1599        points = ImmutableSet([p1, p2, p1])
1600        gsd = Geospatial_data(data_points=list(points),
1601                              points_are_lats_longs=True)
1602
1603        fn = tempfile.mktemp('.urs')
1604        gsd.export_points_file(fn)
1605        handle = open(fn)
1606        lines = handle.readlines()
1607        assert lines[0], '2'
1608        assert lines[1], '-34.0002778 150.0 0'
1609        assert lines[2], '-34.0 150.5 1'
1610        handle.close()
1611        os.remove(fn)
1612
1613    def test_lat_long_set(self):
1614        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1615        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1616
1617        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1618        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1619        p1 = (lat_gong, lon_gong)
1620        p2 = (lat_2, lon_2)
1621        points = ImmutableSet([p1, p2, p1])
1622        gsd = Geospatial_data(data_points=list(points),
1623                              points_are_lats_longs=True)
1624
1625        points = gsd.get_data_points(absolute=True)
1626        #Note the order is unknown, due to using sets
1627        # and it changes from windows to linux
1628        try:
1629            assert num.allclose(points[1][0], 308728.009)
1630            assert num.allclose(points[1][1], 6180432.601)
1631            assert num.allclose(points[0][0], 222908.705)
1632            assert num.allclose(points[0][1], 6233785.284)
1633        except AssertionError:
1634            assert num.allclose(points[0][0], 308728.009)
1635            assert num.allclose(points[0][1], 6180432.601)
1636            assert num.allclose(points[1][0], 222908.705)
1637            assert num.allclose(points[1][1], 6233785.284)
1638
1639        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1640                        'Bad zone error!')
1641        points = gsd.get_data_points(as_lat_long=True)
1642        try:
1643            assert num.allclose(points[0][0], -34)
1644            assert num.allclose(points[0][1], 150)
1645        except AssertionError:
1646            assert num.allclose(points[1][0], -34)
1647            assert num.allclose(points[1][1], 150)
1648
1649    def test_len(self):
1650        points = [[1.0, 2.1], [3.0, 5.3]]
1651        G = Geospatial_data(points)
1652        self.failUnless(2 == len(G), 'Len error!')
1653
1654        points = [[1.0, 2.1]]
1655        G = Geospatial_data(points)
1656        self.failUnless(1 == len(G), 'Len error!')
1657
1658        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1659        G = Geospatial_data(points)
1660        self.failUnless(4 == len(G), 'Len error!')
1661
1662    def test_split(self):
1663        """test if the results from spilt are disjoin sets"""
1664
1665        # below is a workaround until randint works on cyclones compute nodes
1666        if get_host_name()[8:9] != '0':
1667            points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
1668                      [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
1669                      [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0],
1670                      [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
1671                      [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
1672            attributes = {'depth': [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                          'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
1676                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
1677                                    21, 22, 23, 24, 25]}
1678            G = Geospatial_data(points, attributes)
1679
1680            factor = 0.21
1681
1682            # will return G1 with 10% of points and G2 with 90%
1683            G1, G2  = G.split(factor, 100)
1684            assert num.allclose(len(G), len(G1)+len(G2))
1685            assert num.allclose(round(len(G)*factor), len(G1))
1686
1687            P = G1.get_data_points(absolute=False)
1688            expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0],
1689                        [3.0, 1.0], [2.0, 3.0]]
1690            msg = 'Expected %s, but\nP=%s' % (str(expected), str(P))
1691            assert num.allclose(P, expected), msg
1692
1693            A = G1.get_attributes()
1694            expected = [24, 18, 17, 11, 8]
1695            msg = 'expected=%s, but A=%s' % (str(expected), str(A))
1696            assert num.allclose(A, expected), msg
1697
1698    def test_split1(self):
1699        """test if the results from split are disjoint sets"""
1700
1701        # below is a workaround until randint works on cyclones compute nodes
1702        if get_host_name()[8:9] != '0':
1703            from RandomArray import randint,seed
1704
1705            seed(100,100)
1706            a_points = randint(0, 999999, (10,2))
1707            points = a_points.tolist()
1708
1709            G = Geospatial_data(points)
1710
1711            factor = 0.1
1712
1713            # will return G1 with 10% of points and G2 with 90%
1714            G1, G2  = G.split(factor, 100)
1715            assert num.allclose(len(G), len(G1)+len(G2))
1716            assert num.allclose(round(len(G)*factor), len(G1))
1717
1718            P = G1.get_data_points(absolute=False)
1719            expected = [[982420., 28233.]]
1720            msg = 'expected=%s, but\nP=%s' % (str(expected), str(P))
1721            assert num.allclose(P, expected), msg
1722
1723    def test_find_optimal_smoothing_parameter(self):
1724        """
1725        Creates a elevation file representing a hill (sort of) and runs
1726        find_optimal_smoothing_parameter for 3 different alphas,
1727
1728        NOTE the random number seed is provided to control the results
1729        """
1730
1731        from cmath import cos
1732
1733        # below is a workaround until randint works on cyclones compute nodes
1734        if get_host_name()[8:9]!='0':
1735            filename = tempfile.mktemp('.csv')
1736            file = open(filename, 'w')
1737            file.write('x,y,elevation \n')
1738
1739            for i in range(-5, 6):
1740                for j in range(-5, 6):
1741                    # this equation makes a surface like a circle ripple
1742                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1743                    file.write("%s, %s, %s\n" % (i, j, z))
1744
1745            file.close()
1746
1747            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1748                                                 alpha_list=[0.0001, 0.01, 1],
1749                                                 mesh_file=None,
1750                                                 mesh_resolution=3,
1751                                                 north_boundary=5,
1752                                                 south_boundary=-5,
1753                                                 east_boundary=5,
1754                                                 west_boundary=-5,
1755                                                 plot_name=None,
1756                                                 seed_num=100000,
1757                                                 verbose=False)
1758
1759            os.remove(filename)
1760
1761            # print value, alpha
1762            assert(alpha == 0.01)
1763
1764    def test_find_optimal_smoothing_parameter1(self):
1765        """
1766        Creates a elevation file representing  a hill (sort of) and
1767        Then creates a mesh file and passes the mesh file and the elevation
1768        file to find_optimal_smoothing_parameter for 3 different alphas,
1769
1770        NOTE the random number seed is provided to control the results
1771        """
1772
1773        # below is a workaround until randint works on cyclones compute nodes
1774        if get_host_name()[8:9]!='0':
1775            from cmath import cos
1776            from anuga.pmesh.mesh_interface import create_mesh_from_regions
1777
1778            filename = tempfile.mktemp('.csv')
1779            file = open(filename, 'w')
1780            file.write('x,y,elevation \n')
1781
1782            for i in range(-5 ,6):
1783                for j in range(-5, 6):
1784                    # this equation makes a surface like a circle ripple
1785                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1786                    file.write('%s, %s, %s\n' % (i, j, z))
1787
1788            file.close()
1789
1790            poly=[[5,5], [5,-5], [-5,-5], [-5,5]]
1791            internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]]
1792            mesh_filename= tempfile.mktemp('.msh')
1793
1794            create_mesh_from_regions(poly,
1795                                     boundary_tags={'back': [2],
1796                                                    'side': [1,3],
1797                                                    'ocean': [0]},
1798                                     maximum_triangle_area=3,
1799                                     interior_regions=internal_poly,
1800                                     filename=mesh_filename,
1801                                     use_cache=False,
1802                                     verbose=False)
1803
1804            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1805                                                 alpha_list=[0.0001, 0.01, 1],
1806                                                 mesh_file=mesh_filename,
1807                                                 plot_name=None,
1808                                                 seed_num=174,
1809                                                 verbose=False)
1810
1811            os.remove(filename)
1812            os.remove(mesh_filename)
1813
1814            assert(alpha == 0.01)
1815
1816    def test_find_optimal_smoothing_parameter2(self):
1817        '''Tests requirement that mesh file must exist or IOError is thrown
1818
1819        NOTE the random number seed is provided to control the results
1820        '''
1821
1822        from cmath import cos
1823        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1824
1825        filename = tempfile.mktemp('.csv')
1826        mesh_filename= tempfile.mktemp('.msh')
1827
1828        try:
1829            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1830                                             alpha_list=[0.0001, 0.01, 1],
1831                                             mesh_file=mesh_filename,
1832                                             plot_name=None,
1833                                             seed_num=174,
1834                                             verbose=False)
1835        except IOError:
1836            pass
1837        else:
1838            self.fail('Error not thrown error!')
1839
1840################################################################################
1841
1842if __name__ == "__main__":
1843    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1844    runner = unittest.TextTestRunner() #verbosity=2)
1845    runner.run(suite)
1846
1847
Note: See TracBrowser for help on using the repository browser.