/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (19 months, 2 weeks ago) by jfenwick
File MIME type: text/x-python
File size: 17056 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26