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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4145 - (hide annotations)
Fri Jan 18 00:51:49 2013 UTC (6 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 12635 byte(s)
Checkpoint.

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

  ViewVC Help
Powered by ViewVC 1.1.26