/[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 4343 - (show annotations)
Wed Mar 27 06:44:36 2013 UTC (6 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 14423 byte(s)
 
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 """Domain construction from survey data for inversions"""
17
18 __copyright__="""Copyright (c) 2003-2013 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 from esys.escript import 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.__sources=[]
56 self.__background_magnetic_field=None
57 self.setPadding()
58 self.setVerticalExtents()
59 self.fixDensityBelow()
60 self.fixSusceptibilityBelow()
61
62 def addSource(self, source):
63 """
64 Adds a survey data provider to the domain builder.
65 An exception is raised if the domain has already been built or if the
66 UTM zone of `source` does not match the UTM zone of sources already
67 added to the domain builder.
68
69 :param source: The data source to be added
70 :type source: `DataSource`
71 """
72 if self.__domain is not None:
73 raise RuntimeError("Invalid call to addSource(). Domain is already built.")
74 if not isinstance(source, DataSource):
75 raise TypeError("source is not a DataSource")
76
77 if len(self.__sources)>0:
78 if self.__sources[0].getUtmZone() != source.getUtmZone():
79 raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")
80
81 self.__sources.append(source)
82
83 def setFractionalPadding(self, pad_x=None, pad_y=None):
84 """
85 Sets the amount of padding around the dataset as a fraction of the
86 dataset side lengths.
87
88 For example, calling ``setFractionalPadding(0.2, 0.1)`` with a data
89 source of size 10x20 will result in the padded data set size
90 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
91
92 :param pad_x: Padding per side in x direction (default: no padding)
93 :type pad_x: ``float``
94 :param pad_y: Padding per side in y direction (default: no padding)
95 :type pad_y: ``float``
96 :note: `pad_y` is ignored for 2-dimensional domains.
97 """
98 if self.__domain is not None:
99 raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
100 if pad_x is not None:
101 if pad_x < 0:
102 raise ValueError("setFractionalPadding: Arguments must be non-negative")
103 if pad_x > 10:
104 raise ValueError("setFractionalPadding: Argument too large")
105 if pad_y is not None:
106 if pad_y < 0:
107 raise ValueError("setFractionalPadding: Arguments must be non-negative")
108 if pad_y > 10:
109 raise ValueError("setFractionalPadding: Argument too large")
110 self._padding = [pad_x,pad_y], 'f'
111
112 def setPadding(self, pad_x=None, pad_y=None):
113 """
114 Sets the amount of padding around the dataset in absolute length units.
115
116 The final domain size will be the length in x (in y) of the dataset
117 plus twice the value of `pad_x` (`pad_y`). The arguments must be
118 non-negative.
119
120 :param pad_x: Padding per side in x direction (default: no padding)
121 :type pad_x: ``float``
122 :param pad_y: Padding per side in y direction (default: no padding)
123 :type pad_y: ``float``
124 :note: `pad_y` is ignored for 2-dimensional domains.
125 """
126 if self.__domain is not None:
127 raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
128 if pad_x is not None:
129 if pad_x < 0:
130 raise ValueError("setPadding: Arguments must be non-negative")
131 if pad_y is not None:
132 if pad_y < 0:
133 raise ValueError("setPadding: Arguments must be non-negative")
134 self._padding = [pad_x,pad_y], 'l'
135
136 def setElementPadding(self, pad_x=None, pad_y=None):
137 """
138 Sets the amount of padding around the dataset in number of elements
139 (cells).
140
141 When the domain is constructed `pad_x` (`pad_y`) elements are added
142 on each side of the x- (y-) dimension. The arguments must be
143 non-negative.
144
145 :param pad_x: Padding per side in x direction (default: no padding)
146 :type pad_x: ``int``
147 :param pad_y: Padding per side in y direction (default: no padding)
148 :type pad_y: ``int``
149 :note: `pad_y` is ignored for 2-dimensional datasets.
150 """
151 if self.__domain is not None:
152 raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
153 if pad_x is not None:
154 if type(pad_x) is not int:
155 raise TypeError("setElementPadding expects integer arguments")
156 if pad_x < 0:
157 raise ValueError("setElementPadding: Arguments must be non-negative")
158 if pad_y is not None:
159 if type(pad_y) is not int:
160 raise TypeError("setElementPadding expects integer arguments")
161 if pad_y < 0:
162 raise ValueError("setElementPadding: Arguments must be non-negative")
163 self._padding = [pad_x,pad_y], 'e'
164
165 def getGravitySurveys(self):
166 """
167 Returns a list of gravity surveys, see `getSurveys` for details.
168 """
169 return self.getSurveys(DataSource.GRAVITY)
170
171 def getMagneticSurveys(self):
172 """
173 Returns a list of magnetic surveys, see `getSurveys` for details.
174 """
175 return self.getSurveys(DataSource.MAGNETIC)
176
177 def fixDensityBelow(self, depth=None):
178 """
179 Defines the depth below which the density anomaly is set to a given value.
180 if no value is given zero is assumed.
181
182 :param depth: depth below which the density is fixed. If not set , no constraint
183 at depth is applied.
184 :type depth: ``float``
185 """
186 self.__fix_density_below=depth
187
188 def fixSusceptibilityBelow(self, depth=None):
189 """
190 Defines the depth below which the susceptibility anomaly is set to a given value.
191 if no value is given zero is assumed.
192
193 :param depth: depth below which the susceptibility is fixed. If not set , no constraint
194 at depth is applied.
195 :type depth: ``float``
196 """
197 self.__fix_susceptibility_below=depth
198
199 def getSurveys(self, datatype):
200 """
201 Returns a list of `Data` objects for all surveys of type `datatype`
202 available to this domain builder.
203
204 :return: List of surveys which are tuples (anomaly,error).
205 :rtype: ``list``
206 """
207 surveys=[]
208 for src in self.__sources:
209 if src.getDataType()==datatype:
210 surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
211 return surveys
212
213 def setBackgroundMagneticFluxDensity(self, B):
214 """
215 Sets the background magnetic flux density B=(B_North, B_East, B_Vertical)
216 """
217 self.__background_magnetic_field=B
218
219 def getBackgroundMagneticFluxDensity(self):
220 """
221 Returns the background magnetic flux density.
222 """
223 B = self.__background_magnetic_field
224 # this is for Cartesian (FIXME ?)
225 if self.__dim<3:
226 return np.array([B[0], B[2]])
227 else:
228 return np.array(B)
229
230 def getSetDensityMask(self):
231 """
232 Returns the density mask data object which is non-zero for cells
233 whose density value is fixed, zero otherwise.
234 """
235 z=self.getDomain().getX()[self.__dim-1]
236 m = whereNonNegative(z)
237 if self.__fix_density_below:
238 m += whereNonPositive(z+self.__fix_density_below)
239 return m
240
241 def getSetSusceptibilityMask(self):
242 """
243 Returns the susceptibility mask data object which is non-zero for
244 cells whose susceptibility value is fixed, zero otherwise.
245 """
246 z=self.getDomain().getX()[self.__dim-1]
247 m = whereNonNegative(z)
248 if self.__fix_susceptibility_below:
249 m += whereNonPositive(z+self.__fix_susceptibility_below)
250 return m
251
252 def getDomain(self):
253 """
254 Returns a domain that spans the data area plus padding.
255
256 The domain is created the first time this method is called,
257 subsequent calls return the same domain so anything that affects
258 the domain (such as padding) needs to be set beforehand.
259
260 :return: The escript domain for this data source
261 :rtype: `esys.escript.Domain`
262 """
263 if self.__domain is None:
264 self.__domain=self.__createDomain()
265 return self.__domain
266
267 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
268 """
269 This method sets the target domain parameters for the vertical
270 dimension.
271
272 :param depth: Depth of the domain (in meters)
273 :type depth: ``float``
274 :param air_layer: Depth of the layer above sea level (in meters)
275 :type air_layer: ``float``
276 :param num_cells: Number of domain elements for the entire vertical
277 dimension
278 :type num_cells: ``int``
279 """
280 if self.__domain is not None:
281 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
282 self._v_depth=depth
283 self._v_air_layer=air_layer
284 self._v_num_cells=num_cells
285
286 def __getTotalExtentsWithPadding(self):
287 """
288 Helper method that computes origin and number of elements
289 after adding padding to the bounding box of all available survey data.
290 """
291 X0, NX, DX = self.__getTotalExtents()
292 DIM=len(X0)
293 frac=[]
294 # padding is applied to each side so multiply by 2 to get the total
295 # amount of padding per dimension
296 pad, pt = self._padding
297 for i in range(DIM):
298 if pad[i] is None:
299 frac.append(0.)
300 continue
301 if pt == 'f': # fraction of side length
302 frac.append(2.*pad[i])
303 elif pt == 'e': # number of elements
304 frac.append(2.*pad[i]/float(NX[i]))
305 else: # absolute length
306 f=pad[i]/DX[i]
307 frac.append(2.*f/float(NX[i]))
308
309 # calculate new number of elements
310 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
311 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
312 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
313 return X0_padded, NX_padded, DX
314
315 def __getTotalExtents(self):
316 """
317 Helper method that computes the origin, number of elements and
318 minimal element spacing taking into account all available survey data.
319 """
320 if len(self.__sources)==0:
321 raise ValueError("No data")
322 X0, NE, DX = self.__sources[0].getDataExtents()
323 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
324
325 for src in self.__sources[1:]:
326 d_x0, d_ne, d_dx = src.getDataExtents()
327 for i in range(len(d_x0)):
328 X0[i]=min(X0[i], d_x0[i])
329 for i in range(len(d_dx)):
330 DX[i]=min(DX[i], d_dx[i])
331 for i in range(len(d_ne)):
332 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
333 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
334 return X0, NE, DX
335
336 def __createDomain(self):
337 """
338 Creates and returns an escript domain that spans the entire area of
339 available data plus a padding zone. This method is called only once
340 the first time `getDomain()` is invoked.
341
342 :return: The escript domain
343 :rtype: `esys.escript.Domain`
344 """
345 X0, NX, DX = self.__getTotalExtentsWithPadding()
346
347 # number of domain elements
348 NE = NX + [self._v_num_cells]
349
350 # origin of domain
351 origin = X0 + [-self._v_depth]
352 self._dom_origin = [np.floor(oi) for oi in origin]
353
354 # cell size / point spacing
355 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
356 self._spacing = [float(np.floor(si)) for si in spacing]
357
358 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
359 if self.__dim==3:
360 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
361 else:
362 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
363
364 # ripley may internally adjust NE and length, so recompute
365 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
366 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
367
368 self.logger.debug("Domain size: "+str(self._dom_NE))
369 self.logger.debug(" length: "+str(self._dom_len))
370 self.logger.debug(" origin: "+str(self._dom_origin))
371 return dom
372

  ViewVC Help
Powered by ViewVC 1.1.26