/[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 4213 - (show annotations)
Tue Feb 19 01:16:29 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 6735 byte(s)
Some cleanup and more consistent logging.

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
32 # this is mainly to avoid warning messages
33 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
34
35 try:
36 TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
37 except KeyError:
38 TEST_DATA_ROOT='ref_data'
39
40 try:
41 WORKDIR=os.environ['DOWNUNDER_WORKDIR']
42 except KeyError:
43 WORKDIR='.'
44
45
46 ERS_DATA = os.path.join(TEST_DATA_ROOT, 'ermapper_test.ers')
47 ERS_REF = os.path.join(TEST_DATA_ROOT, 'ermapper_test.csv')
48 ERS_NULL = -99999 * 1e-6
49 ERS_SIZE = [20,15]
50 ERS_ORIGIN = [309097.0, 6319002.0]
51 NC_DATA = os.path.join(TEST_DATA_ROOT, 'netcdf_test.nc')
52 NC_REF = os.path.join(TEST_DATA_ROOT, 'netcdf_test.csv')
53 NC_NULL = 0.
54 NC_SIZE = [20,15]
55 NC_ORIGIN = [403320.91466610413, 6414860.942530109]
56 VMIN=-10000.
57 VMAX=10000
58 NE_V=15
59 ALT=0.
60 PAD_X=3
61 PAD_Y=2
62
63 class TestErMapperData(unittest.TestCase):
64 def test_ers_with_padding(self):
65 source = ErMapperData(DataSource.GRAVITY, headerfile=ERS_DATA,
66 altitude=ALT)
67 domainbuilder=DomainBuilder()
68 domainbuilder.addSource(source)
69 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
70 domainbuilder.setElementPadding(PAD_X,PAD_Y)
71 dom=domainbuilder.getDomain()
72 g,s=domainbuilder.getGravitySurveys()[0]
73
74 outfn=os.path.join(WORKDIR, '_ersdata.csv')
75 saveDataCSV(outfn, g=g, s=s)
76
77 X0,NP,DX=source.getDataExtents()
78 DV=(VMAX-VMIN)/NE_V
79
80 # check metadata
81 self.assertEqual(NP, ERS_SIZE, msg="Wrong number of data points")
82 # this test only works if gdal is available
83 try:
84 import osgeo.osr
85 for i in range(len(ERS_ORIGIN)):
86 self.assertAlmostEqual(X0[i], ERS_ORIGIN[i], msg="Data origin wrong")
87 except ImportError:
88 print("Skipping test of data origin since gdal is not installed.")
89
90 # check data
91 nx=NP[0]+2*PAD_X
92 ny=NP[1]+2*PAD_Y
93 nz=NE_V
94 z_data=int(np.round((ALT-VMIN)/DV)-1)
95
96 ref=np.genfromtxt(ERS_REF, delimiter=',', dtype=float)
97 g_ref=ref[:,0].reshape((NP[1],NP[0]))
98 s_ref=ref[:,1].reshape((NP[1],NP[0]))
99
100 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=float)
101 # recompute nz since ripley might have adjusted number of elements
102 nz=len(out)/(nx*ny)
103 g_out=out[:,0].reshape(nz,ny,nx)
104 s_out=out[:,1].reshape(nz,ny,nx)
105 self.assertAlmostEqual(np.abs(
106 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-g_ref).max(),
107 0., msg="Difference in gravity data area")
108
109 self.assertAlmostEqual(np.abs(
110 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-s_ref).max(),
111 0., msg="Difference in error data area")
112
113 # overwrite data -> should only be padding value left
114 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=ERS_NULL
115 self.assertAlmostEqual(np.abs(g_out-ERS_NULL).max(), 0.,
116 msg="Wrong values in padding area")
117
118 class TestNetCdfData(unittest.TestCase):
119 def test_cdf_with_padding(self):
120 source = NetCdfData(DataSource.GRAVITY, NC_DATA, ALT)
121 domainbuilder=DomainBuilder()
122 domainbuilder.addSource(source)
123 domainbuilder.setVerticalExtents(depth=-VMIN, air_layer=VMAX, num_cells=NE_V)
124 domainbuilder.setElementPadding(PAD_X,PAD_Y)
125 dom=domainbuilder.getDomain()
126 g,s=domainbuilder.getGravitySurveys()[0]
127
128 outfn=os.path.join(WORKDIR, '_ncdata.csv')
129 saveDataCSV(outfn, g=g, s=s)
130
131 X0,NP,DX=source.getDataExtents()
132 DV=(VMAX-VMIN)/NE_V
133
134 # check metadata
135 self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
136 # this only works if gdal is available
137 try:
138 import osgeo.osr
139 for i in range(len(NC_ORIGIN)):
140 self.assertAlmostEqual(X0[i], NC_ORIGIN[i], msg="Data origin wrong")
141 except ImportError:
142 print("Skipping test of data origin since gdal is not installed.")
143
144 # check data
145 nx=NP[0]+2*PAD_X
146 ny=NP[1]+2*PAD_Y
147 nz=NE_V
148 z_data=int(np.round((ALT-VMIN)/DV)-1)
149
150 ref=np.genfromtxt(NC_REF, delimiter=',', dtype=float)
151 g_ref=ref[:,0].reshape((NP[1],NP[0]))
152 s_ref=ref[:,1].reshape((NP[1],NP[0]))
153
154 out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=float)
155 # recompute nz since ripley might have adjusted number of elements
156 nz=len(out)/(nx*ny)
157 g_out=out[:,0].reshape(nz,ny,nx)
158 s_out=out[:,1].reshape(nz,ny,nx)
159
160 self.assertAlmostEqual(np.abs(
161 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-g_ref).max(),
162 0., msg="Difference in gravity data area")
163
164 self.assertAlmostEqual(np.abs(
165 s_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]-s_ref).max(),
166 0., msg="Difference in error data area")
167
168 # overwrite data -> should only be padding value left
169 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=NC_NULL
170 self.assertAlmostEqual(np.abs(g_out-NC_NULL).max(), 0.,
171 msg="Wrong values in padding area")
172
173 if __name__ == "__main__":
174 suite = unittest.TestSuite()
175 if getMPISizeWorld()==1:
176 suite.addTest(unittest.makeSuite(TestErMapperData))
177 if 'NetCdfData' in dir():
178 suite.addTest(unittest.makeSuite(TestNetCdfData))
179 else:
180 print("Skipping netCDF data source test since netCDF is not installed")
181 else:
182 print("Skipping data source tests since MPI size > 1")
183 s=unittest.TextTestRunner(verbosity=2).run(suite)
184 if not s.wasSuccessful(): sys.exit(1)
185

  ViewVC Help
Powered by ViewVC 1.1.26