/[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 4109 - (hide annotations)
Fri Dec 14 01:08:52 2012 UTC (6 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 6815 byte(s)
Fixed tests.

1 caltinay 3985
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 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 caltinay 3985
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 caltinay 4016 TEST_DATA_ROOT='ref_data'
43 caltinay 3985
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 caltinay 4014 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 caltinay 4016 NC_ORIGIN = [403320.91466610413, 6414860.942530109]
60 caltinay 3985 VMIN=-10000.
61     VMAX=10000
62     NE_V=15
63     ALT=0.
64 caltinay 4016 PAD_X=3
65     PAD_Y=2
66 caltinay 3985
67 caltinay 4061 class TestErMapperData(unittest.TestCase):
68 caltinay 3985 def test_ers_with_padding(self):
69 caltinay 4061 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 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
75 caltinay 4061 dom=domainbuilder.getDomain()
76     g,s=domainbuilder.getGravitySurveys()[0]
77 caltinay 3985
78     outfn=os.path.join(WORKDIR, '_ersdata.csv')
79 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
80 caltinay 3985
81     X0,NP,DX=source.getDataExtents()
82 caltinay 4061 DV=(VMAX-VMIN)/NE_V
83 caltinay 3985
84     # check metadata
85     self.assertEqual(NP, ERS_SIZE, msg="Wrong number of data points")
86 caltinay 4014 # this test only works if gdal is available
87     try:
88     import osgeo.osr
89 jfenwick 4019 for i in range(len(ERS_ORIGIN)):
90 caltinay 4014 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 caltinay 3985
94     # check data
95 caltinay 4005 nx=NP[0]+2*PAD_X
96     ny=NP[1]+2*PAD_Y
97     nz=NE_V
98 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
99 caltinay 3985
100     ref=np.genfromtxt(ERS_REF, delimiter=',', dtype=float)
101 caltinay 4014 g_ref=ref[:,0].reshape((NP[1],NP[0]))
102     s_ref=ref[:,1].reshape((NP[1],NP[0]))
103 caltinay 3985
104     out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=float)
105 caltinay 3987 # recompute nz since ripley might have adjusted number of elements
106     nz=len(out)/(nx*ny)
107 caltinay 3985 g_out=out[:,0].reshape(nz,ny,nx)
108     s_out=out[:,1].reshape(nz,ny,nx)
109 caltinay 4014 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 caltinay 3985
113     self.assertAlmostEqual(np.abs(
114 caltinay 4014 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 caltinay 3985
117     # overwrite data -> should only be padding value left
118 caltinay 4005 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=ERS_NULL
119 caltinay 3985 self.assertAlmostEqual(np.abs(g_out-ERS_NULL).max(), 0.,
120     msg="Wrong values in padding area")
121    
122 caltinay 4061 class TestNetCdfData(unittest.TestCase):
123 caltinay 4014 def test_cdf_with_padding(self):
124 caltinay 4061 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 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
129 caltinay 4061 dom=domainbuilder.getDomain()
130     g,s=domainbuilder.getGravitySurveys()[0]
131 caltinay 4014
132     outfn=os.path.join(WORKDIR, '_ncdata.csv')
133 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
134 caltinay 4014
135     X0,NP,DX=source.getDataExtents()
136 caltinay 4061 DV=(VMAX-VMIN)/NE_V
137 caltinay 4014
138     # check metadata
139     self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
140     # this only works if gdal is available
141 caltinay 4016 try:
142     import osgeo.osr
143 jfenwick 4019 for i in range(len(NC_ORIGIN)):
144 caltinay 4016 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 caltinay 4014
148     # check data
149     nx=NP[0]+2*PAD_X
150     ny=NP[1]+2*PAD_Y
151     nz=NE_V
152 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
153 caltinay 4014
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 caltinay 3985 if __name__ == "__main__":
178     suite = unittest.TestSuite()
179 caltinay 4016 if getMPISizeWorld()==1:
180 caltinay 4061 suite.addTest(unittest.makeSuite(TestErMapperData))
181     if 'NetCdfData' in dir():
182     suite.addTest(unittest.makeSuite(TestNetCdfData))
183 caltinay 4016 else:
184     print("Skipping netCDF data source test since netCDF is not installed")
185 caltinay 4014 else:
186 caltinay 4016 print("Skipping data source tests since MPI size > 1")
187 caltinay 3985 s=unittest.TextTestRunner(verbosity=2).run(suite)
188     if not s.wasSuccessful(): sys.exit(1)
189    

  ViewVC Help
Powered by ViewVC 1.1.26