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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4394 - (hide annotations)
Tue May 7 04:56:59 2013 UTC (6 years, 4 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 caltinay 3985
2     ##############################################################################
3     #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 caltinay 3985 # 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 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 caltinay 3985 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 caltinay 4016 from esys.escript import inf,sup,saveDataCSV,getMPISizeWorld
29 caltinay 3985 from esys.downunder.datasources import *
30 caltinay 4061 from esys.downunder.domainbuilder import DomainBuilder
31 gross 4373 from esys.downunder.coordinates import WGS84ReferenceSystem
32 caltinay 3985
33 caltinay 4364 try:
34     import pyproj
35     haveProj=True
36     except ImportError:
37     haveProj=False
38     mpisize = getMPISizeWorld()
39    
40 caltinay 3985 # this is mainly to avoid warning messages
41 caltinay 4213 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
42 caltinay 3985
43     try:
44     TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
45     except KeyError:
46 caltinay 4016 TEST_DATA_ROOT='ref_data'
47 caltinay 3985
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 caltinay 4014 ERS_SIZE = [20,15]
58 caltinay 4270 ERS_ORIGIN = [309241.0, 6318655.0]
59 caltinay 4014 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 caltinay 4016 NC_ORIGIN = [403320.91466610413, 6414860.942530109]
64 gross 4373 NC_ORIGIN_WGS84 = [115.97200299999747, -32.399100000001042]
65 caltinay 4362 NUMPY_NULL = -123.4
66 caltinay 3985 VMIN=-10000.
67 caltinay 4362 VMAX=10000.
68 caltinay 3985 NE_V=15
69     ALT=0.
70 caltinay 4016 PAD_X=3
71     PAD_Y=2
72 caltinay 3985
73 caltinay 4362 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 caltinay 4364 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
87 caltinay 4362 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 caltinay 4365 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
134 caltinay 4362 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 caltinay 4364 @unittest.skipIf(not haveProj, 'pyproj not available')
183 caltinay 4061 class TestErMapperData(unittest.TestCase):
184 caltinay 4365 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
185 caltinay 3985 def test_ers_with_padding(self):
186 caltinay 4061 source = ErMapperData(DataSource.GRAVITY, headerfile=ERS_DATA,
187 caltinay 4381 altitude=ALT, scale_factor=1e-6)
188 caltinay 4061 domainbuilder=DomainBuilder()
189     domainbuilder.addSource(source)
190     domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
191 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
192 caltinay 4061 dom=domainbuilder.getDomain()
193     g,s=domainbuilder.getGravitySurveys()[0]
194 caltinay 3985
195     outfn=os.path.join(WORKDIR, '_ersdata.csv')
196 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
197 caltinay 3985
198     X0,NP,DX=source.getDataExtents()
199 caltinay 4061 DV=(VMAX-VMIN)/NE_V
200 caltinay 3985
201     # check metadata
202     self.assertEqual(NP, ERS_SIZE, msg="Wrong number of data points")
203 caltinay 4014 # this test only works if gdal is available
204     try:
205     import osgeo.osr
206 jfenwick 4019 for i in range(len(ERS_ORIGIN)):
207 caltinay 4014 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 caltinay 3985
211     # check data
212 caltinay 4005 nx=NP[0]+2*PAD_X
213     ny=NP[1]+2*PAD_Y
214     nz=NE_V
215 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
216 caltinay 3985
217     ref=np.genfromtxt(ERS_REF, delimiter=',', dtype=float)
218 caltinay 4014 g_ref=ref[:,0].reshape((NP[1],NP[0]))
219     s_ref=ref[:,1].reshape((NP[1],NP[0]))
220 caltinay 3985
221     out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=float)
222 caltinay 3987 # recompute nz since ripley might have adjusted number of elements
223     nz=len(out)/(nx*ny)
224 caltinay 3985 g_out=out[:,0].reshape(nz,ny,nx)
225     s_out=out[:,1].reshape(nz,ny,nx)
226 caltinay 4014 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 caltinay 3985
230     self.assertAlmostEqual(np.abs(
231 caltinay 4014 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 caltinay 3985
234     # overwrite data -> should only be padding value left
235 caltinay 4005 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=ERS_NULL
236 caltinay 3985 self.assertAlmostEqual(np.abs(g_out-ERS_NULL).max(), 0.,
237     msg="Wrong values in padding area")
238    
239 caltinay 4364 @unittest.skipIf('NetCdfData' not in dir(), 'netCDF not available')
240 caltinay 4061 class TestNetCdfData(unittest.TestCase):
241 caltinay 4365 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
242 caltinay 4014 def test_cdf_with_padding(self):
243 caltinay 4381 source = NetCdfData(DataSource.GRAVITY, NC_DATA, ALT, scale_factor=1e-6)
244 caltinay 4061 domainbuilder=DomainBuilder()
245     domainbuilder.addSource(source)
246     domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
247 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
248 caltinay 4061 dom=domainbuilder.getDomain()
249     g,s=domainbuilder.getGravitySurveys()[0]
250 caltinay 4014
251     outfn=os.path.join(WORKDIR, '_ncdata.csv')
252 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
253 caltinay 4014
254     X0,NP,DX=source.getDataExtents()
255 caltinay 4061 DV=(VMAX-VMIN)/NE_V
256 caltinay 4014
257     # check metadata
258     self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
259     # this only works if gdal is available
260 gross 4373
261 caltinay 4016 try:
262     import osgeo.osr
263 jfenwick 4019 for i in range(len(NC_ORIGIN)):
264 caltinay 4016 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 caltinay 4014
268     # check data
269     nx=NP[0]+2*PAD_X
270     ny=NP[1]+2*PAD_Y
271     nz=NE_V
272 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
273 gross 4373
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 caltinay 4014
278 gross 4373 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 caltinay 4381
298 gross 4373 @unittest.skipIf(mpisize>1, "more than 1 MPI rank")
299     def test_cdf_with_padding_ellipsoid(self):
300     ref=WGS84ReferenceSystem()
301 caltinay 4394
302 caltinay 4381 source = NetCdfData(DataSource.GRAVITY, NC_DATA, ALT,
303     reference_system=ref, scale_factor=1e-6)
304 gross 4373 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 caltinay 4394 self.assertAlmostEqual(X0[i], NC_ORIGIN_WGS84[i], msg="Data origin wrong")
322 gross 4373
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 caltinay 4014 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 gross 4373
335 caltinay 4014 # 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 gross 4373
353 caltinay 3985 if __name__ == "__main__":
354     suite = unittest.TestSuite()
355 caltinay 4364 suite.addTest(unittest.makeSuite(TestNumpyData))
356     suite.addTest(unittest.makeSuite(TestErMapperData))
357     suite.addTest(unittest.makeSuite(TestNetCdfData))
358 caltinay 3985 s=unittest.TextTestRunner(verbosity=2).run(suite)
359     if not s.wasSuccessful(): sys.exit(1)
360    

  ViewVC Help
Powered by ViewVC 1.1.26