/[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 4343 - (hide 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 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 gross 4343 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 caltinay 4145 """
186 gross 4120 self.__fix_density_below=depth
187 caltinay 4213
188 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
189 gross 4120 """
190 gross 4343 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 gross 4120 """
197     self.__fix_susceptibility_below=depth
198    
199 gross 4115 def getSurveys(self, datatype):
200 caltinay 4060 """
201 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
202     available to this domain builder.
203 caltinay 4060
204 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
205 caltinay 4060 :rtype: ``list``
206     """
207 gross 4124 surveys=[]
208 caltinay 4145 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 gross 4124 return surveys
212 caltinay 4060
213 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
214 caltinay 4060 """
215 caltinay 4287 Sets the background magnetic flux density B=(B_North, B_East, B_Vertical)
216 caltinay 4060 """
217 gross 4106 self.__background_magnetic_field=B
218 caltinay 4108
219 gross 4115 def getBackgroundMagneticFluxDensity(self):
220 gross 4106 """
221 caltinay 4145 Returns the background magnetic flux density.
222 caltinay 4108 """
223 gross 4106 B = self.__background_magnetic_field
224     # this is for Cartesian (FIXME ?)
225 caltinay 4145 if self.__dim<3:
226 caltinay 4287 return np.array([B[0], B[2]])
227 caltinay 4060 else:
228 caltinay 4287 return np.array(B)
229 caltinay 4060
230     def getSetDensityMask(self):
231 caltinay 4150 """
232     Returns the density mask data object which is non-zero for cells
233     whose density value is fixed, zero otherwise.
234     """
235 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
236 gross 4120 m = whereNonNegative(z)
237     if self.__fix_density_below:
238 caltinay 4145 m += whereNonPositive(z+self.__fix_density_below)
239 gross 4120 return m
240 caltinay 4145
241 caltinay 4060 def getSetSusceptibilityMask(self):
242 caltinay 4150 """
243     Returns the susceptibility mask data object which is non-zero for
244     cells whose susceptibility value is fixed, zero otherwise.
245     """
246 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
247 gross 4120 m = whereNonNegative(z)
248     if self.__fix_susceptibility_below:
249 caltinay 4145 m += whereNonPositive(z+self.__fix_susceptibility_below)
250 gross 4120 return m
251 caltinay 4060
252     def getDomain(self):
253     """
254     Returns a domain that spans the data area plus padding.
255    
256 caltinay 4145 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 caltinay 4060 :rtype: `esys.escript.Domain`
262     """
263 caltinay 4145 if self.__domain is None:
264     self.__domain=self.__createDomain()
265     return self.__domain
266 caltinay 4060
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 caltinay 4145 :param depth: Depth of the domain (in meters)
273 caltinay 4060 :type depth: ``float``
274 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
275 caltinay 4060 :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 caltinay 4150 if self.__domain is not None:
281     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
282 caltinay 4060 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 caltinay 4108 pad, pt = self._padding
297 caltinay 4060 for i in range(DIM):
298 azadeh 4127 if pad[i] is None:
299 caltinay 4145 frac.append(0.)
300     continue
301 caltinay 4108 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 caltinay 4060
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 caltinay 4145 if len(self.__sources)==0:
321 caltinay 4060 raise ValueError("No data")
322 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
323 caltinay 4060 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
324    
325 caltinay 4145 for src in self.__sources[1:]:
326 caltinay 4060 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 caltinay 4145 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 caltinay 4060 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 caltinay 4108
364 caltinay 4060 # ripley may internally adjust NE and length, so recompute
365 caltinay 4145 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 caltinay 4060
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