/[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 4394 - (show annotations)
Tue May 7 04:56:59 2013 UTC (6 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 13831 byte(s)
Fixed rounding with lat/lon coordinates and removed stray print statements.

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

  ViewVC Help
Powered by ViewVC 1.1.26