/[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 4061 - (show annotations)
Fri Nov 9 04:05:28 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 6801 byte(s)
Updated data source test cases.

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

  ViewVC Help
Powered by ViewVC 1.1.26