source: trunk/anuga_core/source/anuga/geospatial_data/test_geospatial_data.py @ 7779

Last change on this file since 7779 was 7779, checked in by hudson, 14 years ago

Fixed wrong exception module import paths.

File size: 68.1 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4import os
5import tempfile
6from math import sqrt, pi
7
8import numpy as num
9
10from anuga.geospatial_data.geospatial_data import *
11from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
12from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
13from anuga.anuga_exceptions import ANUGAError
14from anuga.utilities.system_tools import get_host_name
15from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
16from anuga.config import netcdf_float
17
18
19class Test_Geospatial_data(unittest.TestCase):
20    def setUp(self):
21        pass
22
23    def tearDown(self):
24        pass
25
26    def test_0(self):
27        #Basic points
28        from anuga.coordinate_transforms.geo_reference import Geo_reference
29
30        points = [[1.0, 2.1], [3.0, 5.3]]
31        G = Geospatial_data(points)
32        assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
33
34        # Check __repr__
35        # FIXME (Ole): Is this really machine independent?
36        rep = `G`
37        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
38        msg = 'Representation %s is not equal to %s' % (rep, ref)
39        assert rep == ref, msg
40
41        #Check getter
42        assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
43
44        #Check defaults
45        assert G.attributes is None
46        assert G.geo_reference.zone == Geo_reference().zone
47        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
48        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
49
50    def test_1(self):
51        points = [[1.0, 2.1], [3.0, 5.3]]
52        attributes = [2, 4]
53        G = Geospatial_data(points, attributes)
54        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
55        assert num.allclose(G.attributes.values()[0], [2, 4])
56
57    def test_2(self):
58        from anuga.coordinate_transforms.geo_reference import Geo_reference
59
60        points = [[1.0, 2.1], [3.0, 5.3]]
61        attributes = [2, 4]
62        G = Geospatial_data(points, attributes,
63                            geo_reference=Geo_reference(56, 100, 200))
64
65        assert G.geo_reference.zone == 56
66        assert G.geo_reference.xllcorner == 100
67        assert G.geo_reference.yllcorner == 200
68
69    def test_get_attributes_1(self):
70        from anuga.coordinate_transforms.geo_reference import Geo_reference
71
72        points = [[1.0, 2.1], [3.0, 5.3]]
73        attributes = [2, 4]
74        G = Geospatial_data(points, attributes,
75                            geo_reference=Geo_reference(56, 100, 200))
76
77        P = G.get_data_points(absolute=False)
78        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
79
80        P = G.get_data_points(absolute=True)
81        assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])
82
83        V = G.get_attributes() #Simply get them
84        assert num.allclose(V, [2, 4])
85
86        V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name
87        assert num.allclose(V, [2, 4])
88
89    def test_get_attributes_2(self):
90        #Multiple attributes
91        from anuga.coordinate_transforms.geo_reference import Geo_reference
92
93        points = [[1.0, 2.1], [3.0, 5.3]]
94        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
95        G = Geospatial_data(points, attributes,
96                            geo_reference=Geo_reference(56, 100, 200),
97                            default_attribute_name='a1')
98
99        P = G.get_data_points(absolute=False)
100        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
101
102        V = G.get_attributes() #Get default attribute
103        assert num.allclose(V, [2, 4])
104
105        V = G.get_attributes('a0') #Get by name
106        assert num.allclose(V, [0, 0])
107
108        V = G.get_attributes('a1') #Get by name
109        assert num.allclose(V, [2, 4])
110
111        V = G.get_attributes('a2') #Get by name
112        assert num.allclose(V, [79.4, -7])
113
114                # FIXME: use failUnlessRaises()
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), axis=0))    #??default#
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), axis=0)) #??default#
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        msg = ('new_points=\n%s\nab_points=\n%s'
1415               % (str(new_points), str(ab_points)))
1416        assert num.allclose(new_points, ab_points), msg
1417
1418        geo_reference = Geo_reference(56, 100, 200)
1419        ab_points = [[1.0, 2.1], [3.0, 5.3]]
1420        points = geo_reference.change_points_geo_ref(ab_points)
1421        attributes = [2, 4]
1422        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
1423        new_geospatial = ensure_geospatial(G)
1424        new_points = new_geospatial.get_data_points(absolute=True)
1425        assert num.allclose(new_points, ab_points)
1426
1427    def test_isinstance(self):
1428        fileName = tempfile.mktemp('.csv')
1429        file = open(fileName, 'w')
1430        file.write('x,y,  elevation ,  speed \n\
14311.0, 0.0, 10.0, 0.0\n\
14320.0, 1.0, 0.0, 10.0\n\
14331.0, 0.0, 10.4, 40.0\n')
1434        file.close()
1435
1436        results = Geospatial_data(fileName)
1437        assert num.allclose(results.get_data_points(absolute=True),
1438                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
1439        assert num.allclose(results.get_attributes(attribute_name='elevation'),
1440                            [10.0, 0.0, 10.4])
1441        assert num.allclose(results.get_attributes(attribute_name='speed'),
1442                            [0.0, 10.0, 40.0])
1443
1444        os.remove(fileName)
1445
1446    def test_no_constructors(self):
1447        try:
1448            G = Geospatial_data()
1449        except ValueError:
1450            pass
1451        else:
1452            msg = 'Instance must have a filename or data points'
1453            raise Exception, msg
1454
1455    def test_load_csv_lat_long(self):
1456        '''comma delimited'''
1457
1458        fileName = tempfile.mktemp('.csv')
1459        file = open(fileName, 'w')
1460        file.write('long,lat, elevation, yeah \n\
1461150.916666667,-34.50,452.688000, 10\n\
1462150.0,-34,459.126000, 10\n')
1463        file.close()
1464
1465        results = Geospatial_data(fileName)
1466        os.remove(fileName)
1467        points = results.get_data_points()
1468
1469        assert num.allclose(points[0][0], 308728.009)
1470        assert num.allclose(points[0][1], 6180432.601)
1471        assert num.allclose(points[1][0],  222908.705)
1472        assert num.allclose(points[1][1], 6233785.284)
1473
1474    def test_load_csv_lat_longII(self):
1475        '''comma delimited'''
1476
1477        fileName = tempfile.mktemp('.csv')
1478        file = open(fileName, 'w')
1479        file.write('Lati,LONG,z \n\
1480-34.50,150.916666667,452.688000\n\
1481-34,150.0,459.126000\n')
1482        file.close()
1483
1484        results = Geospatial_data(fileName)
1485        os.remove(fileName)
1486        points = results.get_data_points()
1487
1488        assert num.allclose(points[0][0], 308728.009)
1489        assert num.allclose(points[0][1], 6180432.601)
1490        assert num.allclose(points[1][0], 222908.705)
1491        assert num.allclose(points[1][1], 6233785.284)
1492
1493    def test_load_csv_lat_long_bad(self):
1494        '''comma delimited'''
1495
1496        fileName = tempfile.mktemp('.csv')
1497        file = open(fileName, 'w')
1498        file.write('Lati,LONG,z \n\
1499-25.0,180.0,452.688000\n\
1500-34,150.0,459.126000\n')
1501        file.close()
1502
1503        try:
1504            results = Geospatial_data(fileName)
1505        except ANUGAError:
1506            pass
1507        else:
1508            msg = 'Different zones in Geo references not caught.'
1509            raise Exception, msg
1510
1511        os.remove(fileName)
1512
1513    def test_lat_long(self):
1514        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1515        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1516
1517        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1518        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1519
1520        lats = [lat_gong, lat_2]
1521        longs = [lon_gong, lon_2]
1522        gsd = Geospatial_data(latitudes=lats, longitudes=longs)
1523
1524        points = gsd.get_data_points(absolute=True)
1525
1526        assert num.allclose(points[0][0], 308728.009)
1527        assert num.allclose(points[0][1], 6180432.601)
1528        assert num.allclose(points[1][0], 222908.705)
1529        assert num.allclose(points[1][1], 6233785.284)
1530        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1531                        'Bad zone error!')
1532
1533        # use self.failUnlessRaises(ValueError, Geospatial_data(latitudes=lats))
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(latitudes=lats)
1542        except ValueError:
1543            pass
1544        else:
1545            self.fail('Error not thrown error!')
1546        try:
1547            results = Geospatial_data(longitudes=lats)
1548        except ValueError:
1549            pass
1550        else:
1551            self.fail('Error not thrown error!')
1552        try:
1553            results = Geospatial_data(latitudes=lats, longitudes=longs,
1554                                      geo_reference="p")
1555        except ValueError:
1556            pass
1557        else:
1558            self.fail('Error not thrown error!')
1559
1560        try:
1561            results = Geospatial_data(latitudes=lats, longitudes=longs,
1562                                      data_points=12)
1563        except ValueError:
1564            pass
1565        else:
1566            self.fail('Error not thrown error!')
1567
1568    def test_lat_long2(self):
1569        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1570        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1571
1572        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1573        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1574
1575        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1576        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1577
1578        points = gsd.get_data_points(absolute=True)
1579
1580        assert num.allclose(points[0][0], 308728.009)
1581        assert num.allclose(points[0][1], 6180432.601)
1582        assert num.allclose(points[1][0], 222908.705)
1583        assert num.allclose(points[1][1], 6233785.284)
1584        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1585                        'Bad zone error!')
1586
1587        try:
1588            results = Geospatial_data(points_are_lats_longs=True)
1589        except ValueError:
1590            pass
1591        else:
1592            self.fail('Error not thrown error!')
1593
1594    def test_write_urs_file(self):
1595        lat_gong = degminsec2decimal_degrees(-34, 00, 0)
1596        lon_gong = degminsec2decimal_degrees(150, 30, 0.)
1597
1598        lat_2 = degminsec2decimal_degrees(-34, 00, 1)
1599        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1600        p1 = (lat_gong, lon_gong)
1601        p2 = (lat_2, lon_2)
1602        points = frozenset([p1, p2, p1])
1603        gsd = Geospatial_data(data_points=list(points),
1604                              points_are_lats_longs=True)
1605
1606        fn = tempfile.mktemp('.urs')
1607        gsd.export_points_file(fn)
1608        handle = open(fn)
1609        lines = handle.readlines()
1610        assert lines[0], '2'
1611        assert lines[1], '-34.0002778 150.0 0'
1612        assert lines[2], '-34.0 150.5 1'
1613        handle.close()
1614        os.remove(fn)
1615
1616    def test_lat_long_set(self):
1617        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1618        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1619
1620        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1621        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1622        p1 = (lat_gong, lon_gong)
1623        p2 = (lat_2, lon_2)
1624        points = frozenset([p1, p2, p1])
1625        gsd = Geospatial_data(data_points=list(points),
1626                              points_are_lats_longs=True)
1627
1628        points = gsd.get_data_points(absolute=True)
1629        #Note the order is unknown, due to using sets
1630        # and it changes from windows to linux
1631        try:
1632            assert num.allclose(points[1][0], 308728.009)
1633            assert num.allclose(points[1][1], 6180432.601)
1634            assert num.allclose(points[0][0], 222908.705)
1635            assert num.allclose(points[0][1], 6233785.284)
1636        except AssertionError:
1637            assert num.allclose(points[0][0], 308728.009)
1638            assert num.allclose(points[0][1], 6180432.601)
1639            assert num.allclose(points[1][0], 222908.705)
1640            assert num.allclose(points[1][1], 6233785.284)
1641
1642        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1643                        'Bad zone error!')
1644        points = gsd.get_data_points(as_lat_long=True)
1645        try:
1646            assert num.allclose(points[0][0], -34)
1647            assert num.allclose(points[0][1], 150)
1648        except AssertionError:
1649            assert num.allclose(points[1][0], -34)
1650            assert num.allclose(points[1][1], 150)
1651
1652    def test_len(self):
1653        points = [[1.0, 2.1], [3.0, 5.3]]
1654        G = Geospatial_data(points)
1655        self.failUnless(2 == len(G), 'Len error!')
1656
1657        points = [[1.0, 2.1]]
1658        G = Geospatial_data(points)
1659        self.failUnless(1 == len(G), 'Len error!')
1660
1661        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1662        G = Geospatial_data(points)
1663        self.failUnless(4 == len(G), 'Len error!')
1664
1665    def test_split(self):
1666        """test if the results from spilt are disjoin sets"""
1667
1668        # below is a workaround until randint works on cyclones compute nodes
1669        if get_host_name()[8:9] != '0':
1670            points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
1671                      [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
1672                      [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0],
1673                      [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
1674                      [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
1675            attributes = {'depth': [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                          'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
1679                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
1680                                    21, 22, 23, 24, 25]}
1681            G = Geospatial_data(points, attributes)
1682
1683            factor = 0.21
1684
1685            # will return G1 with 10% of points and G2 with 90%
1686            G1, G2  = G.split(factor, 100)
1687            assert num.allclose(len(G), len(G1)+len(G2))
1688            assert num.allclose(round(len(G)*factor), len(G1))
1689
1690            P = G1.get_data_points(absolute=False)
1691            expected = [[5.,4.], [4.,1.], [2.,4.], [2.,3.], [1.,4.]]
1692            #            expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0],
1693                                     #            [3.0, 1.0], [2.0, 3.0]]
1694            msg = 'Expected %s, but\nP=%s' % (str(expected), str(P))
1695            assert num.allclose(P, expected), msg
1696
1697            A = G1.get_attributes()
1698            expected = [24, 16, 9, 8, 4]
1699            msg = 'expected=%s, but A=%s' % (str(expected), str(A))
1700            assert num.allclose(A, expected), msg
1701
1702    def test_split1(self):
1703        """test if the results from split are disjoint sets"""
1704
1705        # below is a workaround until randint works on cyclones compute nodes
1706        if get_host_name()[8:9] != '0':
1707            from numpy.random import randint, seed
1708
1709            seed((100, 100))
1710            a_points = randint(0, 999999, (10,2))
1711            points = a_points.tolist()
1712
1713            G = Geospatial_data(points)
1714
1715            factor = 0.1
1716
1717            # will return G1 with 10% of points and G2 with 90%
1718            G1, G2  = G.split(factor, 100)
1719            assert num.allclose(len(G), len(G1)+len(G2))
1720            assert num.allclose(round(len(G)*factor), len(G1))
1721
1722            P = G1.get_data_points(absolute=False)
1723            expected = [[425835., 137518.]]
1724            msg = 'expected=%s, but\nP=%s' % (str(expected), str(P))
1725            assert num.allclose(P, expected), msg
1726
1727    def test_find_optimal_smoothing_parameter(self):
1728        """
1729        Creates a elevation file representing a hill (sort of) and runs
1730        find_optimal_smoothing_parameter for 3 different alphas,
1731
1732        NOTE the random number seed is provided to control the results
1733        """
1734
1735        from cmath import cos
1736
1737        # below is a workaround until randint works on cyclones compute nodes
1738        if get_host_name()[8:9]!='0':
1739            filename = tempfile.mktemp('.csv')
1740            file = open(filename, 'w')
1741            file.write('x,y,elevation \n')
1742
1743            for i in range(-5, 6):
1744                for j in range(-5, 6):
1745                    # this equation makes a surface like a circle ripple
1746                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1747                    file.write("%s, %s, %s\n" % (i, j, z))
1748
1749            file.close()
1750
1751            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1752                                                 alpha_list=[0.0001, 0.01, 1],
1753                                                 mesh_file=None,
1754                                                 mesh_resolution=3,
1755                                                 north_boundary=5,
1756                                                 south_boundary=-5,
1757                                                 east_boundary=5,
1758                                                 west_boundary=-5,
1759                                                 plot_name=None,
1760                                                 seed_num=100000,
1761                                                 verbose=False)
1762
1763            os.remove(filename)
1764
1765            # print value, alpha
1766            assert(alpha == 0.01)
1767
1768    def test_find_optimal_smoothing_parameter1(self):
1769        """
1770        Creates a elevation file representing  a hill (sort of) and
1771        Then creates a mesh file and passes the mesh file and the elevation
1772        file to find_optimal_smoothing_parameter for 3 different alphas,
1773
1774        NOTE the random number seed is provided to control the results
1775        """
1776
1777        # below is a workaround until randint works on cyclones compute nodes
1778        if get_host_name()[8:9]!='0':
1779            from cmath import cos
1780            from anuga.pmesh.mesh_interface import create_mesh_from_regions
1781
1782            filename = tempfile.mktemp('.csv')
1783            file = open(filename, 'w')
1784            file.write('x,y,elevation \n')
1785
1786            for i in range(-5 ,6):
1787                for j in range(-5, 6):
1788                    # this equation makes a surface like a circle ripple
1789                    z = abs(cos(((i*i) + (j*j))*.1)*2)
1790                    file.write('%s, %s, %s\n' % (i, j, z))
1791
1792            file.close()
1793
1794            poly=[[5,5], [5,-5], [-5,-5], [-5,5]]
1795            internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]]
1796            mesh_filename= tempfile.mktemp('.msh')
1797
1798            create_mesh_from_regions(poly,
1799                                     boundary_tags={'back': [2],
1800                                                    'side': [1,3],
1801                                                    'ocean': [0]},
1802                                     maximum_triangle_area=3,
1803                                     interior_regions=internal_poly,
1804                                     filename=mesh_filename,
1805                                     use_cache=False,
1806                                     verbose=False)
1807
1808            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1809                                                 alpha_list=[0.0001, 0.01, 1],
1810                                                 mesh_file=mesh_filename,
1811                                                 plot_name=None,
1812                                                 seed_num=174,
1813                                                 verbose=False)
1814
1815            os.remove(filename)
1816            os.remove(mesh_filename)
1817
1818            msg = 'alpha=%s' % str(alpha)
1819            # 0.01 was expected with Numeric.RandomArray RNG
1820            assert alpha==1.0, msg
1821
1822    def test_find_optimal_smoothing_parameter2(self):
1823        '''Tests requirement that mesh file must exist or IOError is thrown
1824
1825        NOTE the random number seed is provided to control the results
1826        '''
1827
1828        from cmath import cos
1829        from anuga.pmesh.mesh_interface import create_mesh_from_regions
1830
1831        filename = tempfile.mktemp('.csv')
1832        mesh_filename= tempfile.mktemp('.msh')
1833
1834        try:
1835            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
1836                                             alpha_list=[0.0001, 0.01, 1],
1837                                             mesh_file=mesh_filename,
1838                                             plot_name=None,
1839                                             seed_num=174,
1840                                             verbose=False)
1841        except IOError:
1842            pass
1843        else:
1844            self.fail('Error not thrown error!')
1845
1846################################################################################
1847
1848if __name__ == "__main__":
1849    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1850    runner = unittest.TextTestRunner() #verbosity=2)
1851    runner.run(suite)
1852
1853
Note: See TracBrowser for help on using the repository browser.