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

Last change on this file since 6883 was 6768, checked in by rwilson, 16 years ago

Removed use of Numeric.Randomarray.

File size: 68.1 KB
RevLine 
[6086]1#!/usr/bin/env python
2
3import unittest
4import os
[6304]5import tempfile
[6086]6from math import sqrt, pi
7from sets import ImmutableSet
8
[6304]9import numpy as num
[6153]10
[6086]11from anuga.geospatial_data.geospatial_data import *
12from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
13from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
14from anuga.utilities.anuga_exceptions import ANUGAError
15from anuga.utilities.system_tools import get_host_name
16from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
[6304]17from anuga.config import netcdf_float
[6086]18
[6304]19
[6086]20class Test_Geospatial_data(unittest.TestCase):
21    def setUp(self):
22        pass
23
24    def tearDown(self):
25        pass
26
27    def test_0(self):
28        #Basic points
29        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]30
[6086]31        points = [[1.0, 2.1], [3.0, 5.3]]
32        G = Geospatial_data(points)
[6153]33        assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
[6086]34
35        # Check __repr__
36        # FIXME (Ole): Is this really machine independent?
37        rep = `G`
38        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
[6304]39        msg = 'Representation %s is not equal to %s' % (rep, ref)
[6086]40        assert rep == ref, msg
41
42        #Check getter
[6153]43        assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
[6304]44
[6086]45        #Check defaults
46        assert G.attributes is None
47        assert G.geo_reference.zone == Geo_reference().zone
48        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
49        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
50
51    def test_1(self):
52        points = [[1.0, 2.1], [3.0, 5.3]]
53        attributes = [2, 4]
[6304]54        G = Geospatial_data(points, attributes)
[6086]55        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
[6153]56        assert num.allclose(G.attributes.values()[0], [2, 4])
[6086]57
58    def test_2(self):
59        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]60
[6086]61        points = [[1.0, 2.1], [3.0, 5.3]]
62        attributes = [2, 4]
63        G = Geospatial_data(points, attributes,
64                            geo_reference=Geo_reference(56, 100, 200))
65
66        assert G.geo_reference.zone == 56
67        assert G.geo_reference.xllcorner == 100
68        assert G.geo_reference.yllcorner == 200
69
70    def test_get_attributes_1(self):
71        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]72
[6086]73        points = [[1.0, 2.1], [3.0, 5.3]]
74        attributes = [2, 4]
75        G = Geospatial_data(points, attributes,
76                            geo_reference=Geo_reference(56, 100, 200))
77
78        P = G.get_data_points(absolute=False)
[6304]79        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
[6086]80
81        P = G.get_data_points(absolute=True)
[6304]82        assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])
[6086]83
84        V = G.get_attributes() #Simply get them
[6153]85        assert num.allclose(V, [2, 4])
[6086]86
87        V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name
[6153]88        assert num.allclose(V, [2, 4])
[6086]89
90    def test_get_attributes_2(self):
91        #Multiple attributes
92        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]93
[6086]94        points = [[1.0, 2.1], [3.0, 5.3]]
95        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
96        G = Geospatial_data(points, attributes,
97                            geo_reference=Geo_reference(56, 100, 200),
98                            default_attribute_name='a1')
99
[6304]100        P = G.get_data_points(absolute=False)
101        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
[6086]102
103        V = G.get_attributes() #Get default attribute
[6153]104        assert num.allclose(V, [2, 4])
[6086]105
106        V = G.get_attributes('a0') #Get by name
[6153]107        assert num.allclose(V, [0, 0])
[6086]108
109        V = G.get_attributes('a1') #Get by name
[6153]110        assert num.allclose(V, [2, 4])
[6086]111
112        V = G.get_attributes('a2') #Get by name
[6153]113        assert num.allclose(V, [79.4, -7])
[6086]114
115        try:
116            V = G.get_attributes('hdnoatedu') #Invalid
117        except AssertionError:
118            pass
119        else:
[6304]120            raise Exception, 'Should have raised exception'
[6086]121
122    def test_get_data_points(self):
[6304]123        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
[6086]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)
[6304]128
[6086]129        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
130        results = spatial.get_data_points(absolute=False)
[6153]131        assert num.allclose(results, points_rel)
[6304]132
[6086]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)
[6304]137        results = spatial.get_data_points(geo_reference=geo_ref)
138
[6153]139        assert num.allclose(results, points_rel)
[6086]140
141    def test_get_data_points_lat_long(self):
[6304]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
[6086]148
[6304]149        points_Lat_long = [[-30., 130], [-32, 131]]
150
151        spatial = Geospatial_data(latitudes=[-30, -32.], longitudes=[130, 131])
[6086]152        results = spatial.get_data_points(as_lat_long=True)
[6153]153        assert num.allclose(results, points_Lat_long)
[6304]154
[6086]155    def test_get_data_points_lat_longII(self):
156        # x,y  North,east long,lat
157        boundary_polygon = [[ 250000, 7630000]]
158        zone = 50
[6304]159
[6086]160        geo_reference = Geo_reference(zone=zone)
[6304]161        geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference)
[6086]162        seg_lat_long = geo.get_data_points(as_lat_long=True)
[6304]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
[6086]167
168    def test_get_data_points_lat_longIII(self):
169        # x,y  North,east long,lat
[6304]170        # for northern hemisphere
[6086]171        boundary_polygon = [[419944.8, 918642.4]]
172        zone = 47
[6304]173
[6086]174        geo_reference = Geo_reference(zone=zone)
[6304]175        geo = Geospatial_data(boundary_polygon, geo_reference=geo_reference)
[6086]176        seg_lat_long = geo.get_data_points(as_lat_long=True,
177                                           isSouthHemisphere=False)
[6304]178        lat_result = degminsec2decimal_degrees(8.31, 0, 0)
179        long_result = degminsec2decimal_degrees(98.273, 0, 0)
[6153]180        assert num.allclose(seg_lat_long[0][0], lat_result)#lat
181        assert num.allclose(seg_lat_long[0][1], long_result)#long
[6086]182
[6304]183    def test_set_geo_reference(self):
184        '''test_set_georeference
[6086]185
[6304]186        Test that georeference can be changed without changing the
[6086]187        absolute values.
[6304]188        '''
189
190        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
[6086]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)
[6304]195
[6086]196        # Create without geo_ref properly set
[6304]197        G = Geospatial_data(points_rel)
[6153]198        assert not num.allclose(points_ab, G.get_data_points(absolute=True))
[6304]199
[6086]200        # Create the way it should be
201        G = Geospatial_data(points_rel, geo_reference=geo_ref)
[6153]202        assert num.allclose(points_ab, G.get_data_points(absolute=True))
[6304]203
[6086]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)
[6410]209        msg = ('points_ab=\n%s\nbut G.get_data_points(absolute=True)=\n%s'
[6304]210               % (str(points_ab), str(G.get_data_points(absolute=True))))
211        assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg
[6086]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
[6304]216
[6086]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')
[6304]226        assert points_dict.has_key('attributelist')
[6086]227        assert points_dict.has_key('geo_reference')
228
[6153]229        assert num.allclose( points_dict['pointlist'], points )
[6086]230
231        A = points_dict['attributelist']
232        assert A.has_key('a0')
233        assert A.has_key('a1')
[6304]234        assert A.has_key('a2')
[6086]235
[6153]236        assert num.allclose( A['a0'], [0, 0] )
[6304]237        assert num.allclose( A['a1'], [2, 4] )
[6153]238        assert num.allclose( A['a2'], [79.4, -7] )
[6086]239
240        geo = points_dict['geo_reference']
241        assert geo is G.geo_reference
242
243    def test_conversions_from_points_dict(self):
[6304]244        '''test conversions from points_dict'''
[6086]245
246        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]247
[6086]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)
[6304]258        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
[6086]259
260        V = G.get_attributes('a0') #Get by name
[6153]261        assert num.allclose(V, [0, 0])
[6086]262
263        V = G.get_attributes('a1') #Get by name
[6153]264        assert num.allclose(V, [2, 4])
[6086]265
266        V = G.get_attributes('a2') #Get by name
[6153]267        assert num.allclose(V, [79.4, -7])
[6086]268
269    def test_add(self):
[6304]270        '''test the addition of two geospatical objects
271        no geo_reference see next test
272        '''
273
[6086]274        points = [[1.0, 2.1], [3.0, 5.3]]
[6304]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
[6086]280        G = G1 + G2
281
282        assert G.attributes.has_key('depth')
283        assert G.attributes.has_key('elevation')
[6153]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
[6086]289    def test_addII(self):
[6304]290        '''test the addition of two geospatical objects
291        no geo_reference see next test
292        '''
293
[6086]294        points = [[1.0, 2.1], [3.0, 5.3]]
[6304]295        attributes = {'depth': [2, 4]}
296        G1 = Geospatial_data(points, attributes)
297
[6086]298        points = [[5.0, 2.1], [3.0, 50.3]]
[6304]299        attributes = {'depth': [200, 400]}
[6086]300        G2 = Geospatial_data(points, attributes)
[6304]301
[6086]302        G = G1 + G2
303
[6304]304        assert G.attributes.has_key('depth')
[6086]305        assert G.attributes.keys(), ['depth']
[6153]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]])
[6304]309
[6086]310    def test_add_with_geo (self):
[6304]311        '''Difference in Geo_reference resolved'''
312
[6086]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')
[6304]324
[6086]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)
[6153]330        assert num.allclose(P1, [[2.0, 4.1], [4.0, 7.3]])
[6086]331
332        P2 = G2.get_data_points(absolute=True)
[6304]333        assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]])
334
[6086]335        G = G1 + G2
336
337        # Check absoluteness
[6153]338        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
339        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6086]340
341        P = G.get_data_points(absolute=True)
342
[6428]343        assert num.allclose(P, num.concatenate((P1,P2), axis=0))    #??default#
[6410]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))
[6153]348        assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3],
[6304]349                                [5.1, 9.1], [6.1, 6.3]]), msg
[6086]350
[6304]351    def test_add_with_geo_absolute (self):
352        '''Difference in Geo_reference resolved'''
[6086]353
[6153]354        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
[6304]355        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
[6086]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
[6304]361        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(),
362                                        geo_ref1.get_yllcorner()],
[6086]363                             attributes1, geo_ref1)
[6304]364
365        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(),
366                                        geo_ref2.get_yllcorner()],
[6086]367                             attributes2, geo_ref2)
368
369        #Check that absolute values are as expected
370        P1 = G1.get_data_points(absolute=True)
[6153]371        assert num.allclose(P1, points1)
[6086]372
373        P1 = G1.get_data_points(absolute=False)
[6410]374        msg = ('P1=\n%s\nbut points1 - <...>=\n%s'
[6304]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
[6086]380
381        P2 = G2.get_data_points(absolute=True)
[6153]382        assert num.allclose(P2, points2)
[6086]383
384        P2 = G2.get_data_points(absolute=False)
[6304]385        assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(),
386                                           geo_ref2.get_yllcorner()])
387
[6086]388        G = G1 + G2
389        P = G.get_data_points(absolute=True)
390
[6428]391        assert num.allclose(P, num.concatenate((points1, points2), axis=0)) #??default#
[6086]392
[6304]393    def test_add_with_None(self):
394        '''test that None can be added to a geospatical objects'''
[6086]395
[6153]396        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
[6304]397        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
[6086]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
[6304]407        attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]}
408        attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]}
[6086]409
410        G1 = Geospatial_data(points1, attributes1, geo_ref1)
[6153]411        assert num.allclose(G1.get_geo_reference().get_xllcorner(), 1.0)
412        assert num.allclose(G1.get_geo_reference().get_yllcorner(), 2.0)
[6086]413        assert G1.attributes.has_key('depth')
414        assert G1.attributes.has_key('elevation')
[6153]415        assert num.allclose(G1.attributes['depth'], [2, 4.7])
[6304]416        assert num.allclose(G1.attributes['elevation'], [6.1, 5])
417
[6086]418        G2 = Geospatial_data(points2, attributes2, geo_ref2)
[6153]419        assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1)
420        assert num.allclose(G2.get_geo_reference().get_yllcorner(), 3.0)
[6086]421        assert G2.attributes.has_key('depth')
422        assert G2.attributes.has_key('elevation')
[6153]423        assert num.allclose(G2.attributes['depth'], [-2.3, 4])
[6304]424        assert num.allclose(G2.attributes['elevation'], [2.5, 1])
[6086]425
426        #Check that absolute values are as expected
427        P1 = G1.get_data_points(absolute=True)
[6153]428        assert num.allclose(P1, [[3.0, 6.1], [5.0, 9.3]])
[6086]429
430        P2 = G2.get_data_points(absolute=True)
[6304]431        assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])
[6086]432
433        # Normal add
434        G = G1 + None
435        assert G.attributes.has_key('depth')
436        assert G.attributes.has_key('elevation')
[6153]437        assert num.allclose(G.attributes['depth'], [2, 4.7])
[6304]438        assert num.allclose(G.attributes['elevation'], [6.1, 5])
[6086]439
440        # Points are now absolute.
[6153]441        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
442        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6304]443        P = G.get_data_points(absolute=True)
[6410]444        msg = 'P=\n%s' % str(P)
[6304]445        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg
[6086]446
447        G = G2 + None
448        assert G.attributes.has_key('depth')
449        assert G.attributes.has_key('elevation')
[6153]450        assert num.allclose(G.attributes['depth'], [-2.3, 4])
[6304]451        assert num.allclose(G.attributes['elevation'], [2.5, 1])
[6153]452        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
453        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6304]454
455        P = G.get_data_points(absolute=True)
[6153]456        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
[6086]457
458        # Reverse add
459        G = None + G1
460        assert G.attributes.has_key('depth')
461        assert G.attributes.has_key('elevation')
[6153]462        assert num.allclose(G.attributes['depth'], [2, 4.7])
[6304]463        assert num.allclose(G.attributes['elevation'], [6.1, 5])
[6086]464
465        # Points are now absolute.
[6153]466        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
467        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6086]468
[6304]469        P = G.get_data_points(absolute=True)
470        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
[6086]471
472        G = None + G2
473        assert G.attributes.has_key('depth')
474        assert G.attributes.has_key('elevation')
[6153]475        assert num.allclose(G.attributes['depth'], [-2.3, 4])
[6304]476        assert num.allclose(G.attributes['elevation'], [2.5, 1])
[6086]477
[6153]478        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
479        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6304]480
481        P = G.get_data_points(absolute=True)
[6153]482        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
[6086]483
[6304]484    def test_clip0(self):
485        '''test_clip0(self):
[6086]486
487        Test that point sets can be clipped by a polygon
[6304]488        '''
489
[6086]490        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]491
[6086]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
[6304]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]])
[6086]500
501        # Then a more complex polygon
502        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[6304]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]]
[6086]505        G = Geospatial_data(points)
506
[6153]507        assert num.allclose(G.clip(polygon).get_data_points(),
[6304]508                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
[6086]509
510    def test_clip0_with_attributes(self):
[6304]511        '''test_clip0_with_attributes(self):
512
[6086]513        Test that point sets with attributes can be clipped by a polygon
[6304]514        '''
515
[6086]516        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]517
[6086]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,
[6153]523                    'att2': num.array(attributes)+1}
[6304]524
[6086]525        G = Geospatial_data(points, att_dict)
526
[6304]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]])
[6153]531        assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
[6304]532        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
[6086]533
534        # Then a more complex polygon
535        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[6304]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]]
[6086]538
539        # This time just one attribute
540        attributes = [2, -4, 5, 76, -2, 0.1]
541        G = Geospatial_data(points, attributes)
[6153]542        assert num.allclose(G.clip(polygon).get_data_points(),
[6304]543                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
[6153]544        assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
[6086]545
546    def test_clip1(self):
[6304]547        '''test_clip1(self):
548
[6086]549        Test that point sets can be clipped by a polygon given as
550        another Geospatial dataset
[6304]551        '''
552
[6086]553        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]554
[6086]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,
[6153]559                    'att2': num.array(attributes)+1}
[6086]560        G = Geospatial_data(points, att_dict)
[6304]561
562        # First try the unit square
563        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
[6153]564        assert num.allclose(G.clip(U).get_data_points(),
[6304]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])
[6086]568
569        # Then a more complex polygon
[6304]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]
[6086]573        G = Geospatial_data(points, attributes)
[6304]574        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1],
575                                   [2,1], [0,1]])
[6086]576
[6153]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])
[6086]580
581    def test_clip0_outside(self):
[6304]582        '''test_clip0_outside(self):
583
[6086]584        Test that point sets can be clipped outside of a polygon
[6304]585        '''
586
[6086]587        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]588
[6086]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]]
[6304]591        attributes = [2, -4, 5, 76, -2, 0.1, 3]
[6086]592        G = Geospatial_data(points, attributes)
593
[6304]594        # First try the unit square
[6086]595        U = [[0,0], [1,0], [1,1], [0,1]]
[6153]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]])
[6304]598        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
[6086]599
600        # Then a more complex polygon
601        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
[6304]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]
[6086]605        G = Geospatial_data(points, attributes)
[6153]606        assert num.allclose(G.clip_outside(polygon).get_data_points(),
607                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
[6304]608        assert num.allclose(G.clip_outside(polygon).get_attributes(),
609                            [2, -2, 0.1])
[6086]610
[6304]611    def test_clip1_outside(self):
612        '''test_clip1_outside(self):
[6086]613
614        Test that point sets can be clipped outside of a polygon given as
615        another Geospatial dataset
[6304]616        '''
617
[6086]618        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]619
[6086]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]]
[6304]622        attributes = [2, -4, 5, 76, -2, 0.1, 3]
623        G = Geospatial_data(points, attributes)
[6086]624
[6304]625        # First try the unit square
626        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
[6153]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]])
[6304]629        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
[6086]630
631        # Then a more complex polygon
[6304]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]
[6086]635        G = Geospatial_data(points, attributes)
[6304]636        polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1],
637                                   [2, 1], [0, 1]])
[6153]638        assert num.allclose(G.clip_outside(polygon).get_data_points(),
639                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
[6304]640        assert num.allclose(G.clip_outside(polygon).get_attributes(),
641                            [2, -2, 0.1])
[6086]642
[6304]643    def test_clip1_inside_outside(self):
644        '''test_clip1_inside_outside(self):
[6086]645
646        Test that point sets can be clipped outside of a polygon given as
647        another Geospatial dataset
[6304]648        '''
649
[6086]650        from anuga.coordinate_transforms.geo_reference import Geo_reference
[6304]651
[6086]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]]
[6304]654        attributes = [2, -4, 5, 76, -2, 0.1, 3]
[6086]655        G = Geospatial_data(points, attributes)
656
[6304]657        # First try the unit square
658        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
[6086]659        G1 = G.clip(U)
[6304]660        assert num.allclose(G1.get_data_points(),
661                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
[6153]662        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
[6086]663        G2 = G.clip_outside(U)
[6153]664        assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
665                                                  [3.0, 5.3], [2.4, 3.3]])
[6304]666        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
[6086]667
668        # New ordering
[6304]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]
[6086]672
[6153]673        assert num.allclose((G1+G2).get_data_points(), new_points)
674        assert num.allclose((G1+G2).get_attributes(), new_attributes)
[6086]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
[6304]683        # Check result
684        assert num.allclose(G3.get_data_points(), new_points)
685        assert num.allclose(G3.get_attributes(), new_attributes)
[6086]686
687        os.remove(FN)
688
689    def test_load_csv(self):
[6304]690        fileName = tempfile.mktemp('.csv')
691        file = open(fileName,'w')
692        file.write('x,y,elevation speed \n\
[6086]6931.0 0.0 10.0 0.0\n\
6940.0 1.0 0.0 10.0\n\
[6304]6951.0 0.0 10.4 40.0\n')
[6086]696        file.close()
[6304]697
[6086]698        results = Geospatial_data(fileName)
699        os.remove(fileName)
[6304]700        assert num.allclose(results.get_data_points(),
701                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
[6153]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])
[6086]706
[6304]707################################################################################
708# Test CSV files
709################################################################################
[6086]710
711    def test_load_csv_lat_long_bad_blocking(self):
[6304]712        '''test_load_csv_lat_long_bad_blocking(self):
[6086]713        Different zones in Geo references
[6304]714        '''
715
716        fileName = tempfile.mktemp('.csv')
717        file = open(fileName, 'w')
718        file.write('Lati,LONG,z \n\
[6086]719-25.0,180.0,452.688000\n\
[6304]720-34,150.0,459.126000\n')
[6086]721        file.close()
[6304]722
[6086]723        results = Geospatial_data(fileName, max_read_lines=1,
724                                  load_file_now=False)
[6304]725
[6086]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.'
[6304]733            raise Exception, msg
734
[6086]735        os.remove(fileName)
[6304]736
[6086]737    def test_load_csv(self):
[6304]738        fileName = tempfile.mktemp('.txt')
739        file = open(fileName, 'w')
740        file.write(' x,y, elevation ,  speed \n\
[6086]7411.0, 0.0, 10.0, 0.0\n\
7420.0, 1.0, 0.0, 10.0\n\
[6304]7431.0, 0.0 ,10.4, 40.0\n')
[6086]744        file.close()
745
746        results = Geospatial_data(fileName, max_read_lines=2)
747
[6304]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])
[6086]754
755        # Blocking
756        geo_list = []
757        for i in results:
758            geo_list.append(i)
[6304]759
[6153]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(),
[6304]765                            [[1.0, 0.0]])
[6153]766        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
767                            [10.4])
[6304]768
769        os.remove(fileName)
770
[6086]771    def test_load_csv_bad(self):
[6304]772        '''test_load_csv_bad(self):
[6086]773        header column, body column missmatch
774        (Format error)
[6304]775        '''
776
777        fileName = tempfile.mktemp('.txt')
778        file = open(fileName, 'w')
779        file.write(' elevation ,  speed \n\
[6086]7801.0, 0.0, 10.0, 0.0\n\
7810.0, 1.0, 0.0, 10.0\n\
[6304]7821.0, 0.0 ,10.4, 40.0\n')
[6086]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:
[6304]796            msg = 'Bad file did not raise error!'
797            raise Exception, msg
[6086]798        os.remove(fileName)
799
800    def test_load_csv_badII(self):
[6304]801        '''test_load_csv_bad(self):
[6086]802        header column, body column missmatch
803        (Format error)
[6304]804        '''
805
806        fileName = tempfile.mktemp('.txt')
807        file = open(fileName, 'w')
808        file.write(' x,y,elevation ,  speed \n\
[6086]8091.0, 0.0, 10.0, 0.0\n\
8100.0, 1.0, 0.0, 10\n\
[6304]8111.0, 0.0 ,10.4 yeah\n')
[6086]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:
[6304]825            msg = 'Bad file did not raise error!'
826            raise Exception, msg
[6086]827        os.remove(fileName)
828
829    def test_load_csv_badIII(self):
[6304]830        '''test_load_csv_bad(self):
[6086]831        space delimited
[6304]832        '''
833
834        fileName = tempfile.mktemp('.txt')
835        file = open(fileName, 'w')
836        file.write(' x y elevation   speed \n\
[6086]8371. 0.0 10.0 0.0\n\
8380.0 1.0 0.0 10.0\n\
[6304]8391.0 0.0 10.4 40.0\n')
[6086]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:
[6304]848            msg = 'Bad file did not raise error!'
849            raise Exception, msg
[6086]850        os.remove(fileName)
[6304]851
[6086]852    def test_load_csv_badIV(self):
[6304]853        ''' test_load_csv_bad(self):
[6086]854        header column, body column missmatch
855        (Format error)
[6304]856        '''
857
858        fileName = tempfile.mktemp('.txt')
859        file = open(fileName, 'w')
860        file.write(' x,y,elevation ,  speed \n\
[6086]8611.0, 0.0, 10.0, wow\n\
8620.0, 1.0, 0.0, ha\n\
[6304]8631.0, 0.0 ,10.4, yeah\n')
[6086]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:
[6304]877            msg = 'Bad file did not raise error!'
878            raise Exception, msg
[6086]879        os.remove(fileName)
880
881    def test_load_pts_blocking(self):
882        #This is pts!
[6304]883        fileName = tempfile.mktemp('.txt')
884        file = open(fileName, 'w')
885        file.write(' x,y, elevation ,  speed \n\
[6086]8861.0, 0.0, 10.0, 0.0\n\
8870.0, 1.0, 0.0, 10.0\n\
[6304]8881.0, 0.0 ,10.4, 40.0\n')
[6086]889        file.close()
890
[6304]891        pts_file = tempfile.mktemp('.pts')
892
[6086]893        convert = Geospatial_data(fileName)
894        convert.export_points_file(pts_file)
895        results = Geospatial_data(pts_file, max_read_lines=2)
896
[6304]897        assert num.allclose(results.get_data_points(),
898                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
[6153]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])
[6086]903
904        # Blocking
905        geo_list = []
906        for i in results:
[6304]907            geo_list.append(i)
[6153]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(),
[6304]913                            [[1.0, 0.0]])
[6153]914        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
915                            [10.4])
[6086]916
[6304]917        os.remove(fileName)
918        os.remove(pts_file)
919
[6086]920    def verbose_test_load_pts_blocking(self):
[6304]921        fileName = tempfile.mktemp('.txt')
922        file = open(fileName, 'w')
923        file.write(' x,y, elevation ,  speed \n\
[6086]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\
[6304]9701.0, 0.0 ,10.4, 40.0\n')
[6086]971        file.close()
972
[6304]973        pts_file = tempfile.mktemp('.pts')
[6086]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:
[6304]981            geo_list.append(i)
[6153]982        assert num.allclose(geo_list[0].get_data_points(),
[6304]983                            [[1.0, 0.0], [0.0, 1.0]])
[6153]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(),
[6304]987                            [[1.0, 0.0],[0.0, 1.0] ])
[6153]988        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
989                            [10.0, 0.0])
[6304]990
991        os.remove(fileName)
[6086]992        os.remove(pts_file)
993
994    def test_new_export_pts_file(self):
995        att_dict = {}
[6153]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])
[6304]999
1000        fileName = tempfile.mktemp('.pts')
[6086]1001        G = Geospatial_data(pointlist, att_dict)
1002        G.export_points_file(fileName)
1003        results = Geospatial_data(file_name = fileName)
[6304]1004        os.remove(fileName)
[6086]1005
[6304]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])
[6086]1010        answer = [10.0, 1.0, 10.4]
[6304]1011        assert num.allclose(results.get_attributes(attribute_name='brightness'),
1012                            answer)
[6086]1013
1014    def test_new_export_absolute_pts_file(self):
1015        att_dict = {}
[6304]1016        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]1017        att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
1018        att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
[6086]1019        geo_ref = Geo_reference(50, 25, 55)
[6304]1020
1021        fileName = tempfile.mktemp('.pts')
[6086]1022        G = Geospatial_data(pointlist, att_dict, geo_ref)
1023        G.export_points_file(fileName, absolute=True)
1024        results = Geospatial_data(file_name = fileName)
[6304]1025        os.remove(fileName)
[6086]1026
[6410]1027        msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s'
[6304]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
[6086]1035        answer = [10.0, 1.0, 10.4]
[6304]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
[6086]1042
1043    def test_loadpts(self):
1044        from Scientific.IO.NetCDF import NetCDFFile
1045
[6304]1046        fileName = tempfile.mktemp('.pts')
[6086]1047        # NetCDF file definition
1048        outfile = NetCDFFile(fileName, netcdf_mode_w)
[6304]1049
[6086]1050        # dimension definitions
[6304]1051        outfile.createDimension('number_of_points', 3)
1052        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1053
[6086]1054        # variable definitions
[6304]1055        outfile.createVariable('points', netcdf_float, ('number_of_points',
1056                                                        'number_of_dimensions'))
1057        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1058
[6086]1059        # Get handles to the variables
1060        points = outfile.variables['points']
1061        elevation = outfile.variables['elevation']
[6304]1062
[6086]1063        points[0, :] = [1.0,0.0]
[6304]1064        elevation[0] = 10.0
[6086]1065        points[1, :] = [0.0,1.0]
[6304]1066        elevation[1] = 0.0
[6086]1067        points[2, :] = [1.0,0.0]
[6304]1068        elevation[2] = 10.4
[6086]1069
1070        outfile.close()
[6304]1071
[6086]1072        results = Geospatial_data(file_name = fileName)
1073        os.remove(fileName)
[6304]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
[6086]1081    def test_writepts(self):
[6304]1082        '''Test that storage of x,y,attributes works'''
1083
[6086]1084        att_dict = {}
[6304]1085        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]1086        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1087        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
[6304]1088        geo_reference=Geo_reference(56, 1.9, 1.9)
[6086]1089
1090        # Test pts format
[6304]1091        fileName = tempfile.mktemp('.pts')
[6086]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
[6304]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])
[6086]1101        answer = [10.0, 0.0, 10.4]
[6153]1102        assert num.allclose(results.get_attributes('brightness'), answer)
[6086]1103        self.failUnless(geo_reference == geo_reference,
[6304]1104                        'test_writepts failed. Test geo_reference')
[6086]1105
1106    def test_write_csv_attributes(self):
[6304]1107        '''Test that storage of x,y,attributes works'''
1108
[6086]1109        att_dict = {}
[6304]1110        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]1111        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1112        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
[6304]1113        geo_reference=Geo_reference(56, 0, 0)
1114
[6086]1115        # Test txt format
[6304]1116        fileName = tempfile.mktemp('.txt')
[6086]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)
[6304]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])
[6086]1126        answer = [10.0, 0.0, 10.4]
[6153]1127        assert num.allclose(results.get_attributes('brightness'), answer)
[6304]1128
[6086]1129    def test_write_csv_attributes_lat_long(self):
[6304]1130        '''Test that storage of x,y,attributes works'''
1131
[6086]1132        att_dict = {}
[6304]1133        pointlist = num.array([[-21.5, 114.5], [-21.6, 114.5], [-21.7, 114.5]])
[6153]1134        att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
1135        att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
[6304]1136
[6086]1137        # Test txt format
[6304]1138        fileName = tempfile.mktemp('.txt')
[6086]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)
[6304]1143
[6153]1144        assert num.allclose(results.get_data_points(False, as_lat_long=True),
[6304]1145                            pointlist)
1146        assert num.allclose(results.get_attributes('elevation'),
1147                            [10.0, 0.0, 10.4])
[6086]1148        answer = [10.0, 0.0, 10.4]
[6153]1149        assert num.allclose(results.get_attributes('brightness'), answer)
[6304]1150
[6086]1151    def test_writepts_no_attributes(self):
[6304]1152        '''Test that storage of x,y alone works'''
[6086]1153
1154        att_dict = {}
[6304]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)
[6086]1157
1158        # Test pts format
[6304]1159        fileName = tempfile.mktemp('.pts')
[6086]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
[6304]1165        assert num.allclose(results.get_data_points(False),
1166                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6086]1167        self.failUnless(geo_reference == geo_reference,
[6304]1168                        'test_writepts failed. Test geo_reference')
1169
[6086]1170    def test_write_csv_no_attributes(self):
[6304]1171        '''Test that storage of x,y alone works'''
1172
[6086]1173        att_dict = {}
[6304]1174        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6086]1175        geo_reference=Geo_reference(56,0,0)
[6304]1176
[6086]1177        # Test format
[6304]1178        fileName = tempfile.mktemp('.txt')
[6086]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
[6304]1184        assert num.allclose(results.get_data_points(False),
1185                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6086]1186
[6304]1187################################################################################
1188# Check bad PTS files.
1189################################################################################
1190
[6086]1191    def test_load_bad_no_file_pts(self):
[6304]1192        fileName = tempfile.mktemp('.pts')
[6086]1193        try:
[6304]1194            results = Geospatial_data(file_name=fileName)
[6086]1195        except IOError:
1196            pass
1197        else:
1198            msg = 'imaginary file did not raise error!'
[6304]1199            raise Exception, msg
[6086]1200
1201    def test_create_from_pts_file(self):
1202        from Scientific.IO.NetCDF import NetCDFFile
1203
[6304]1204        # NetCDF file definition
[6086]1205        FN = 'test_points.pts'
1206        outfile = NetCDFFile(FN, netcdf_mode_w)
[6304]1207
[6086]1208        # dimension definitions
[6304]1209        outfile.createDimension('number_of_points', 3)
1210        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1211
[6086]1212        # variable definitions
[6304]1213        outfile.createVariable('points', netcdf_float, ('number_of_points',
1214                                                        'number_of_dimensions'))
1215        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1216
[6086]1217        # Get handles to the variables
1218        points = outfile.variables['points']
1219        elevation = outfile.variables['elevation']
[6304]1220
[6086]1221        points[0, :] = [1.0,0.0]
[6304]1222        elevation[0] = 10.0
[6086]1223        points[1, :] = [0.0,1.0]
[6304]1224        elevation[1] = 0.0
[6086]1225        points[2, :] = [1.0,0.0]
[6304]1226        elevation[2] = 10.4
[6086]1227
1228        outfile.close()
1229
1230        G = Geospatial_data(file_name = FN)
1231
[6153]1232        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
1233        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
[6304]1234        assert num.allclose(G.get_data_points(),
1235                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]1236        assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4])
[6086]1237        os.remove(FN)
1238
1239    def test_create_from_pts_file_with_geo(self):
[6304]1240        '''Test if Geospatial data is correctly instantiated from a pts file.'''
1241
[6086]1242        from Scientific.IO.NetCDF import NetCDFFile
1243
[6304]1244        # NetCDF file definition
[6086]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
[6304]1255        outfile.createDimension('number_of_points', 3)
1256        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
1257
[6086]1258        # variable definitions
[6304]1259        outfile.createVariable('points', netcdf_float, ('number_of_points',
1260                                                        'number_of_dimensions'))
1261        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
1262
[6086]1263        # Get handles to the variables
1264        points = outfile.variables['points']
1265        elevation = outfile.variables['elevation']
1266
1267        points[0, :] = [1.0,0.0]
[6304]1268        elevation[0] = 10.0
[6086]1269        points[1, :] = [0.0,1.0]
[6304]1270        elevation[1] = 0.0
[6086]1271        points[2, :] = [1.0,0.0]
[6304]1272        elevation[2] = 10.4
[6086]1273
1274        outfile.close()
1275
1276        G = Geospatial_data(file_name = FN)
1277
[6153]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])
[6304]1284
[6086]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
[6304]1290        checking results are correct
[6086]1291        '''
[6304]1292
[6086]1293        # create files
1294        att_dict1 = {}
[6304]1295        pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]1296        att_dict1['elevation'] = num.array([-10.0, 0.0, 10.4])
1297        att_dict1['brightness'] = num.array([10.0, 0.0, 10.4])
[6086]1298        geo_reference1 = Geo_reference(56, 2.0, 1.0)
[6304]1299
[6086]1300        att_dict2 = {}
[6304]1301        pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0], [2.0, 1.0]])
[6153]1302        att_dict2['elevation'] = num.array([1.0, 15.0, 1.4])
1303        att_dict2['brightness'] = num.array([14.0, 1.0, -12.4])
[6304]1304        geo_reference2 = Geo_reference(56, 1.0, 2.0)
[6086]1305
1306        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
1307        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
1308
[6304]1309        fileName1 = tempfile.mktemp('.txt')
1310        fileName2 = tempfile.mktemp('.pts')
1311
1312        # makes files
[6086]1313        G1.export_points_file(fileName1)
1314        G2.export_points_file(fileName2)
[6304]1315
[6086]1316        # add files
[6304]1317        G3 = Geospatial_data(file_name=fileName1)
1318        G4 = Geospatial_data(file_name=fileName2)
[6086]1319        G = G3 + G4
1320
1321        #read results
[6153]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])
[6086]1328        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
[6304]1329        assert num.allclose(G.get_attributes(attribute_name='brightness'),
1330                            answer)
[6086]1331        self.failUnless(G.get_geo_reference() == geo_reference1,
[6304]1332                        'test_writepts failed. Test geo_reference')
1333
[6086]1334        os.remove(fileName1)
1335        os.remove(fileName2)
[6304]1336
[6086]1337    def test_ensure_absolute(self):
[6304]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]]
[6086]1341        new_points = ensure_absolute(points)
[6304]1342
[6153]1343        assert num.allclose(new_points, points)
[6086]1344
[6304]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]])
[6086]1348        new_points = ensure_absolute(points)
[6304]1349
[6153]1350        assert num.allclose(new_points, points)
[6086]1351
[6304]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
[6086]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]
[6304]1363        new_points = ensure_absolute(data_points, geo_reference=mesh_origin)
1364
[6153]1365        assert num.allclose(new_points, ab_points)
[6086]1366
1367        geo = Geo_reference(56,67,-56)
[6304]1368        data_points = geo.change_points_geo_ref(ab_points)
1369        new_points = ensure_absolute(data_points, geo_reference=geo)
[6086]1370
[6153]1371        assert num.allclose(new_points, ab_points)
[6086]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]
[6304]1377        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
[6086]1378        new_points = ensure_absolute(G)
[6304]1379
[6153]1380        assert num.allclose(new_points, ab_points)
[6086]1381
1382    def test_ensure_geospatial(self):
[6304]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]]
[6086]1386        new_points = ensure_geospatial(points)
[6304]1387        assert num.allclose(new_points.get_data_points(absolute=True), points)
[6086]1388
[6304]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]])
[6086]1392        new_points = ensure_geospatial(points)
[6304]1393        assert num.allclose(new_points.get_data_points(absolute=True), points)
1394
[6153]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]])
[6304]1398        mesh_origin = (56, 290000, 618000)      # zone, easting, northing
1399        data_points = num.zeros((ab_points.shape), num.float)
[6086]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,
[6304]1406                                           geo_reference=mesh_origin)
[6086]1407        new_points = new_geospatial.get_data_points(absolute=True)
[6153]1408        assert num.allclose(new_points, ab_points)
[6086]1409
[6304]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)
[6086]1413        new_points = new_geospatial.get_data_points(absolute=True)
[6441]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
[6086]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]
[6304]1422        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
1423        new_geospatial = ensure_geospatial(G)
[6086]1424        new_points = new_geospatial.get_data_points(absolute=True)
[6153]1425        assert num.allclose(new_points, ab_points)
[6304]1426
[6086]1427    def test_isinstance(self):
[6304]1428        fileName = tempfile.mktemp('.csv')
1429        file = open(fileName, 'w')
1430        file.write('x,y,  elevation ,  speed \n\
[6086]14311.0, 0.0, 10.0, 0.0\n\
14320.0, 1.0, 0.0, 10.0\n\
[6304]14331.0, 0.0, 10.4, 40.0\n')
[6086]1434        file.close()
1435
1436        results = Geospatial_data(fileName)
[6153]1437        assert num.allclose(results.get_data_points(absolute=True),
[6304]1438                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
[6153]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])
[6086]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'
[6304]1453            raise Exception, msg
[6086]1454
1455    def test_load_csv_lat_long(self):
[6304]1456        '''comma delimited'''
[6086]1457
[6304]1458        fileName = tempfile.mktemp('.csv')
1459        file = open(fileName, 'w')
1460        file.write('long,lat, elevation, yeah \n\
[6086]1461150.916666667,-34.50,452.688000, 10\n\
[6304]1462150.0,-34,459.126000, 10\n')
[6086]1463        file.close()
[6304]1464
[6086]1465        results = Geospatial_data(fileName)
1466        os.remove(fileName)
1467        points = results.get_data_points()
[6304]1468
[6153]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)
[6304]1473
[6086]1474    def test_load_csv_lat_longII(self):
[6304]1475        '''comma delimited'''
[6086]1476
[6304]1477        fileName = tempfile.mktemp('.csv')
1478        file = open(fileName, 'w')
1479        file.write('Lati,LONG,z \n\
[6086]1480-34.50,150.916666667,452.688000\n\
[6304]1481-34,150.0,459.126000\n')
[6086]1482        file.close()
[6304]1483
[6086]1484        results = Geospatial_data(fileName)
1485        os.remove(fileName)
1486        points = results.get_data_points()
[6304]1487
[6153]1488        assert num.allclose(points[0][0], 308728.009)
1489        assert num.allclose(points[0][1], 6180432.601)
[6304]1490        assert num.allclose(points[1][0], 222908.705)
[6153]1491        assert num.allclose(points[1][1], 6233785.284)
[6086]1492
1493    def test_load_csv_lat_long_bad(self):
[6304]1494        '''comma delimited'''
[6086]1495
[6304]1496        fileName = tempfile.mktemp('.csv')
1497        file = open(fileName, 'w')
1498        file.write('Lati,LONG,z \n\
[6086]1499-25.0,180.0,452.688000\n\
[6304]1500-34,150.0,459.126000\n')
[6086]1501        file.close()
[6304]1502
[6086]1503        try:
1504            results = Geospatial_data(fileName)
1505        except ANUGAError:
1506            pass
1507        else:
1508            msg = 'Different zones in Geo references not caught.'
[6304]1509            raise Exception, msg
1510
[6086]1511        os.remove(fileName)
[6304]1512
[6086]1513    def test_lat_long(self):
[6304]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
[6086]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)
[6304]1525
[6153]1526        assert num.allclose(points[0][0], 308728.009)
1527        assert num.allclose(points[0][1], 6180432.601)
[6304]1528        assert num.allclose(points[1][0], 222908.705)
[6153]1529        assert num.allclose(points[1][1], 6233785.284)
[6086]1530        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1531                        'Bad zone error!')
[6304]1532
[6086]1533        try:
1534            results = Geospatial_data(latitudes=lats)
1535        except ValueError:
1536            pass
1537        else:
[6304]1538            self.fail('Error not thrown error!')
[6086]1539        try:
1540            results = Geospatial_data(latitudes=lats)
1541        except ValueError:
1542            pass
1543        else:
[6304]1544            self.fail('Error not thrown error!')
[6086]1545        try:
1546            results = Geospatial_data(longitudes=lats)
1547        except ValueError:
1548            pass
1549        else:
[6304]1550            self.fail('Error not thrown error!')
[6086]1551        try:
1552            results = Geospatial_data(latitudes=lats, longitudes=longs,
1553                                      geo_reference="p")
1554        except ValueError:
1555            pass
1556        else:
[6304]1557            self.fail('Error not thrown error!')
1558
[6086]1559        try:
1560            results = Geospatial_data(latitudes=lats, longitudes=longs,
1561                                      data_points=12)
1562        except ValueError:
1563            pass
1564        else:
[6304]1565            self.fail('Error not thrown error!')
[6086]1566
1567    def test_lat_long2(self):
[6304]1568        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1569        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1570
1571        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1572        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
1573
[6086]1574        points = [[lat_gong, lon_gong], [lat_2, lon_2]]
1575        gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
1576
1577        points = gsd.get_data_points(absolute=True)
[6304]1578
[6153]1579        assert num.allclose(points[0][0], 308728.009)
1580        assert num.allclose(points[0][1], 6180432.601)
[6304]1581        assert num.allclose(points[1][0], 222908.705)
[6153]1582        assert num.allclose(points[1][1], 6233785.284)
[6086]1583        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1584                        'Bad zone error!')
1585
1586        try:
1587            results = Geospatial_data(points_are_lats_longs=True)
1588        except ValueError:
1589            pass
1590        else:
[6304]1591            self.fail('Error not thrown error!')
[6086]1592
[6304]1593    def test_write_urs_file(self):
1594        lat_gong = degminsec2decimal_degrees(-34, 00, 0)
1595        lon_gong = degminsec2decimal_degrees(150, 30, 0.)
[6086]1596
[6304]1597        lat_2 = degminsec2decimal_degrees(-34, 00, 1)
1598        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
[6086]1599        p1 = (lat_gong, lon_gong)
1600        p2 = (lat_2, lon_2)
1601        points = ImmutableSet([p1, p2, p1])
1602        gsd = Geospatial_data(data_points=list(points),
1603                              points_are_lats_longs=True)
[6304]1604
1605        fn = tempfile.mktemp('.urs')
[6086]1606        gsd.export_points_file(fn)
1607        handle = open(fn)
1608        lines = handle.readlines()
[6304]1609        assert lines[0], '2'
1610        assert lines[1], '-34.0002778 150.0 0'
1611        assert lines[2], '-34.0 150.5 1'
[6086]1612        handle.close()
1613        os.remove(fn)
[6304]1614
[6086]1615    def test_lat_long_set(self):
[6304]1616        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
1617        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
1618
1619        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
1620        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
[6086]1621        p1 = (lat_gong, lon_gong)
1622        p2 = (lat_2, lon_2)
1623        points = ImmutableSet([p1, p2, p1])
1624        gsd = Geospatial_data(data_points=list(points),
1625                              points_are_lats_longs=True)
1626
1627        points = gsd.get_data_points(absolute=True)
1628        #Note the order is unknown, due to using sets
1629        # and it changes from windows to linux
1630        try:
[6153]1631            assert num.allclose(points[1][0], 308728.009)
1632            assert num.allclose(points[1][1], 6180432.601)
[6304]1633            assert num.allclose(points[0][0], 222908.705)
[6153]1634            assert num.allclose(points[0][1], 6233785.284)
[6086]1635        except AssertionError:
[6153]1636            assert num.allclose(points[0][0], 308728.009)
1637            assert num.allclose(points[0][1], 6180432.601)
[6304]1638            assert num.allclose(points[1][0], 222908.705)
[6153]1639            assert num.allclose(points[1][1], 6233785.284)
[6304]1640
[6086]1641        self.failUnless(gsd.get_geo_reference().get_zone() == 56,
1642                        'Bad zone error!')
1643        points = gsd.get_data_points(as_lat_long=True)
1644        try:
[6153]1645            assert num.allclose(points[0][0], -34)
1646            assert num.allclose(points[0][1], 150)
[6086]1647        except AssertionError:
[6153]1648            assert num.allclose(points[1][0], -34)
1649            assert num.allclose(points[1][1], 150)
[6086]1650
1651    def test_len(self):
1652        points = [[1.0, 2.1], [3.0, 5.3]]
1653        G = Geospatial_data(points)
[6304]1654        self.failUnless(2 == len(G), 'Len error!')
1655
[6086]1656        points = [[1.0, 2.1]]
1657        G = Geospatial_data(points)
[6304]1658        self.failUnless(1 == len(G), 'Len error!')
[6086]1659
1660        points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
1661        G = Geospatial_data(points)
[6304]1662        self.failUnless(4 == len(G), 'Len error!')
1663
[6086]1664    def test_split(self):
1665        """test if the results from spilt are disjoin sets"""
[6304]1666
1667        # below is a workaround until randint works on cyclones compute nodes
1668        if get_host_name()[8:9] != '0':
[6086]1669            points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
1670                      [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
1671                      [3.0, 1.0], [3.0, 2.0],[3.0, 3.0], [3.0, 4.0], [3.0, 5.0],
1672                      [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
1673                      [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
[6304]1674            attributes = {'depth': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
1675                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
1676                                    21, 22, 23, 24, 25],
1677                          'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
1678                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
1679                                    21, 22, 23, 24, 25]}
[6086]1680            G = Geospatial_data(points, attributes)
[6304]1681
[6086]1682            factor = 0.21
[6304]1683
1684            # will return G1 with 10% of points and G2 with 90%
1685            G1, G2  = G.split(factor, 100)
[6153]1686            assert num.allclose(len(G), len(G1)+len(G2))
1687            assert num.allclose(round(len(G)*factor), len(G1))
[6304]1688
[6086]1689            P = G1.get_data_points(absolute=False)
[6768]1690            expected = [[5.,4.], [4.,1.], [2.,4.], [2.,3.], [1.,4.]]
1691            #            expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0],
1692                                     #            [3.0, 1.0], [2.0, 3.0]]
[6304]1693            msg = 'Expected %s, but\nP=%s' % (str(expected), str(P))
1694            assert num.allclose(P, expected), msg
1695
[6086]1696            A = G1.get_attributes()
[6768]1697            expected = [24, 16, 9, 8, 4]
[6304]1698            msg = 'expected=%s, but A=%s' % (str(expected), str(A))
1699            assert num.allclose(A, expected), msg
1700
[6086]1701    def test_split1(self):
[6304]1702        """test if the results from split are disjoint sets"""
[6086]1703
[6304]1704        # below is a workaround until randint works on cyclones compute nodes
1705        if get_host_name()[8:9] != '0':
[6768]1706            from numpy.random import randint, seed
[6304]1707
[6768]1708            seed((100, 100))
[6304]1709            a_points = randint(0, 999999, (10,2))
[6086]1710            points = a_points.tolist()
[6304]1711
[6086]1712            G = Geospatial_data(points)
[6304]1713
[6086]1714            factor = 0.1
[6304]1715
1716            # will return G1 with 10% of points and G2 with 90%
1717            G1, G2  = G.split(factor, 100)
[6153]1718            assert num.allclose(len(G), len(G1)+len(G2))
1719            assert num.allclose(round(len(G)*factor), len(G1))
[6304]1720
[6086]1721            P = G1.get_data_points(absolute=False)
[6768]1722            expected = [[425835., 137518.]]
[6304]1723            msg = 'expected=%s, but\nP=%s' % (str(expected), str(P))
1724            assert num.allclose(P, expected), msg
[6086]1725
1726    def test_find_optimal_smoothing_parameter(self):
1727        """
[6304]1728        Creates a elevation file representing a hill (sort of) and runs
[6086]1729        find_optimal_smoothing_parameter for 3 different alphas,
[6304]1730
[6086]1731        NOTE the random number seed is provided to control the results
1732        """
[6304]1733
[6086]1734        from cmath import cos
1735
[6304]1736        # below is a workaround until randint works on cyclones compute nodes
[6086]1737        if get_host_name()[8:9]!='0':
[6304]1738            filename = tempfile.mktemp('.csv')
1739            file = open(filename, 'w')
1740            file.write('x,y,elevation \n')
[6086]1741
[6304]1742            for i in range(-5, 6):
1743                for j in range(-5, 6):
1744                    # this equation makes a surface like a circle ripple
[6086]1745                    z = abs(cos(((i*i) + (j*j))*.1)*2)
[6304]1746                    file.write("%s, %s, %s\n" % (i, j, z))
1747
[6086]1748            file.close()
[6304]1749
1750            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
[6086]1751                                                 alpha_list=[0.0001, 0.01, 1],
1752                                                 mesh_file=None,
1753                                                 mesh_resolution=3,
1754                                                 north_boundary=5,
1755                                                 south_boundary=-5,
1756                                                 east_boundary=5,
1757                                                 west_boundary=-5,
1758                                                 plot_name=None,
1759                                                 seed_num=100000,
1760                                                 verbose=False)
[6304]1761
[6086]1762            os.remove(filename)
[6304]1763
[6086]1764            # print value, alpha
[6304]1765            assert(alpha == 0.01)
[6086]1766
1767    def test_find_optimal_smoothing_parameter1(self):
1768        """
[6304]1769        Creates a elevation file representing  a hill (sort of) and
[6086]1770        Then creates a mesh file and passes the mesh file and the elevation
1771        file to find_optimal_smoothing_parameter for 3 different alphas,
[6304]1772
[6086]1773        NOTE the random number seed is provided to control the results
1774        """
[6304]1775
1776        # below is a workaround until randint works on cyclones compute nodes
[6086]1777        if get_host_name()[8:9]!='0':
1778            from cmath import cos
1779            from anuga.pmesh.mesh_interface import create_mesh_from_regions
[6304]1780
1781            filename = tempfile.mktemp('.csv')
1782            file = open(filename, 'w')
1783            file.write('x,y,elevation \n')
1784
1785            for i in range(-5 ,6):
1786                for j in range(-5, 6):
1787                    # this equation makes a surface like a circle ripple
[6086]1788                    z = abs(cos(((i*i) + (j*j))*.1)*2)
[6304]1789                    file.write('%s, %s, %s\n' % (i, j, z))
1790
[6086]1791            file.close()
[6304]1792
1793            poly=[[5,5], [5,-5], [-5,-5], [-5,5]]
1794            internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]]
1795            mesh_filename= tempfile.mktemp('.msh')
1796
[6086]1797            create_mesh_from_regions(poly,
[6304]1798                                     boundary_tags={'back': [2],
1799                                                    'side': [1,3],
1800                                                    'ocean': [0]},
1801                                     maximum_triangle_area=3,
1802                                     interior_regions=internal_poly,
1803                                     filename=mesh_filename,
1804                                     use_cache=False,
1805                                     verbose=False)
1806
1807            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
[6086]1808                                                 alpha_list=[0.0001, 0.01, 1],
1809                                                 mesh_file=mesh_filename,
1810                                                 plot_name=None,
1811                                                 seed_num=174,
1812                                                 verbose=False)
[6304]1813
[6086]1814            os.remove(filename)
1815            os.remove(mesh_filename)
1816
[6768]1817            msg = 'alpha=%s' % str(alpha)
1818            # 0.01 was expected with Numeric.RandomArray RNG
1819            assert alpha==1.0, msg
[6304]1820
[6086]1821    def test_find_optimal_smoothing_parameter2(self):
[6304]1822        '''Tests requirement that mesh file must exist or IOError is thrown
1823
[6086]1824        NOTE the random number seed is provided to control the results
[6304]1825        '''
1826
[6086]1827        from cmath import cos
1828        from anuga.pmesh.mesh_interface import create_mesh_from_regions
[6304]1829
1830        filename = tempfile.mktemp('.csv')
1831        mesh_filename= tempfile.mktemp('.msh')
1832
[6086]1833        try:
[6304]1834            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
[6086]1835                                             alpha_list=[0.0001, 0.01, 1],
1836                                             mesh_file=mesh_filename,
1837                                             plot_name=None,
1838                                             seed_num=174,
1839                                             verbose=False)
1840        except IOError:
1841            pass
1842        else:
[6304]1843            self.fail('Error not thrown error!')
1844
1845################################################################################
1846
[6086]1847if __name__ == "__main__":
1848    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
1849    runner = unittest.TextTestRunner() #verbosity=2)
1850    runner.run(suite)
1851
[6304]1852
Note: See TracBrowser for help on using the repository browser.