/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Contents of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4060 - (show annotations)
Fri Nov 9 03:50:19 2012 UTC (8 years ago) by caltinay
File MIME type: text/x-python
File size: 9393 byte(s)
-Overhaul of inversion data sources. Domain is now generated and managed by separate class.
-Removed UBCDataSource which was only used for testing.
-Updated example scripts.

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 """Domain construction from survey data for inversions"""
17
18 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = ['DomainBuilder']
26
27 import logging
28 import numpy as np
29 from esys.escript.util import *
30 import esys.escript.unitsSI as U
31 from esys.ripley import Brick, Rectangle
32 from .datasources import DataSource
33
34 class DomainBuilder(object):
35 """
36 This class is responsible for constructing an escript Domain object with
37 suitable extents and resolution for survey data (`DataSource` objects)
38 that is added to it.
39 """
40
41 def __init__(self, dim=3):
42 """
43 Constructor.
44
45 :param dim: Dimensionality (2 or 3) of the target domain.
46 This has implications for the survey data than can be
47 added. By default a 3D domain is created.
48 :type dim: ``int``
49 """
50 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
51 if dim not in (2,3):
52 raise ValueError("Number of dimensions must be 2 or 3")
53 self._domain=None
54 self._dim=dim
55 self._gravity_surveys=[]
56 self._magnetic_surveys=[]
57 self._sources=[]
58 self.setPadding()
59 self.setVerticalExtents()
60
61 def addSource(self, source):
62 """
63 Adds a survey data provider to the domain builder.
64
65 :param source: The data source to be added.
66 :type source: `DataSource`
67 """
68 if not isinstance(source, DataSource):
69 raise TypeError("source is not a DataSource")
70 self._sources.append(source)
71
72 def setPadding(self, pad_x=10, pad_y=10):
73 """
74 Sets the amount of padding around the dataset. If `pad_x`/`pad_y`
75 is >=1 the value is treated as number of elements to be added to the
76 domain (per side).
77 If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the
78 dataset size. For example, calling ``setPadding(3, 0.1)`` to a data
79 source with size 10x20 will result in the padded data set size
80 16x24 (10+2*3, 20*(1+2*0.1))
81
82 :param pad_x: Padding per side in x direction (default: 10 elements)
83 :type pad_x: ``int`` or ``float``
84 :param pad_y: Padding per side in y direction (default: 10 elements).
85 This value is only used for 3-dimensional datasets
86 :type pad_y: ``int`` or ``float``
87 """
88 self._padding=[pad_x,pad_y]
89
90 def getGravitySurveys(self):
91 """
92 Returns a list of `Data` objects for all gravity surveys available to
93 this domain builder.
94
95 :return: List of gravity surveys which are tuples (anomaly,error).
96 :rtype: ``list``
97 """
98 if len(self._gravity_surveys)==0:
99 for src in self._sources:
100 if src.getDataType()==DataSource.GRAVITY:
101 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
102 self._gravity_surveys.append(survey)
103 return self._gravity_surveys
104
105 def getMagneticSurveys(self):
106 """
107 Returns a list of `Data` objects for all magnetic surveys available to
108 this domain builder.
109
110 :return: List of magnetic surveys which are tuples (anomaly,error).
111 :rtype: ``list``
112 """
113 if len(self._magnetic_surveys)==0:
114 for src in self._sources:
115 if src.getDataType()==DataSource.MAGNETIC:
116 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
117 self._magnetic_surveys.append(survey)
118 return self._magnetic_surveys
119
120 def getBackgroundMagneticField(self):
121 #FIXME:
122 latitude=-28.5
123 theta = (90-latitude)/180.*np.pi
124 B_0 = U.Mu_0*U.Magnetic_Dipole_Moment_Earth / (4*np.pi*U.R_Earth**3)
125 B_theta = B_0*sin(theta)
126 B_r = 2*B_0*cos(theta)
127 if self._dim < 3:
128 return np.array([0., -B_r])
129 else:
130 return np.array([-B_theta, 0., -B_r])
131
132 def getSetDensityMask(self):
133 return wherePositive(self.getDomain().getX()[self._dim-1])
134
135 def getSetSusceptibilityMask(self):
136 return wherePositive(self.getDomain().getX()[self._dim-1])
137
138 def getDomain(self):
139 """
140 Returns a domain that spans the data area plus padding.
141 The domain is created the first time this method is called, subsequent
142 calls return the same domain so anything that affects the domain
143 (such as padding) needs to be set beforehand.
144
145 :return: The escript domain for this data source.
146 :rtype: `esys.escript.Domain`
147 """
148 if self._domain is None:
149 self._domain=self.__createDomain()
150 return self._domain
151
152 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
153 """
154 This method sets the target domain parameters for the vertical
155 dimension.
156
157 :param depth: Depth in meters of the domain.
158 :type depth: ``float``
159 :param air_layer: Depth of the layer above sea level in meters
160 :type air_layer: ``float``
161 :param num_cells: Number of domain elements for the entire vertical
162 dimension
163 :type num_cells: ``int``
164 """
165 self._v_depth=depth
166 self._v_air_layer=air_layer
167 self._v_num_cells=num_cells
168
169 def __getTotalExtentsWithPadding(self):
170 """
171 Helper method that computes origin and number of elements
172 after adding padding to the bounding box of all available survey data.
173 """
174 X0, NX, DX = self.__getTotalExtents()
175 DIM=len(X0)
176 frac=[]
177 # padding is applied to each side so multiply by 2 to get the total
178 # amount of padding per dimension
179 for i in range(DIM):
180 if self._padding[i]>=0 and self._padding[i]<1:
181 frac.append(2.*self._padding[i])
182 else:
183 frac.append(2.*self._padding[i]/float(NX[i]))
184
185 # calculate new number of elements
186 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
187 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
188 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
189 return X0_padded, NX_padded, DX
190
191 def __getTotalExtents(self):
192 """
193 Helper method that computes the origin, number of elements and
194 minimal element spacing taking into account all available survey data.
195 """
196 if len(self._sources)==0:
197 raise ValueError("No data")
198 X0, NE, DX = self._sources[0].getDataExtents()
199 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
200
201 for src in self._sources[1:]:
202 d_x0, d_ne, d_dx = src.getDataExtents()
203 for i in range(len(d_x0)):
204 X0[i]=min(X0[i], d_x0[i])
205 for i in range(len(d_dx)):
206 DX[i]=min(DX[i], d_dx[i])
207 for i in range(len(d_ne)):
208 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
209 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
210 return X0, NE, DX
211
212 def __createDomain(self):
213 """
214 Creates and returns an escript domain that spans the entire area of
215 available data plus a padding zone. This method is called only once
216 the first time `getDomain()` is invoked.
217
218 :return: The escript domain
219 :rtype: `esys.escript.Domain`
220 """
221 X0, NX, DX = self.__getTotalExtentsWithPadding()
222
223 # number of domain elements
224 NE = NX + [self._v_num_cells]
225
226 # origin of domain
227 origin = X0 + [-self._v_depth]
228 self._dom_origin = [np.floor(oi) for oi in origin]
229
230 # cell size / point spacing
231 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
232 self._spacing = [float(np.floor(si)) for si in spacing]
233
234 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
235 if self._dim==3:
236 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
237 else:
238 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
239
240 # ripley may internally adjust NE and length, so recompute
241 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
242 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
243
244 self.logger.debug("Domain size: "+str(self._dom_NE))
245 self.logger.debug(" length: "+str(self._dom_len))
246 self.logger.debug(" origin: "+str(self._dom_origin))
247 return dom
248

  ViewVC Help
Powered by ViewVC 1.1.26