/[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 4106 - (hide 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 caltinay 4060
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 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     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 gross 4106 self.__background_magnetic_field = None
61 caltinay 4060
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 gross 4106 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 caltinay 4060 def getBackgroundMagneticField(self):
128 gross 4106 """
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 caltinay 4060 else:
136 gross 4106 return np.array([-B[1], -B[2], -B[0]])
137 caltinay 4060
138     def getSetDensityMask(self):
139 gross 4102 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 caltinay 4060
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