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

Last change on this file since 6488 was 6488, checked in by ole, 15 years ago

Added new test that revealed another fitting problem where
matrix AtA isn't being built. This test was derived from
the JJKelly example by Petar Milevsky.

The test has been disabled with a preciding X

See ticket:314 for more info

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