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

  ViewVC Help
Powered by ViewVC 1.1.26