/[escript]/trunk/downunder/test/python/run_datasources.py
ViewVC logotype

Contents of /trunk/downunder/test/python/run_datasources.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4938 - (show annotations)
Wed May 14 01:13:23 2014 UTC (5 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 16942 byte(s)
Modify unit tests to read their classes from
esys.escriptcore.utestselect

Change the line in that file to switch between unittest and unittest2


1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 import logging
25 import numpy as np
26 import os
27 import sys
28 import esys.escriptcore.utestselect as unittest
29 from esys.escript import inf,sup,saveDataCSV,getMPISizeWorld
30 from esys.downunder.datasources import *
31 from esys.downunder.domainbuilder import DomainBuilder
32 from esys.downunder.coordinates import WGS84ReferenceSystem
33
34 try:
35 import pyproj
36 haveProj=True
37 except ImportError:
38 haveProj=False
39 mpisize = getMPISizeWorld()
40
41 # this is mainly to avoid warning messages
42 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
43
44 try:
45 TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
46 except KeyError:
47 TEST_DATA_ROOT='ref_data'
48
49 try:
50 WORKDIR=os.environ['DOWNUNDER_WORKDIR']
51 except KeyError:
52 WORKDIR='.'
53
54
55 ERS32_DATA = os.path.join(TEST_DATA_ROOT, 'ermapper32_test.ers')
56 ERS64_DATA = os.path.join(TEST_DATA_ROOT, 'ermapper64_test.ers')
57 ERS_REF = os.path.join(TEST_DATA_ROOT, 'ermapper_test.csv')
58 ERS_NULL = -99999 * 1e-6
59 ERS_SIZE = [20,15]
60 ERS_ORIGIN = [309241.0, 6318655.0]
61 NC_DATA = os.path.join(TEST_DATA_ROOT, 'netcdf_test.nc')
62 NC_REF = os.path.join(TEST_DATA_ROOT, 'netcdf_test.csv')
63 NC_NULL = 0.
64 NC_SIZE = [20,15]
65 NC_ORIGIN = [403320.91466610413, 6414860.942530109]
66 NC_ORIGIN_WGS84 = [115.97200299999747, -32.399100000001042]
67 NUMPY_NULL = -123.4
68 VMIN=-10000.
69 VMAX=10000.
70 NE_V=15
71 ALT=0.
72 PAD_X=3
73 PAD_Y=2
74
75 class TestNumpyData(unittest.TestCase):
76 def test_numpy_argument_check(self):
77 # invalid data type
78 self.assertRaises(ValueError, NumpyData, '_mydatatype_', [1,2])
79 # invalid shape of data
80 self.assertRaises(ValueError, NumpyData, DataSource.GRAVITY, 42)
81 # invalid shape of data
82 self.assertRaises(ValueError, NumpyData, DataSource.GRAVITY, np.zeros((2,2,2,2)))
83 # invalid shape of error
84 self.assertRaises(ValueError, NumpyData, DataSource.GRAVITY, [1,2], [1,2,3])
85 # invalid shape of length
86 self.assertRaises(ValueError, NumpyData, DataSource.GRAVITY, [1,2], [1,2], [2,3,2])
87
88 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
89 def test_numpy_data_1d(self):
90 DIM=1
91 testdata = np.arange(20)
92 error = 1.*np.ones(testdata.shape)
93 source = NumpyData(DataSource.GRAVITY, testdata, null_value=NUMPY_NULL)
94 X0,NP,DX=source.getDataExtents()
95 for i in range(DIM):
96 self.assertAlmostEqual(X0[i], 0., msg="Data origin wrong")
97 self.assertEqual(NP[i], testdata.shape[DIM-i-1], msg="Wrong number of data points")
98 self.assertAlmostEqual(DX[i], 1000./testdata.shape[DIM-i-1], msg="Wrong cell size")
99
100 domainbuilder=DomainBuilder(dim=2)
101 domainbuilder.addSource(source)
102 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
103 domainbuilder.setElementPadding(PAD_X)
104 dom=domainbuilder.getDomain()
105 g,s=domainbuilder.getGravitySurveys()[0]
106
107 outfn=os.path.join(WORKDIR, '_npdata1d.csv')
108 saveDataCSV(outfn, g=g, s=s)
109
110 DV=(VMAX-VMIN)/NE_V
111
112 # check data
113 nx=NP[0]+2*PAD_X
114 nz=NE_V
115 z_data=int(np.round((ALT-VMIN)/DV)-1)
116
117 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=np.float64)
118 # recompute nz since ripley might have adjusted number of elements
119 nz=len(out)/nx
120 g_out=out[:,0].reshape(nz,nx)
121 s_out=out[:,1].reshape(nz,nx)
122 self.assertAlmostEqual(np.abs(
123 g_out[z_data, PAD_X:PAD_X+NP[0]]-testdata).max(),
124 0., msg="Difference in gravity data area")
125
126 self.assertAlmostEqual(np.abs(
127 s_out[z_data, PAD_X:PAD_X+NP[0]]-error).max(),
128 0., msg="Difference in error data area")
129
130 # overwrite data -> should only be padding value left
131 g_out[z_data, PAD_X:PAD_X+NP[0]]=NUMPY_NULL
132 self.assertAlmostEqual(np.abs(g_out-NUMPY_NULL).max(), 0.,
133 msg="Wrong values in padding area")
134
135 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
136 def test_numpy_data_2d(self):
137 DIM=2
138 testdata = np.arange(20*21).reshape(20,21)
139 error = 1.*np.ones(testdata.shape)
140 source = NumpyData(DataSource.GRAVITY, testdata, null_value=NUMPY_NULL)
141 X0,NP,DX=source.getDataExtents()
142 for i in range(DIM):
143 self.assertAlmostEqual(X0[i], 0., msg="Data origin wrong")
144 self.assertEqual(NP[i], testdata.shape[DIM-i-1], msg="Wrong number of data points")
145 self.assertAlmostEqual(DX[i], 1000./testdata.shape[DIM-i-1], msg="Wrong cell size")
146
147 domainbuilder=DomainBuilder(dim=3)
148 domainbuilder.addSource(source)
149 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
150 domainbuilder.setElementPadding(PAD_X, PAD_Y)
151 dom=domainbuilder.getDomain()
152 g,s=domainbuilder.getGravitySurveys()[0]
153
154 outfn=os.path.join(WORKDIR, '_npdata2d.csv')
155 saveDataCSV(outfn, g=g, s=s)
156
157 DV=(VMAX-VMIN)/NE_V
158
159 # check data
160 nx=NP[0]+2*PAD_X
161 ny=NP[1]+2*PAD_Y
162 nz=NE_V
163 z_data=int(np.round((ALT-VMIN)/DV)-1)
164
165 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=np.float64)
166 # recompute nz since ripley might have adjusted number of elements
167 nz=len(out)/(nx*ny)
168 g_out=out[:,0].reshape(nz,ny,nx)
169 s_out=out[:,1].reshape(nz,ny,nx)
170 self.assertAlmostEqual(np.abs(
171 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-testdata).max(),
172 0., msg="Difference in gravity data area")
173
174 self.assertAlmostEqual(np.abs(
175 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-error).max(),
176 0., msg="Difference in error data area")
177
178 # overwrite data -> should only be padding value left
179 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=NUMPY_NULL
180 self.assertAlmostEqual(np.abs(g_out-NUMPY_NULL).max(), 0.,
181 msg="Wrong values in padding area")
182
183
184 @unittest.skipIf(not haveProj, 'pyproj not available')
185 class TestErMapperData(unittest.TestCase):
186 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
187 def test_ers32_with_padding(self):
188 self._ers_tester(ERS32_DATA)
189
190 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
191 def test_ers64_with_padding(self):
192 self._ers_tester(ERS64_DATA)
193
194 def _ers_tester(self, filename):
195 source = ErMapperData(DataSource.GRAVITY, headerfile=filename,
196 altitude=ALT, scale_factor=1e-6)
197 domainbuilder=DomainBuilder()
198 domainbuilder.addSource(source)
199 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
200 domainbuilder.setElementPadding(PAD_X,PAD_Y)
201 dom=domainbuilder.getDomain()
202 g,s=domainbuilder.getGravitySurveys()[0]
203
204 outfn=os.path.join(WORKDIR, '_ersdata.csv')
205 saveDataCSV(outfn, g=g, s=s)
206
207 X0,NP,DX=source.getDataExtents()
208 DV=(VMAX-VMIN)/NE_V
209
210 # check metadata
211 self.assertEqual(NP, ERS_SIZE, msg="Wrong number of data points")
212 # this test only works if gdal is available
213 try:
214 import osgeo.osr
215 for i in range(len(ERS_ORIGIN)):
216 self.assertAlmostEqual(X0[i], ERS_ORIGIN[i], msg="Data origin wrong")
217 except ImportError:
218 print("Skipping test of data origin since gdal is not installed.")
219
220 # check data
221 nx=NP[0]+2*PAD_X
222 ny=NP[1]+2*PAD_Y
223 nz=NE_V
224 z_data=int(np.round((ALT-VMIN)/DV)-1)
225
226 ref=np.genfromtxt(ERS_REF, delimiter=',', dtype=np.float64)
227 g_ref=ref[:,0].reshape((NP[1],NP[0]))
228 s_ref=ref[:,1].reshape((NP[1],NP[0]))
229
230 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=np.float64)
231 # recompute nz since ripley might have adjusted number of elements
232 nz=len(out)/(nx*ny)
233 g_out=out[:,0].reshape(nz,ny,nx)
234 s_out=out[:,1].reshape(nz,ny,nx)
235 self.assertAlmostEqual(np.abs(
236 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-g_ref).max(),
237 0., msg="Difference in gravity data area")
238
239 self.assertAlmostEqual(np.abs(
240 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-s_ref).max(),
241 0., msg="Difference in error data area")
242
243 # overwrite data -> should only be padding value left
244 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=ERS_NULL
245 self.assertAlmostEqual(np.abs(g_out-ERS_NULL).max(), 0.,
246 msg="Wrong values in padding area")
247
248 @unittest.skipIf('NetCdfData' not in dir(), 'netCDF not available')
249 class TestNetCdfData(unittest.TestCase):
250 @unittest.skipIf(not haveProj, 'pyproj not available')
251 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
252 @unittest.skipIf(not haveProj, 'pyproj not available')
253 def test_cdf_with_padding(self):
254 source = NetCdfData(DataSource.GRAVITY, NC_DATA, ALT, scale_factor=1e-6)
255 domainbuilder=DomainBuilder()
256 domainbuilder.addSource(source)
257 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
258 domainbuilder.setElementPadding(PAD_X,PAD_Y)
259 dom=domainbuilder.getDomain()
260 g,s=domainbuilder.getGravitySurveys()[0]
261
262 outfn=os.path.join(WORKDIR, '_ncdata.csv')
263 saveDataCSV(outfn, g=g, s=s)
264
265 X0,NP,DX=source.getDataExtents()
266 DV=(VMAX-VMIN)/NE_V
267
268 # check metadata
269 self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
270 # this only works if gdal is available
271
272 try:
273 import osgeo.osr
274 for i in range(len(NC_ORIGIN)):
275 self.assertAlmostEqual(X0[i], NC_ORIGIN[i], msg="Data origin wrong")
276 except ImportError:
277 print("Skipping test of data origin since gdal is not installed.")
278
279 # check data
280 nx=NP[0]+2*PAD_X
281 ny=NP[1]+2*PAD_Y
282 nz=NE_V
283 z_data=int(np.round((ALT-VMIN)/DV)-1)
284
285 ref=np.genfromtxt(NC_REF, delimiter=',', dtype=np.float64)
286 g_ref=ref[:,0].reshape((NP[1],NP[0]))
287 s_ref=ref[:,1].reshape((NP[1],NP[0]))
288
289 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=np.float64)
290 # recompute nz since ripley might have adjusted number of elements
291
292 nz=len(out)/(nx*ny)
293 g_out=out[:,0].reshape(nz,ny,nx)
294 s_out=out[:,1].reshape(nz,ny,nx)
295
296 self.assertAlmostEqual(np.abs(
297 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-g_ref).max(),
298 0., msg="Difference in gravity data area")
299
300 self.assertAlmostEqual(np.abs(
301 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-s_ref).max(),
302 0., msg="Difference in error data area")
303
304 # overwrite data -> should only be padding value left
305 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=NC_NULL
306 self.assertAlmostEqual(np.abs(g_out-NC_NULL).max(), 0.,
307 msg="Wrong values in padding area")
308
309 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
310 def test_cdf_with_padding_ellipsoid(self):
311 ref=WGS84ReferenceSystem()
312
313 source = NetCdfData(DataSource.GRAVITY, NC_DATA, ALT,
314 reference_system=ref, scale_factor=1e-6)
315 domainbuilder=DomainBuilder(reference_system=ref)
316 domainbuilder.addSource(source)
317 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
318 domainbuilder.setElementPadding(PAD_X,PAD_Y)
319 dom=domainbuilder.getDomain()
320 g,s=domainbuilder.getGravitySurveys()[0]
321
322 outfn=os.path.join(WORKDIR, '_ncdata.csv')
323 saveDataCSV(outfn, g=g, s=s)
324
325 X0,NP,DX=source.getDataExtents()
326 DV=(VMAX-VMIN)/NE_V
327
328 # check metadata
329 self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
330
331 for i in range(len(NC_ORIGIN)):
332 self.assertAlmostEqual(X0[i], NC_ORIGIN_WGS84[i], msg="Data origin wrong")
333
334 # check data
335 nx=NP[0]+2*PAD_X
336 ny=NP[1]+2*PAD_Y
337 nz=NE_V
338 z_data=int(np.round((ALT-VMIN)/DV)-1)
339
340 ref=np.genfromtxt(NC_REF, delimiter=',', dtype=np.float64)
341 g_ref=ref[:,0].reshape((NP[1],NP[0]))
342 s_ref=ref[:,1].reshape((NP[1],NP[0]))
343
344 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=np.float64)
345
346 # recompute nz since ripley might have adjusted number of elements
347 nz=len(out)/(nx*ny)
348 g_out=out[:,0].reshape(nz,ny,nx)
349 s_out=out[:,1].reshape(nz,ny,nx)
350
351 self.assertAlmostEqual(np.abs(
352 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-g_ref).max(),
353 0., msg="Difference in gravity data area")
354
355 self.assertAlmostEqual(np.abs(
356 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-s_ref).max(),
357 0., msg="Difference in error data area")
358
359 # overwrite data -> should only be padding value left
360 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=NC_NULL
361 self.assertAlmostEqual(np.abs(g_out-NC_NULL).max(), 0.,
362 msg="Wrong values in padding area")
363 class TestSeimicSource(unittest.TestCase):
364 def test_seismic_source(self):
365 ss= SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,6.], elevation=7.)
366 self.assertEqual(ss.getLocation(), (1.,2) )
367 self.assertEqual(ss.getFrequency(), 3.)
368 self.assertEqual(ss.getPower(), complex(4.,-4.) )
369 self.assertEqual(ss.getOrientation(), [5.,6.])
370 self.assertEqual(ss.getElevation(), 7.)
371
372 self.assertTrue( ss == SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
373 self.assertFalse( ss == SeismicSource(x=-1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
374 self.assertFalse( ss == SeismicSource(x=1., y=-2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
375 self.assertFalse( ss == SeismicSource(x=1., y=2., omega=-3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
376 self.assertFalse( ss == SeismicSource(x=1., y=2., omega=3., power = complex(-4.,-4.), orientation=[5.,6.],elevation=7.) )
377 self.assertFalse( ss == SeismicSource(x=1., y=2., omega=3., power = complex(4.,4.), orientation=[5.,6.],elevation=7.) )
378 self.assertFalse( ss == SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,-6.],elevation=7.) )
379 self.assertFalse( ss == SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,-6.],elevation=-7.) )
380
381 self.assertFalse( ss != SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
382 self.assertTrue( ss != SeismicSource(x=-1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
383 self.assertTrue( ss != SeismicSource(x=1., y=-2., omega=3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
384 self.assertTrue( ss != SeismicSource(x=1., y=2., omega=-3., power = complex(4.,-4.), orientation=[5.,6.],elevation=7.) )
385 self.assertTrue( ss != SeismicSource(x=1., y=2., omega=3., power = complex(-4.,-4.), orientation=[5.,6.],elevation=7.) )
386 self.assertTrue( ss != SeismicSource(x=1., y=2., omega=3., power = complex(4.,4.), orientation=[5.,6.],elevation=7.) )
387 self.assertTrue( ss != SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,-6.],elevation=7.) )
388 self.assertTrue( ss != SeismicSource(x=1., y=2., omega=3., power = complex(4.,-4.), orientation=[5.,-6.],elevation=-7.) )
389
390 if __name__ == "__main__":
391 suite = unittest.TestSuite()
392 suite.addTest(unittest.makeSuite(TestNumpyData))
393 suite.addTest(unittest.makeSuite(TestErMapperData))
394 suite.addTest(unittest.makeSuite(TestNetCdfData))
395 suite.addTest(unittest.makeSuite(TestSeimicSource))
396 s=unittest.TextTestRunner(verbosity=2).run(suite)
397 if not s.wasSuccessful(): sys.exit(1)
398

  ViewVC Help
Powered by ViewVC 1.1.26