/[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 4984 - (show annotations)
Mon Jun 2 02:50:34 2014 UTC (5 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 16679 byte(s)
revamping testrunners, now uses automated discovery and allows running specific tests without modifying files (see escriptcore/py_src/testing.py for more info/examples)

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

  ViewVC Help
Powered by ViewVC 1.1.26