/[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 4729 - (hide annotations)
Sat Mar 8 00:56:37 2014 UTC (5 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 20943 byte(s)
a more general for of the inversioncostfunction has required for waveform inversin
1 caltinay 4060
2     ##############################################################################
3     #
4 jfenwick 4657 # Copyright (c) 2003-2014 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 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 caltinay 4060 #
15     ##############################################################################
16    
17     """Domain construction from survey data for inversions"""
18    
19 jfenwick 4657 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
20 caltinay 4060 http://www.uq.edu.au
21     Primary Business: Queensland, Australia"""
22     __license__="""Licensed under the Open Software License version 3.0
23     http://www.opensource.org/licenses/osl-3.0.php"""
24     __url__="https://launchpad.net/escript-finley"
25    
26     __all__ = ['DomainBuilder']
27    
28     import logging
29     import numpy as np
30     from esys.escript.util import *
31 caltinay 4097 from esys.escript import unitsSI as U
32 caltinay 4060 from esys.ripley import Brick, Rectangle
33     from .datasources import DataSource
34 gross 4373 from .coordinates import ReferenceSystem, CartesianReferenceSystem
35 caltinay 4060
36 gross 4729
37 caltinay 4060 class DomainBuilder(object):
38     """
39     This class is responsible for constructing an escript Domain object with
40     suitable extents and resolution for survey data (`DataSource` objects)
41 caltinay 4435 that are added to it.
42    
43     The domain covers a region above and below the Earth surface. The
44     East-West direction is used as the x- or longitudinal or x[0] direction,
45     the North-South direction is used as the y- or latitudinal or x[1]
46     direction, the vertical direction is denoted by z or radial or x[2]
47     direction. The corresponding terms are used synonymously.
48 caltinay 4060 """
49    
50 gross 4373 def __init__(self, dim=3, reference_system=None):
51 caltinay 4060 """
52     Constructor.
53    
54     :param dim: Dimensionality (2 or 3) of the target domain.
55     This has implications for the survey data than can be
56     added. By default a 3D domain is created.
57     :type dim: ``int``
58 gross 4373 :param reference_system: reference coordinate system. By default the
59     Cartesian coordinate system is used.
60     :type reference_system: `ReferenceSystem`
61 caltinay 4060 """
62     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
63     if dim not in (2,3):
64     raise ValueError("Number of dimensions must be 2 or 3")
65 caltinay 4394 if not reference_system:
66     self.__reference_system=CartesianReferenceSystem()
67     else:
68     self.__reference_system=reference_system
69 caltinay 4435
70 caltinay 4436 if self.__reference_system.isCartesian():
71     self.__v_scale=1.
72     else:
73     self.__v_scale=1./self.getReferenceSystem().getHeightUnit()
74    
75 caltinay 4145 self.__domain=None
76     self.__dim=dim
77     self.__sources=[]
78     self.__background_magnetic_field=None
79 gross 4729 self.__tags=[] # list of all tags being used by all data sources beeing attached:
80 gross 4373 self.setElementPadding()
81 caltinay 4060 self.setVerticalExtents()
82 gross 4120 self.fixDensityBelow()
83     self.fixSusceptibilityBelow()
84 gross 4729 self.fixVelocityBelow()
85 gross 4373
86     def getReferenceSystem(self):
87     """
88     returns the reference coordinate system
89    
90     :rtype: `ReferenceSystem`
91     """
92     return self.__reference_system
93 gross 4729 def getTags(self):
94     """
95     returns a list of all tags beeing used by the attached data sources.
96     the list may be empty.
97     """
98     return self.__tags
99    
100 caltinay 4060 def addSource(self, source):
101     """
102     Adds a survey data provider to the domain builder.
103 caltinay 4201 An exception is raised if the domain has already been built or if the
104     UTM zone of `source` does not match the UTM zone of sources already
105 caltinay 4362 added to the domain builder (see Inversion Cookbook for more
106     information). An exception is also raised if the dimensionality of the
107     data source is incompatible with this domain builder. That is, the
108     dimensionality of the data must be one less than the dimensionality
109     of the domain (specified in the constructor).
110 caltinay 4060
111 gross 4373 :param source: The data source to be added. Its reference system needs
112     to match the reference system of the DomainBuilder.
113 caltinay 4060 :type source: `DataSource`
114     """
115 caltinay 4150 if self.__domain is not None:
116     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
117 caltinay 4060 if not isinstance(source, DataSource):
118     raise TypeError("source is not a DataSource")
119 caltinay 4394 if not source.getReferenceSystem() == self.getReferenceSystem():
120     raise ValueError("source reference system does not match.")
121 caltinay 4201
122 caltinay 4362 DATA_DIM = len(source.getDataExtents()[0])
123     if DATA_DIM != self.__dim-1:
124     raise ValueError("Data must be %d-dimensional."%(self.__dim-1))
125 caltinay 4201 if len(self.__sources)>0:
126     if self.__sources[0].getUtmZone() != source.getUtmZone():
127     raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")
128    
129 caltinay 4145 self.__sources.append(source)
130 gross 4729 if source.getTags(): self.__tags=list(set(self.__tags + source.getTags()))
131 caltinay 4060
132 gross 4373 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
133 caltinay 4060 """
134 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
135     dataset side lengths.
136    
137 caltinay 4150 For example, calling ``setFractionalPadding(0.2, 0.1)`` with a data
138     source of size 10x20 will result in the padded data set size
139 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
140 caltinay 4060
141 caltinay 4108 :param pad_x: Padding per side in x direction (default: no padding)
142     :type pad_x: ``float``
143     :param pad_y: Padding per side in y direction (default: no padding)
144     :type pad_y: ``float``
145 gross 4373 :param pad_lat: Padding per side in latitudinal direction (default: no padding)
146     :type pad_lat: ``float``
147     :param pad_lon: Padding per side in longitudinal direction (default: no padding)
148     :type pad_lon: ``float``
149     :note: `pad_y` is ignored for 2-dimensional domains.
150 caltinay 4060 """
151 gross 4373 if not pad_lat == None:
152 caltinay 4394 if not pad_x == None:
153     raise ValueError("Either pad_lat or pad_x can be set.")
154     else:
155     pad_x = pad_lat
156 gross 4373 if not pad_lon == None:
157 caltinay 4394 if not pad_y == None:
158     raise ValueError("Either pad_lon or pad_y can be set.")
159     else:
160     pad_y = pad_lan
161 caltinay 4150 if self.__domain is not None:
162     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
163 caltinay 4108 if pad_x is not None:
164     if pad_x < 0:
165     raise ValueError("setFractionalPadding: Arguments must be non-negative")
166     if pad_x > 10:
167     raise ValueError("setFractionalPadding: Argument too large")
168     if pad_y is not None:
169     if pad_y < 0:
170     raise ValueError("setFractionalPadding: Arguments must be non-negative")
171     if pad_y > 10:
172     raise ValueError("setFractionalPadding: Argument too large")
173     self._padding = [pad_x,pad_y], 'f'
174 caltinay 4213
175 gross 4376 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
176 caltinay 4108 """
177     Sets the amount of padding around the dataset in absolute length units.
178 caltinay 4060
179 caltinay 4108 The final domain size will be the length in x (in y) of the dataset
180     plus twice the value of `pad_x` (`pad_y`). The arguments must be
181     non-negative.
182    
183     :param pad_x: Padding per side in x direction (default: no padding)
184 gross 4373 :type pad_x: ``float`` in units of length (meter)
185 caltinay 4108 :param pad_y: Padding per side in y direction (default: no padding)
186 gross 4373 :type pad_y: ``float`` in units of length (meter)
187 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
188 gross 4373 :note: this function can only be used if the reference system is Cartesian
189 caltinay 4108 """
190 gross 4373 if not self.getReferenceSystem().isCartesian():
191 caltinay 4394 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
192 caltinay 4150 if self.__domain is not None:
193     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
194 caltinay 4108 if pad_x is not None:
195     if pad_x < 0:
196     raise ValueError("setPadding: Arguments must be non-negative")
197     if pad_y is not None:
198     if pad_y < 0:
199     raise ValueError("setPadding: Arguments must be non-negative")
200     self._padding = [pad_x,pad_y], 'l'
201 gross 4373
202     def setGeoPadding(self, pad_lat=None, pad_lon=None):
203     """
204     Sets the amount of padding around the dataset in longitude and latitude.
205 caltinay 4108
206 caltinay 4435 The final domain size will be the extent in the latitudinal (in
207     longitudinal) direction of the dataset plus twice the value of
208     `pad_lat` (`pad_lon`). The arguments must be non-negative.
209 gross 4373
210 caltinay 4435 :param pad_lat: Padding per side in latitudinal direction (default: 0)
211 gross 4373 :type pad_lat: ``float`` in units of degree
212 caltinay 4435 :param pad_lon: Padding per side in longitudinal direction (default: 0)
213 gross 4373 :type pad_lon: ``float`` in units of degree
214     :note: `pad_lon` is ignored for 2-dimensional domains.
215     :note: this function can only be used if the reference system is not Cartesian
216 caltinay 4108 """
217 gross 4373 if self.getReferenceSystem().isCartesian():
218 caltinay 4394 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
219 gross 4373 if self.__domain is not None:
220     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
221     if pad_lat is not None:
222     if pad_lat < 0:
223     raise ValueError("setPadding: Arguments must be non-negative")
224     if pad_lon is not None:
225     if pad_lon < 0:
226     raise ValueError("setPadding: Arguments must be non-negative")
227     self._padding = [pad_lat,pad_lon], 'd'
228    
229     def setElementPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
230     """
231 caltinay 4145 Sets the amount of padding around the dataset in number of elements
232     (cells).
233 caltinay 4108
234     When the domain is constructed `pad_x` (`pad_y`) elements are added
235     on each side of the x- (y-) dimension. The arguments must be
236     non-negative.
237    
238     :param pad_x: Padding per side in x direction (default: no padding)
239     :type pad_x: ``int``
240     :param pad_y: Padding per side in y direction (default: no padding)
241     :type pad_y: ``int``
242     :note: `pad_y` is ignored for 2-dimensional datasets.
243     """
244 gross 4373 if not pad_lat == None:
245 caltinay 4394 if not pad_x == None:
246     raise ValueError("Either pad_lat or pad_x can be set.")
247     else:
248     pad_x = pad_lat
249 gross 4373 if not pad_lon == None:
250 caltinay 4394 if not pad_y == None:
251     raise ValueError("Either pad_lon or pad_y can be set.")
252     else:
253     pad_y = pad_lan
254    
255 caltinay 4150 if self.__domain is not None:
256     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
257 caltinay 4108 if pad_x is not None:
258     if type(pad_x) is not int:
259     raise TypeError("setElementPadding expects integer arguments")
260     if pad_x < 0:
261     raise ValueError("setElementPadding: Arguments must be non-negative")
262     if pad_y is not None:
263     if type(pad_y) is not int:
264     raise TypeError("setElementPadding expects integer arguments")
265     if pad_y < 0:
266     raise ValueError("setElementPadding: Arguments must be non-negative")
267     self._padding = [pad_x,pad_y], 'e'
268 caltinay 4213
269 caltinay 4060 def getGravitySurveys(self):
270 caltinay 4145 """
271     Returns a list of gravity surveys, see `getSurveys` for details.
272     """
273     return self.getSurveys(DataSource.GRAVITY)
274 caltinay 4213
275 gross 4115 def getMagneticSurveys(self):
276 gross 4120 """
277 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
278 gross 4120 """
279 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
280 caltinay 4213
281 gross 4729
282 caltinay 4145 def fixDensityBelow(self, depth=None):
283     """
284 caltinay 4435 Defines the depth below which the density anomaly is set to a given
285     value. If no value is given zero is assumed.
286 gross 4343
287 caltinay 4435 :param depth: depth below which the density is fixed. If not set, no
288     constraint at depth is applied.
289 gross 4343 :type depth: ``float``
290 caltinay 4145 """
291 gross 4120 self.__fix_density_below=depth
292 caltinay 4213
293 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
294 gross 4120 """
295 caltinay 4435 Defines the depth below which the susceptibility anomaly is set to a
296     given value. If no value is given zero is assumed.
297 gross 4343
298 caltinay 4435 :param depth: depth below which the susceptibility is fixed. If not
299     set, no constraint at depth is applied.
300 gross 4343 :type depth: ``float``
301 gross 4120 """
302     self.__fix_susceptibility_below=depth
303    
304 gross 4729 def fixVelocityBelow(self, depth=None):
305 caltinay 4060 """
306 gross 4729 Defines the depth below which the velocity and Q index is set to a
307     given value. If no value is given zero is assumed.
308    
309     :param depth: depth below which the velocity is fixed. If not
310     set, no constraint at depth is applied.
311     :type depth: ``float``
312     """
313     self.__fix_velocity_below=depth
314    
315    
316     def getSurveys(self, datatype, tags=None):
317     """
318 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
319 gross 4729 available to this domain builder. If a list of `tags` is given
320     only data sources whose tag matching the tag list are returned
321 caltinay 4060
322 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
323 caltinay 4060 :rtype: ``list``
324     """
325 gross 4124 surveys=[]
326 caltinay 4145 for src in self.__sources:
327     if src.getDataType()==datatype:
328 gross 4729 if tags is None or ( src.getTags() is not None and all( [ t in tags for t in src.getTags() ] ) ) :
329     surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
330 gross 4124 return surveys
331 caltinay 4060
332 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
333 caltinay 4060 """
334 gross 4433 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
335 caltinay 4060 """
336 gross 4106 self.__background_magnetic_field=B
337 caltinay 4108
338 gross 4115 def getBackgroundMagneticFluxDensity(self):
339 gross 4106 """
340 caltinay 4145 Returns the background magnetic flux density.
341 caltinay 4108 """
342 gross 4106 B = self.__background_magnetic_field
343 caltinay 4462 if B is None:
344     raise ValueError("No background magnetic flux density set!")
345    
346 gross 4433 if self.__dim < 3 :
347 caltinay 4287 return np.array([B[0], B[2]])
348 caltinay 4060 else:
349 caltinay 4287 return np.array(B)
350 caltinay 4060
351     def getSetDensityMask(self):
352 caltinay 4150 """
353     Returns the density mask data object which is non-zero for cells
354     whose density value is fixed, zero otherwise.
355     """
356 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
357 gross 4120 m = whereNonNegative(z)
358     if self.__fix_density_below:
359 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
360 gross 4120 return m
361 caltinay 4145
362 caltinay 4060 def getSetSusceptibilityMask(self):
363 caltinay 4150 """
364     Returns the susceptibility mask data object which is non-zero for
365     cells whose susceptibility value is fixed, zero otherwise.
366     """
367 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
368 gross 4120 m = whereNonNegative(z)
369     if self.__fix_susceptibility_below:
370 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
371 gross 4120 return m
372 caltinay 4060
373     def getDomain(self):
374     """
375     Returns a domain that spans the data area plus padding.
376    
377 caltinay 4145 The domain is created the first time this method is called,
378     subsequent calls return the same domain so anything that affects
379     the domain (such as padding) needs to be set beforehand.
380    
381     :return: The escript domain for this data source
382 caltinay 4060 :rtype: `esys.escript.Domain`
383     """
384 caltinay 4145 if self.__domain is None:
385     self.__domain=self.__createDomain()
386     return self.__domain
387 caltinay 4060
388     def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
389     """
390     This method sets the target domain parameters for the vertical
391     dimension.
392    
393 caltinay 4145 :param depth: Depth of the domain (in meters)
394 caltinay 4060 :type depth: ``float``
395 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
396 caltinay 4060 :type air_layer: ``float``
397     :param num_cells: Number of domain elements for the entire vertical
398     dimension
399     :type num_cells: ``int``
400     """
401 caltinay 4150 if self.__domain is not None:
402     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
403 caltinay 4060 self._v_depth=depth
404     self._v_air_layer=air_layer
405     self._v_num_cells=num_cells
406    
407     def __getTotalExtentsWithPadding(self):
408     """
409 caltinay 4362 Helper method that computes origin and number of data elements
410 caltinay 4060 after adding padding to the bounding box of all available survey data.
411     """
412     X0, NX, DX = self.__getTotalExtents()
413 caltinay 4362 DATA_DIM=len(X0)
414 caltinay 4060 frac=[]
415     # padding is applied to each side so multiply by 2 to get the total
416     # amount of padding per dimension
417 caltinay 4108 pad, pt = self._padding
418 caltinay 4362 for i in range(DATA_DIM):
419 azadeh 4127 if pad[i] is None:
420 caltinay 4145 frac.append(0.)
421     continue
422 gross 4373 if pt == 'f' : # fraction of side length
423 caltinay 4108 frac.append(2.*pad[i])
424     elif pt == 'e': # number of elements
425     frac.append(2.*pad[i]/float(NX[i]))
426     else: # absolute length
427     f=pad[i]/DX[i]
428     frac.append(2.*f/float(NX[i]))
429 caltinay 4060
430     # calculate new number of elements
431 caltinay 4362 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DATA_DIM)]
432     NXdiff=[NX_padded[i]-NX[i] for i in range(DATA_DIM)]
433     X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DATA_DIM)]
434 caltinay 4060 return X0_padded, NX_padded, DX
435    
436     def __getTotalExtents(self):
437     """
438     Helper method that computes the origin, number of elements and
439     minimal element spacing taking into account all available survey data.
440     """
441 caltinay 4145 if len(self.__sources)==0:
442 caltinay 4060 raise ValueError("No data")
443 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
444 caltinay 4613 # do not mess with the values if only one source used
445     if len(self.__sources)>1:
446     XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
447 caltinay 4060
448 caltinay 4613 for src in self.__sources[1:]:
449     d_x0, d_ne, d_dx = src.getDataExtents()
450     for i in range(len(d_x0)):
451     X0[i]=min(X0[i], d_x0[i])
452     for i in range(len(d_dx)):
453     DX[i]=min(DX[i], d_dx[i])
454     for i in range(len(d_ne)):
455     XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
456     # FIXME: should this be rounded up instead?
457     NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
458 caltinay 4060 return X0, NE, DX
459    
460     def __createDomain(self):
461     """
462     Creates and returns an escript domain that spans the entire area of
463     available data plus a padding zone. This method is called only once
464     the first time `getDomain()` is invoked.
465    
466     :return: The escript domain
467     :rtype: `esys.escript.Domain`
468     """
469     X0, NX, DX = self.__getTotalExtentsWithPadding()
470    
471     # number of domain elements
472     NE = NX + [self._v_num_cells]
473    
474     # origin of domain
475 caltinay 4435 origin = X0 + [-self._v_depth*self.__v_scale]
476    
477 caltinay 4394 if self.getReferenceSystem().isCartesian():
478     # rounding will give us about meter-accuracy with UTM coordinates
479     self._dom_origin = [np.floor(oi) for oi in origin]
480     else:
481     # this should give us about meter-accuracy with lat/lon coords
482     self._dom_origin = [1e-5*np.floor(oi*1e5) for oi in origin]
483 caltinay 4435
484 caltinay 4060 # cell size / point spacing
485 caltinay 4435 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
486 gross 4373 #self._spacing = [float(np.floor(si)) for si in spacing]
487     self._spacing = spacing
488 gross 4433
489 caltinay 4435 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
490 caltinay 4613
491 caltinay 4145 if self.__dim==3:
492 caltinay 4060 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
493     else:
494 caltinay 4434 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
495 caltinay 4108
496 caltinay 4060 # ripley may internally adjust NE and length, so recompute
497 caltinay 4145 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
498 caltinay 4435 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
499 caltinay 4060
500     self.logger.debug("Domain size: "+str(self._dom_NE))
501     self.logger.debug(" length: "+str(self._dom_len))
502     self.logger.debug(" origin: "+str(self._dom_origin))
503     return dom
504    

  ViewVC Help
Powered by ViewVC 1.1.26