/[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 4362 - (show annotations)
Tue Apr 16 04:37:13 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 14933 byte(s)
Implemented NumpyData downunder source (w/ tests).
Added check of data dimensionality to domainbuilder since data_dim must equal
domain_dim-1 at the moment.

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

  ViewVC Help
Powered by ViewVC 1.1.26