/[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 4106 - (show annotations)
Thu Dec 13 05:08:44 2012 UTC (6 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 9718 byte(s)
some clarification on how to handle back ground magnetic field plus renaming of SyntheticData
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 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-2012 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._gravity_surveys=[]
56 self._magnetic_surveys=[]
57 self._sources=[]
58 self.setPadding()
59 self.setVerticalExtents()
60 self.__background_magnetic_field = None
61
62 def addSource(self, source):
63 """
64 Adds a survey data provider to the domain builder.
65
66 :param source: The data source to be added.
67 :type source: `DataSource`
68 """
69 if not isinstance(source, DataSource):
70 raise TypeError("source is not a DataSource")
71 self._sources.append(source)
72
73 def setPadding(self, pad_x=10, pad_y=10):
74 """
75 Sets the amount of padding around the dataset. If `pad_x`/`pad_y`
76 is >=1 the value is treated as number of elements to be added to the
77 domain (per side).
78 If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the
79 dataset size. For example, calling ``setPadding(3, 0.1)`` to a data
80 source with size 10x20 will result in the padded data set size
81 16x24 (10+2*3, 20*(1+2*0.1))
82
83 :param pad_x: Padding per side in x direction (default: 10 elements)
84 :type pad_x: ``int`` or ``float``
85 :param pad_y: Padding per side in y direction (default: 10 elements).
86 This value is only used for 3-dimensional datasets
87 :type pad_y: ``int`` or ``float``
88 """
89 self._padding=[pad_x,pad_y]
90
91 def getGravitySurveys(self):
92 """
93 Returns a list of `Data` objects for all gravity surveys available to
94 this domain builder.
95
96 :return: List of gravity surveys which are tuples (anomaly,error).
97 :rtype: ``list``
98 """
99 if len(self._gravity_surveys)==0:
100 for src in self._sources:
101 if src.getDataType()==DataSource.GRAVITY:
102 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
103 self._gravity_surveys.append(survey)
104 return self._gravity_surveys
105
106 def getMagneticSurveys(self):
107 """
108 Returns a list of `Data` objects for all magnetic surveys available to
109 this domain builder.
110
111 :return: List of magnetic surveys which are tuples (anomaly,error).
112 :rtype: ``list``
113 """
114 if len(self._magnetic_surveys)==0:
115 for src in self._sources:
116 if src.getDataType()==DataSource.MAGNETIC:
117 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
118 self._magnetic_surveys.append(survey)
119 return self._magnetic_surveys
120
121 def setBackgroundMagneticField(self, B):
122 """
123 sets the back ground magnetic field B=(B_r,B_theta, B_phi)
124 """
125 self.__background_magnetic_field=B
126
127 def getBackgroundMagneticField(self):
128 """
129 returns the back ground magnetic field.
130 """
131 B = self.__background_magnetic_field
132 # this is for Cartesian (FIXME ?)
133 if self._dim<3:
134 return np.array([-B[2], -B[0]])
135 else:
136 return np.array([-B[1], -B[2], -B[0]])
137
138 def getSetDensityMask(self):
139 x=self.getDomain().getX()
140 return wherePositive(x[self._dim-1])+whereZero(x[self._dim-1]-inf(x[self._dim-1]))
141 # \ + whereZero(x[0]-inf(x[0]))+ whereZero(x[0]-sup(x[0]))
142
143 def getSetSusceptibilityMask(self):
144 return wherePositive(self.getDomain().getX()[self._dim-1])
145
146 def getDomain(self):
147 """
148 Returns a domain that spans the data area plus padding.
149 The domain is created the first time this method is called, subsequent
150 calls return the same domain so anything that affects the domain
151 (such as padding) needs to be set beforehand.
152
153 :return: The escript domain for this data source.
154 :rtype: `esys.escript.Domain`
155 """
156 if self._domain is None:
157 self._domain=self.__createDomain()
158 return self._domain
159
160 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
161 """
162 This method sets the target domain parameters for the vertical
163 dimension.
164
165 :param depth: Depth in meters of the domain.
166 :type depth: ``float``
167 :param air_layer: Depth of the layer above sea level in meters
168 :type air_layer: ``float``
169 :param num_cells: Number of domain elements for the entire vertical
170 dimension
171 :type num_cells: ``int``
172 """
173 self._v_depth=depth
174 self._v_air_layer=air_layer
175 self._v_num_cells=num_cells
176
177 def __getTotalExtentsWithPadding(self):
178 """
179 Helper method that computes origin and number of elements
180 after adding padding to the bounding box of all available survey data.
181 """
182 X0, NX, DX = self.__getTotalExtents()
183 DIM=len(X0)
184 frac=[]
185 # padding is applied to each side so multiply by 2 to get the total
186 # amount of padding per dimension
187 for i in range(DIM):
188 if self._padding[i]>=0 and self._padding[i]<1:
189 frac.append(2.*self._padding[i])
190 else:
191 frac.append(2.*self._padding[i]/float(NX[i]))
192
193 # calculate new number of elements
194 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
195 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
196 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
197 return X0_padded, NX_padded, DX
198
199 def __getTotalExtents(self):
200 """
201 Helper method that computes the origin, number of elements and
202 minimal element spacing taking into account all available survey data.
203 """
204 if len(self._sources)==0:
205 raise ValueError("No data")
206 X0, NE, DX = self._sources[0].getDataExtents()
207 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
208
209 for src in self._sources[1:]:
210 d_x0, d_ne, d_dx = src.getDataExtents()
211 for i in range(len(d_x0)):
212 X0[i]=min(X0[i], d_x0[i])
213 for i in range(len(d_dx)):
214 DX[i]=min(DX[i], d_dx[i])
215 for i in range(len(d_ne)):
216 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
217 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
218 return X0, NE, DX
219
220 def __createDomain(self):
221 """
222 Creates and returns an escript domain that spans the entire area of
223 available data plus a padding zone. This method is called only once
224 the first time `getDomain()` is invoked.
225
226 :return: The escript domain
227 :rtype: `esys.escript.Domain`
228 """
229 X0, NX, DX = self.__getTotalExtentsWithPadding()
230
231 # number of domain elements
232 NE = NX + [self._v_num_cells]
233
234 # origin of domain
235 origin = X0 + [-self._v_depth]
236 self._dom_origin = [np.floor(oi) for oi in origin]
237
238 # cell size / point spacing
239 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
240 self._spacing = [float(np.floor(si)) for si in spacing]
241
242 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
243 if self._dim==3:
244 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
245 else:
246 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
247
248 # ripley may internally adjust NE and length, so recompute
249 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
250 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
251
252 self.logger.debug("Domain size: "+str(self._dom_NE))
253 self.logger.debug(" length: "+str(self._dom_len))
254 self.logger.debug(" origin: "+str(self._dom_origin))
255 return dom
256

  ViewVC Help
Powered by ViewVC 1.1.26