source: branches/inundation-numpy-branch/geospatial_data/test_geospatial_data.py @ 6689

Last change on this file since 6689 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 8.9 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5import os
6from numpy import zeros, array, allclose, concatenate
7from math import sqrt, pi
8
9from anuga.geospatial_data.geospatial_data import *
10
11from coordinate_transforms.geo_reference import Geo_reference
12
13class Test_Geospatial_data(unittest.TestCase):
14    def setUp(self):
15        pass
16
17    def tearDown(self):
18        pass
19
20
21    def test_0(self):
22        """Basic points
23        """
24        from coordinate_transforms.geo_reference import Geo_reference
25       
26        points = [[1.0, 2.1], [3.0, 5.3]]
27        G = Geospatial_data(points)
28
29        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
30
31        #Check defaults
32        assert G.attributes is None
33       
34        assert G.geo_reference.zone == Geo_reference().zone
35        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
36        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
37       
38
39    def test_1(self):
40        points = [[1.0, 2.1], [3.0, 5.3]]
41        attributes = [2, 4]
42        G = Geospatial_data(points, attributes)       
43
44        assert G.attributes.keys()[0] == 'attribute'
45        assert allclose(G.attributes.values()[0], [2, 4])
46       
47
48    def test_2(self):
49        from coordinate_transforms.geo_reference import Geo_reference
50        points = [[1.0, 2.1], [3.0, 5.3]]
51        attributes = [2, 4]
52        G = Geospatial_data(points, attributes,
53                            geo_reference = Geo_reference(56, 100, 200))
54
55        assert G.geo_reference.zone == 56
56        assert G.geo_reference.xllcorner == 100
57        assert G.geo_reference.yllcorner == 200
58
59
60    def test_get_attributes_1(self):
61        from coordinate_transforms.geo_reference import Geo_reference
62        points = [[1.0, 2.1], [3.0, 5.3]]
63        attributes = [2, 4]
64        G = Geospatial_data(points, attributes,
65                            geo_reference = Geo_reference(56, 100, 200))
66
67
68        P = G.get_data_points(absolute = False)
69        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
70
71        P = G.get_data_points(absolute = True)
72        assert allclose(P, [[101.0, 202.1], [103.0, 205.3]])       
73
74        V = G.get_attributes() #Simply get them
75        assert allclose(V, [2, 4])
76
77        V = G.get_attributes('attribute') #Get by name
78        assert allclose(V, [2, 4])
79
80    def test_get_attributes_2(self):
81        """Multiple attributes
82        """
83       
84        from coordinate_transforms.geo_reference import Geo_reference
85        points = [[1.0, 2.1], [3.0, 5.3]]
86        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
87        G = Geospatial_data(points, attributes,
88                            geo_reference = Geo_reference(56, 100, 200),
89                            default_attribute_name = 'a1')
90
91
92        P = G.get_data_points(absolute = False)
93        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
94       
95        V = G.get_attributes() #Get default attribute
96        assert allclose(V, [2, 4])
97
98        V = G.get_attributes('a0') #Get by name
99        assert allclose(V, [0, 0])
100
101        V = G.get_attributes('a1') #Get by name
102        assert allclose(V, [2, 4])
103
104        V = G.get_attributes('a2') #Get by name
105        assert allclose(V, [79.4, -7])
106
107        try:
108            V = G.get_attributes('hdnoatedu') #Invalid
109        except AssertionError:
110            pass
111        else:
112            raise 'Should have raised exception' 
113
114
115    def test_conversions_to_points_dict(self):
116        """test conversions to points_dict
117        """
118       
119        from coordinate_transforms.geo_reference import Geo_reference
120        points = [[1.0, 2.1], [3.0, 5.3]]
121        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
122        G = Geospatial_data(points, attributes,
123                            geo_reference = Geo_reference(56, 100, 200),
124                            default_attribute_name = 'a1')
125
126
127        points_dict = geospatial_data2points_dictionary(G)
128
129        assert points_dict.has_key('pointlist')
130        assert points_dict.has_key('attributelist')       
131        assert points_dict.has_key('geo_reference')
132
133        assert allclose( points_dict['pointlist'], points )
134
135        A = points_dict['attributelist']
136        assert A.has_key('a0')
137        assert A.has_key('a1')
138        assert A.has_key('a2')       
139
140        assert allclose( A['a0'], [0, 0] )
141        assert allclose( A['a1'], [2, 4] )       
142        assert allclose( A['a2'], [79.4, -7] )
143
144
145        geo = points_dict['geo_reference']
146        assert geo is G.geo_reference
147
148
149    def test_conversions_from_points_dict(self):
150        """test conversions from points_dict
151        """
152
153        from coordinate_transforms.geo_reference import Geo_reference
154       
155        points = [[1.0, 2.1], [3.0, 5.3]]
156        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
157
158        points_dict = {}
159        points_dict['pointlist'] = points
160        points_dict['attributelist'] = attributes
161        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
162       
163
164        G = points_dictionary2geospatial_data(points_dict)
165
166        P = G.get_data_points(absolute = False)
167        assert allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
168       
169        #V = G.get_attribute_values() #Get default attribute
170        #assert allclose(V, [2, 4])
171
172        V = G.get_attributes('a0') #Get by name
173        assert allclose(V, [0, 0])
174
175        V = G.get_attributes('a1') #Get by name
176        assert allclose(V, [2, 4])
177
178        V = G.get_attributes('a2') #Get by name
179        assert allclose(V, [79.4, -7])
180
181    def test_add(self):
182        """ test the addition of two geospatical objects
183            no geo_reference see next test
184        """
185        points = [[1.0, 2.1], [3.0, 5.3]]
186        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
187        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
188        G1 = Geospatial_data(points, attributes)       
189        G2 = Geospatial_data(points, attributes1) 
190       
191#        print 'G1 depth=', G1.get_attributes('depth')
192#        print 'G1 attrib keys=', G1.attributes.keys()
193        g3 = geospatial_data2points_dictionary(G1)
194#        print 'g3=', g3
195       
196        G = G1 + G2
197
198#        g = geospatial_data2points_dictionary(G)
199#        print 'G=', g
200#        print 'G points =', G.get_data_points()
201#        print 'G attrib keys =', G.attributes.keys()
202#        print 'G test =', G.get_attributes('elevation')
203        assert G.attributes.has_key('depth')
204        assert G.attributes.has_key('elevation')
205        assert allclose(G.attributes['depth'], [2, 4, 2, 4])
206        assert allclose(G.attributes['elevation'], [6.1, 5, 2.5, 1])
207        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
208                                              [1.0, 2.1], [3.0, 5.3]])
209       
210    def test_add1 (self):
211        """
212        difference in Geo_reference resolved
213        """
214        points = [[1.0, 2.1], [3.0, 5.3]]
215        points1 = [[5.0, 6.1], [6.0, 3.3]]
216        attributes = [2, 4]
217        attributes1 = [5, 76]
218        geo_ref = Geo_reference(55, 1.0, 2.0)
219        geo_ref1 = Geo_reference(zone = 55, xllcorner = 0.1,
220                                yllcorner = 3.0, datum = 'wgs84',
221                                projection = 'UTM', units = 'm')
222                               
223        G1 = Geospatial_data(points, attributes, geo_ref)
224        G2 = Geospatial_data(points1, attributes1, geo_ref1)
225
226        G = G1 + G2
227       
228        assert allclose(G.get_geo_reference().get_xllcorner(), 0.1)
229        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
230#        print 'G data points=', G.get_data_points()
231        assert allclose(G.get_data_points(), [[2.0, 4.1], [4.0, 7.3], [5.1, 9.1], [6.1, 6.3]])                             
232
233    def xtest_create_from_file(self):
234        """Check that object can be created from a points file
235        """
236
237        from load_mesh.loadASCII import export_points_file
238       
239
240        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
241        attributes = [2, 4, 5, 76]
242
243        # Use old pointsdict format
244        pointsdict = {'pointlist': points,
245                      'attributelist': {'att1': attributes,
246                                        'att2': array(attributes) + 1}} 
247       
248        # Create points as an xya file
249        FN = 'test_points.xya'
250        export_points_file(FN, pointsdict)
251
252
253        #Create object from file
254        G = Geospatial_data(filename = FN)
255
256        assert allclose(G.get_data_points(), points)
257        assert allclose(G.get_attributes('att1'), attributes)
258        assert allclose(G.get_attributes('att2'), array(attributes) + 1)
259       
260
261       
262        os.remove(FN)
263       
264
265if __name__ == "__main__":
266    suite = unittest.makeSuite(Test_Geospatial_data,'test')
267    runner = unittest.TextTestRunner()
268    runner.run(suite)
269
270   
Note: See TracBrowser for help on using the repository browser.