/[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 4213 - (hide 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 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 caltinay 4201 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 caltinay 4060
69 caltinay 4145 :param source: The data source to be added
70 caltinay 4060 :type source: `DataSource`
71     """
72 caltinay 4150 if self.__domain is not None:
73     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
74 caltinay 4060 if not isinstance(source, DataSource):
75     raise TypeError("source is not a DataSource")
76 caltinay 4201
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 caltinay 4145 self.__sources.append(source)
82 caltinay 4060
83 caltinay 4108 def setFractionalPadding(self, pad_x=None, pad_y=None):
84 caltinay 4060 """
85 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
86     dataset side lengths.
87    
88 caltinay 4150 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 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
91 caltinay 4060
92 caltinay 4108 :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 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
97 caltinay 4060 """
98 caltinay 4150 if self.__domain is not None:
99     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
100 caltinay 4108 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 caltinay 4213
112 caltinay 4108 def setPadding(self, pad_x=None, pad_y=None):
113     """
114     Sets the amount of padding around the dataset in absolute length units.
115 caltinay 4060
116 caltinay 4108 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 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
125 caltinay 4108 """
126 caltinay 4150 if self.__domain is not None:
127     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
128 caltinay 4108 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 caltinay 4145 Sets the amount of padding around the dataset in number of elements
139     (cells).
140 caltinay 4108
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 caltinay 4150 if self.__domain is not None:
152     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
153 caltinay 4108 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 caltinay 4213
165 caltinay 4060 def getGravitySurveys(self):
166 caltinay 4145 """
167     Returns a list of gravity surveys, see `getSurveys` for details.
168     """
169     return self.getSurveys(DataSource.GRAVITY)
170 caltinay 4213
171 gross 4115 def getMagneticSurveys(self):
172 gross 4120 """
173 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
174 gross 4120 """
175 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
176 caltinay 4213
177 caltinay 4145 def fixDensityBelow(self, depth=None):
178     """
179     Defines the depth below which the density anomaly is set to zero.
180     """
181 gross 4120 self.__fix_density_below=depth
182 caltinay 4213
183 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
184 gross 4120 """
185 caltinay 4145 Defines the depth below which the susceptibility anomaly is set to
186     zero.
187 gross 4120 """
188     self.__fix_susceptibility_below=depth
189    
190 gross 4115 def getSurveys(self, datatype):
191 caltinay 4060 """
192 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
193     available to this domain builder.
194 caltinay 4060
195 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
196 caltinay 4060 :rtype: ``list``
197     """
198 gross 4124 surveys=[]
199 caltinay 4145 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 gross 4124 return surveys
203 caltinay 4060
204 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
205 caltinay 4060 """
206 caltinay 4145 Sets the background magnetic flux density B=(B_r, B_theta, B_phi)
207 caltinay 4060 """
208 gross 4106 self.__background_magnetic_field=B
209 caltinay 4108
210 gross 4115 def getBackgroundMagneticFluxDensity(self):
211 gross 4106 """
212 caltinay 4145 Returns the background magnetic flux density.
213 caltinay 4108 """
214 gross 4106 B = self.__background_magnetic_field
215     # this is for Cartesian (FIXME ?)
216 caltinay 4145 if self.__dim<3:
217     return np.array([-B[2], -B[0]])
218 caltinay 4060 else:
219 caltinay 4145 return np.array([-B[1], -B[2], -B[0]])
220 caltinay 4060
221     def getSetDensityMask(self):
222 caltinay 4150 """
223     Returns the density mask data object which is non-zero for cells
224     whose density value is fixed, zero otherwise.
225     """
226 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
227 gross 4120 m = whereNonNegative(z)
228     if self.__fix_density_below:
229 caltinay 4145 m += whereNonPositive(z+self.__fix_density_below)
230 gross 4120 return m
231 caltinay 4145
232 caltinay 4060 def getSetSusceptibilityMask(self):
233 caltinay 4150 """
234     Returns the susceptibility mask data object which is non-zero for
235     cells whose susceptibility value is fixed, zero otherwise.
236     """
237 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
238 gross 4120 m = whereNonNegative(z)
239     if self.__fix_susceptibility_below:
240 caltinay 4145 m += whereNonPositive(z+self.__fix_susceptibility_below)
241 gross 4120 return m
242 caltinay 4060
243     def getDomain(self):
244     """
245     Returns a domain that spans the data area plus padding.
246    
247 caltinay 4145 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 caltinay 4060 :rtype: `esys.escript.Domain`
253     """
254 caltinay 4145 if self.__domain is None:
255     self.__domain=self.__createDomain()
256     return self.__domain
257 caltinay 4060
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 caltinay 4145 :param depth: Depth of the domain (in meters)
264 caltinay 4060 :type depth: ``float``
265 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
266 caltinay 4060 :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 caltinay 4150 if self.__domain is not None:
272     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
273 caltinay 4060 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 caltinay 4108 pad, pt = self._padding
288 caltinay 4060 for i in range(DIM):
289 azadeh 4127 if pad[i] is None:
290 caltinay 4145 frac.append(0.)
291     continue
292 caltinay 4108 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 caltinay 4060
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 caltinay 4145 if len(self.__sources)==0:
312 caltinay 4060 raise ValueError("No data")
313 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
314 caltinay 4060 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
315    
316 caltinay 4145 for src in self.__sources[1:]:
317 caltinay 4060 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 caltinay 4145 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 caltinay 4060 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 caltinay 4108
355 caltinay 4060 # ripley may internally adjust NE and length, so recompute
356 caltinay 4145 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 caltinay 4060
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