/[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 4270 - (hide annotations)
Mon Mar 4 04:01:28 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 6735 byte(s)
Corrected test for modified handling of ER Mapper origin.

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 caltinay 3985
32     # this is mainly to avoid warning messages
33 caltinay 4213 logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
34 caltinay 3985
35     try:
36     TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
37     except KeyError:
38 caltinay 4016 TEST_DATA_ROOT='ref_data'
39 caltinay 3985
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 caltinay 4014 ERS_SIZE = [20,15]
50 caltinay 4270 ERS_ORIGIN = [309241.0, 6318655.0]
51 caltinay 4014 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 caltinay 4016 NC_ORIGIN = [403320.91466610413, 6414860.942530109]
56 caltinay 3985 VMIN=-10000.
57     VMAX=10000
58     NE_V=15
59     ALT=0.
60 caltinay 4016 PAD_X=3
61     PAD_Y=2
62 caltinay 3985
63 caltinay 4061 class TestErMapperData(unittest.TestCase):
64 caltinay 3985 def test_ers_with_padding(self):
65 caltinay 4061 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 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
71 caltinay 4061 dom=domainbuilder.getDomain()
72     g,s=domainbuilder.getGravitySurveys()[0]
73 caltinay 3985
74     outfn=os.path.join(WORKDIR, '_ersdata.csv')
75 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
76 caltinay 3985
77     X0,NP,DX=source.getDataExtents()
78 caltinay 4061 DV=(VMAX-VMIN)/NE_V
79 caltinay 3985
80     # check metadata
81     self.assertEqual(NP, ERS_SIZE, msg="Wrong number of data points")
82 caltinay 4014 # this test only works if gdal is available
83     try:
84     import osgeo.osr
85 jfenwick 4019 for i in range(len(ERS_ORIGIN)):
86 caltinay 4014 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 caltinay 3985
90     # check data
91 caltinay 4005 nx=NP[0]+2*PAD_X
92     ny=NP[1]+2*PAD_Y
93     nz=NE_V
94 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
95 caltinay 3985
96     ref=np.genfromtxt(ERS_REF, delimiter=',', dtype=float)
97 caltinay 4014 g_ref=ref[:,0].reshape((NP[1],NP[0]))
98     s_ref=ref[:,1].reshape((NP[1],NP[0]))
99 caltinay 3985
100     out=np.genfromtxt(outfn, delimiter=',', skip_header=1, dtype=float)
101 caltinay 3987 # recompute nz since ripley might have adjusted number of elements
102     nz=len(out)/(nx*ny)
103 caltinay 3985 g_out=out[:,0].reshape(nz,ny,nx)
104     s_out=out[:,1].reshape(nz,ny,nx)
105 caltinay 4014 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 caltinay 3985
109     self.assertAlmostEqual(np.abs(
110 caltinay 4014 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 caltinay 3985
113     # overwrite data -> should only be padding value left
114 caltinay 4005 g_out[z_data, PAD_Y:PAD_Y+NP[1], PAD_X:PAD_X+NP[0]]=ERS_NULL
115 caltinay 3985 self.assertAlmostEqual(np.abs(g_out-ERS_NULL).max(), 0.,
116     msg="Wrong values in padding area")
117    
118 caltinay 4061 class TestNetCdfData(unittest.TestCase):
119 caltinay 4014 def test_cdf_with_padding(self):
120 caltinay 4061 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 caltinay 4109 domainbuilder.setElementPadding(PAD_X,PAD_Y)
125 caltinay 4061 dom=domainbuilder.getDomain()
126     g,s=domainbuilder.getGravitySurveys()[0]
127 caltinay 4014
128     outfn=os.path.join(WORKDIR, '_ncdata.csv')
129 caltinay 4061 saveDataCSV(outfn, g=g, s=s)
130 caltinay 4014
131     X0,NP,DX=source.getDataExtents()
132 caltinay 4061 DV=(VMAX-VMIN)/NE_V
133 caltinay 4014
134     # check metadata
135     self.assertEqual(NP, NC_SIZE, msg="Wrong number of data points")
136     # this only works if gdal is available
137 caltinay 4016 try:
138     import osgeo.osr
139 jfenwick 4019 for i in range(len(NC_ORIGIN)):
140 caltinay 4016 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 caltinay 4014
144     # check data
145     nx=NP[0]+2*PAD_X
146     ny=NP[1]+2*PAD_Y
147     nz=NE_V
148 caltinay 4061 z_data=int(np.round((ALT-VMIN)/DV)-1)
149 caltinay 4014
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 caltinay 3985 if __name__ == "__main__":
174     suite = unittest.TestSuite()
175 caltinay 4016 if getMPISizeWorld()==1:
176 caltinay 4061 suite.addTest(unittest.makeSuite(TestErMapperData))
177     if 'NetCdfData' in dir():
178     suite.addTest(unittest.makeSuite(TestNetCdfData))
179 caltinay 4016 else:
180     print("Skipping netCDF data source test since netCDF is not installed")
181 caltinay 4014 else:
182 caltinay 4016 print("Skipping data source tests since MPI size > 1")
183 caltinay 3985 s=unittest.TextTestRunner(verbosity=2).run(suite)
184     if not s.wasSuccessful(): sys.exit(1)
185    

  ViewVC Help
Powered by ViewVC 1.1.26