/[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 4213 - (show annotations)
Tue Feb 19 01:16:29 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 13982 byte(s)
Some cleanup and more consistent logging.

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 zero.
180 """
181 self.__fix_density_below=depth
182
183 def fixSusceptibilityBelow(self, depth=None):
184 """
185 Defines the depth below which the susceptibility anomaly is set to
186 zero.
187 """
188 self.__fix_susceptibility_below=depth
189
190 def getSurveys(self, datatype):
191 """
192 Returns a list of `Data` objects for all surveys of type `datatype`
193 available to this domain builder.
194
195 :return: List of surveys which are tuples (anomaly,error).
196 :rtype: ``list``
197 """
198 surveys=[]
199 for src in self.__sources:
200 if src.getDataType()==datatype:
201 surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
202 return surveys
203
204 def setBackgroundMagneticFluxDensity(self, B):
205 """
206 Sets the background magnetic flux density B=(B_r, B_theta, B_phi)
207 """
208 self.__background_magnetic_field=B
209
210 def getBackgroundMagneticFluxDensity(self):
211 """
212 Returns the background magnetic flux density.
213 """
214 B = self.__background_magnetic_field
215 # this is for Cartesian (FIXME ?)
216 if self.__dim<3:
217 return np.array([-B[2], -B[0]])
218 else:
219 return np.array([-B[1], -B[2], -B[0]])
220
221 def getSetDensityMask(self):
222 """
223 Returns the density mask data object which is non-zero for cells
224 whose density value is fixed, zero otherwise.
225 """
226 z=self.getDomain().getX()[self.__dim-1]
227 m = whereNonNegative(z)
228 if self.__fix_density_below:
229 m += whereNonPositive(z+self.__fix_density_below)
230 return m
231
232 def getSetSusceptibilityMask(self):
233 """
234 Returns the susceptibility mask data object which is non-zero for
235 cells whose susceptibility value is fixed, zero otherwise.
236 """
237 z=self.getDomain().getX()[self.__dim-1]
238 m = whereNonNegative(z)
239 if self.__fix_susceptibility_below:
240 m += whereNonPositive(z+self.__fix_susceptibility_below)
241 return m
242
243 def getDomain(self):
244 """
245 Returns a domain that spans the data area plus padding.
246
247 The domain is created the first time this method is called,
248 subsequent calls return the same domain so anything that affects
249 the domain (such as padding) needs to be set beforehand.
250
251 :return: The escript domain for this data source
252 :rtype: `esys.escript.Domain`
253 """
254 if self.__domain is None:
255 self.__domain=self.__createDomain()
256 return self.__domain
257
258 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
259 """
260 This method sets the target domain parameters for the vertical
261 dimension.
262
263 :param depth: Depth of the domain (in meters)
264 :type depth: ``float``
265 :param air_layer: Depth of the layer above sea level (in meters)
266 :type air_layer: ``float``
267 :param num_cells: Number of domain elements for the entire vertical
268 dimension
269 :type num_cells: ``int``
270 """
271 if self.__domain is not None:
272 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
273 self._v_depth=depth
274 self._v_air_layer=air_layer
275 self._v_num_cells=num_cells
276
277 def __getTotalExtentsWithPadding(self):
278 """
279 Helper method that computes origin and number of elements
280 after adding padding to the bounding box of all available survey data.
281 """
282 X0, NX, DX = self.__getTotalExtents()
283 DIM=len(X0)
284 frac=[]
285 # padding is applied to each side so multiply by 2 to get the total
286 # amount of padding per dimension
287 pad, pt = self._padding
288 for i in range(DIM):
289 if pad[i] is None:
290 frac.append(0.)
291 continue
292 if pt == 'f': # fraction of side length
293 frac.append(2.*pad[i])
294 elif pt == 'e': # number of elements
295 frac.append(2.*pad[i]/float(NX[i]))
296 else: # absolute length
297 f=pad[i]/DX[i]
298 frac.append(2.*f/float(NX[i]))
299
300 # calculate new number of elements
301 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
302 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
303 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
304 return X0_padded, NX_padded, DX
305
306 def __getTotalExtents(self):
307 """
308 Helper method that computes the origin, number of elements and
309 minimal element spacing taking into account all available survey data.
310 """
311 if len(self.__sources)==0:
312 raise ValueError("No data")
313 X0, NE, DX = self.__sources[0].getDataExtents()
314 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
315
316 for src in self.__sources[1:]:
317 d_x0, d_ne, d_dx = src.getDataExtents()
318 for i in range(len(d_x0)):
319 X0[i]=min(X0[i], d_x0[i])
320 for i in range(len(d_dx)):
321 DX[i]=min(DX[i], d_dx[i])
322 for i in range(len(d_ne)):
323 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
324 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
325 return X0, NE, DX
326
327 def __createDomain(self):
328 """
329 Creates and returns an escript domain that spans the entire area of
330 available data plus a padding zone. This method is called only once
331 the first time `getDomain()` is invoked.
332
333 :return: The escript domain
334 :rtype: `esys.escript.Domain`
335 """
336 X0, NX, DX = self.__getTotalExtentsWithPadding()
337
338 # number of domain elements
339 NE = NX + [self._v_num_cells]
340
341 # origin of domain
342 origin = X0 + [-self._v_depth]
343 self._dom_origin = [np.floor(oi) for oi in origin]
344
345 # cell size / point spacing
346 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
347 self._spacing = [float(np.floor(si)) for si in spacing]
348
349 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
350 if self.__dim==3:
351 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
352 else:
353 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
354
355 # ripley may internally adjust NE and length, so recompute
356 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
357 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
358
359 self.logger.debug("Domain size: "+str(self._dom_NE))
360 self.logger.debug(" length: "+str(self._dom_len))
361 self.logger.debug(" origin: "+str(self._dom_origin))
362 return dom
363

  ViewVC Help
Powered by ViewVC 1.1.26