/[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 7005 - (show annotations)
Sat Jul 11 00:20:49 2020 UTC (12 months, 2 weeks ago) by aellery
File MIME type: text/x-python
File size: 17209 byte(s)
Fixed a number of broken tests.


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

  ViewVC Help
Powered by ViewVC 1.1.26