/[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 4362 - (hide 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 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 caltinay 4362 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 caltinay 4060
73 caltinay 4145 :param source: The data source to be added
74 caltinay 4060 :type source: `DataSource`
75     """
76 caltinay 4150 if self.__domain is not None:
77     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
78 caltinay 4060 if not isinstance(source, DataSource):
79     raise TypeError("source is not a DataSource")
80 caltinay 4201
81 caltinay 4362 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 caltinay 4201 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 caltinay 4145 self.__sources.append(source)
89 caltinay 4060
90 caltinay 4108 def setFractionalPadding(self, pad_x=None, pad_y=None):
91 caltinay 4060 """
92 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
93     dataset side lengths.
94    
95 caltinay 4150 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 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
98 caltinay 4060
99 caltinay 4108 :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 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
104 caltinay 4060 """
105 caltinay 4150 if self.__domain is not None:
106     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
107 caltinay 4108 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 caltinay 4213
119 caltinay 4108 def setPadding(self, pad_x=None, pad_y=None):
120     """
121     Sets the amount of padding around the dataset in absolute length units.
122 caltinay 4060
123 caltinay 4108 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 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
132 caltinay 4108 """
133 caltinay 4150 if self.__domain is not None:
134     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
135 caltinay 4108 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 caltinay 4145 Sets the amount of padding around the dataset in number of elements
146     (cells).
147 caltinay 4108
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 caltinay 4150 if self.__domain is not None:
159     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
160 caltinay 4108 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 caltinay 4213
172 caltinay 4060 def getGravitySurveys(self):
173 caltinay 4145 """
174     Returns a list of gravity surveys, see `getSurveys` for details.
175     """
176     return self.getSurveys(DataSource.GRAVITY)
177 caltinay 4213
178 gross 4115 def getMagneticSurveys(self):
179 gross 4120 """
180 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
181 gross 4120 """
182 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
183 caltinay 4213
184 caltinay 4145 def fixDensityBelow(self, depth=None):
185     """
186 gross 4343 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 caltinay 4145 """
193 gross 4120 self.__fix_density_below=depth
194 caltinay 4213
195 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
196 gross 4120 """
197 gross 4343 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 gross 4120 """
204     self.__fix_susceptibility_below=depth
205    
206 gross 4115 def getSurveys(self, datatype):
207 caltinay 4060 """
208 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
209     available to this domain builder.
210 caltinay 4060
211 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
212 caltinay 4060 :rtype: ``list``
213     """
214 gross 4124 surveys=[]
215 caltinay 4145 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 gross 4124 return surveys
219 caltinay 4060
220 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
221 caltinay 4060 """
222 caltinay 4287 Sets the background magnetic flux density B=(B_North, B_East, B_Vertical)
223 caltinay 4060 """
224 gross 4106 self.__background_magnetic_field=B
225 caltinay 4108
226 gross 4115 def getBackgroundMagneticFluxDensity(self):
227 gross 4106 """
228 caltinay 4145 Returns the background magnetic flux density.
229 caltinay 4108 """
230 gross 4106 B = self.__background_magnetic_field
231     # this is for Cartesian (FIXME ?)
232 caltinay 4145 if self.__dim<3:
233 caltinay 4287 return np.array([B[0], B[2]])
234 caltinay 4060 else:
235 caltinay 4287 return np.array(B)
236 caltinay 4060
237     def getSetDensityMask(self):
238 caltinay 4150 """
239     Returns the density mask data object which is non-zero for cells
240     whose density value is fixed, zero otherwise.
241     """
242 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
243 gross 4120 m = whereNonNegative(z)
244     if self.__fix_density_below:
245 caltinay 4145 m += whereNonPositive(z+self.__fix_density_below)
246 gross 4120 return m
247 caltinay 4145
248 caltinay 4060 def getSetSusceptibilityMask(self):
249 caltinay 4150 """
250     Returns the susceptibility mask data object which is non-zero for
251     cells whose susceptibility value is fixed, zero otherwise.
252     """
253 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
254 gross 4120 m = whereNonNegative(z)
255     if self.__fix_susceptibility_below:
256 caltinay 4145 m += whereNonPositive(z+self.__fix_susceptibility_below)
257 gross 4120 return m
258 caltinay 4060
259     def getDomain(self):
260     """
261     Returns a domain that spans the data area plus padding.
262    
263 caltinay 4145 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 caltinay 4060 :rtype: `esys.escript.Domain`
269     """
270 caltinay 4145 if self.__domain is None:
271     self.__domain=self.__createDomain()
272     return self.__domain
273 caltinay 4060
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 caltinay 4145 :param depth: Depth of the domain (in meters)
280 caltinay 4060 :type depth: ``float``
281 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
282 caltinay 4060 :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 caltinay 4150 if self.__domain is not None:
288     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
289 caltinay 4060 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 caltinay 4362 Helper method that computes origin and number of data elements
296 caltinay 4060 after adding padding to the bounding box of all available survey data.
297     """
298     X0, NX, DX = self.__getTotalExtents()
299 caltinay 4362 DATA_DIM=len(X0)
300 caltinay 4060 frac=[]
301     # padding is applied to each side so multiply by 2 to get the total
302     # amount of padding per dimension
303 caltinay 4108 pad, pt = self._padding
304 caltinay 4362 for i in range(DATA_DIM):
305 azadeh 4127 if pad[i] is None:
306 caltinay 4145 frac.append(0.)
307     continue
308 caltinay 4108 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 caltinay 4060
316     # calculate new number of elements
317 caltinay 4362 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 caltinay 4060 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 caltinay 4145 if len(self.__sources)==0:
328 caltinay 4060 raise ValueError("No data")
329 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
330 caltinay 4060 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
331    
332 caltinay 4145 for src in self.__sources[1:]:
333 caltinay 4060 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 caltinay 4145 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 caltinay 4060 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 caltinay 4108
371 caltinay 4060 # ripley may internally adjust NE and length, so recompute
372 caltinay 4145 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 caltinay 4060
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