/[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 4729 - (show 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
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 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 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 """Domain construction from survey data for inversions"""
18
19 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
20 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 from esys.escript import unitsSI as U
32 from esys.ripley import Brick, Rectangle
33 from .datasources import DataSource
34 from .coordinates import ReferenceSystem, CartesianReferenceSystem
35
36
37 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 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 """
49
50 def __init__(self, dim=3, reference_system=None):
51 """
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 :param reference_system: reference coordinate system. By default the
59 Cartesian coordinate system is used.
60 :type reference_system: `ReferenceSystem`
61 """
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 if not reference_system:
66 self.__reference_system=CartesianReferenceSystem()
67 else:
68 self.__reference_system=reference_system
69
70 if self.__reference_system.isCartesian():
71 self.__v_scale=1.
72 else:
73 self.__v_scale=1./self.getReferenceSystem().getHeightUnit()
74
75 self.__domain=None
76 self.__dim=dim
77 self.__sources=[]
78 self.__background_magnetic_field=None
79 self.__tags=[] # list of all tags being used by all data sources beeing attached:
80 self.setElementPadding()
81 self.setVerticalExtents()
82 self.fixDensityBelow()
83 self.fixSusceptibilityBelow()
84 self.fixVelocityBelow()
85
86 def getReferenceSystem(self):
87 """
88 returns the reference coordinate system
89
90 :rtype: `ReferenceSystem`
91 """
92 return self.__reference_system
93 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 def addSource(self, source):
101 """
102 Adds a survey data provider to the domain builder.
103 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 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
111 :param source: The data source to be added. Its reference system needs
112 to match the reference system of the DomainBuilder.
113 :type source: `DataSource`
114 """
115 if self.__domain is not None:
116 raise RuntimeError("Invalid call to addSource(). Domain is already built.")
117 if not isinstance(source, DataSource):
118 raise TypeError("source is not a DataSource")
119 if not source.getReferenceSystem() == self.getReferenceSystem():
120 raise ValueError("source reference system does not match.")
121
122 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 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 self.__sources.append(source)
130 if source.getTags(): self.__tags=list(set(self.__tags + source.getTags()))
131
132 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
133 """
134 Sets the amount of padding around the dataset as a fraction of the
135 dataset side lengths.
136
137 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 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
140
141 :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 :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 """
151 if not pad_lat == None:
152 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 if not pad_lon == None:
157 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 if self.__domain is not None:
162 raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
163 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
175 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
176 """
177 Sets the amount of padding around the dataset in absolute length units.
178
179 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 :type pad_x: ``float`` in units of length (meter)
185 :param pad_y: Padding per side in y direction (default: no padding)
186 :type pad_y: ``float`` in units of length (meter)
187 :note: `pad_y` is ignored for 2-dimensional domains.
188 :note: this function can only be used if the reference system is Cartesian
189 """
190 if not self.getReferenceSystem().isCartesian():
191 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
192 if self.__domain is not None:
193 raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
194 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
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
206 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
210 :param pad_lat: Padding per side in latitudinal direction (default: 0)
211 :type pad_lat: ``float`` in units of degree
212 :param pad_lon: Padding per side in longitudinal direction (default: 0)
213 :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 """
217 if self.getReferenceSystem().isCartesian():
218 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
219 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 Sets the amount of padding around the dataset in number of elements
232 (cells).
233
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 if not pad_lat == None:
245 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 if not pad_lon == None:
250 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 if self.__domain is not None:
256 raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
257 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
269 def getGravitySurveys(self):
270 """
271 Returns a list of gravity surveys, see `getSurveys` for details.
272 """
273 return self.getSurveys(DataSource.GRAVITY)
274
275 def getMagneticSurveys(self):
276 """
277 Returns a list of magnetic surveys, see `getSurveys` for details.
278 """
279 return self.getSurveys(DataSource.MAGNETIC)
280
281
282 def fixDensityBelow(self, depth=None):
283 """
284 Defines the depth below which the density anomaly is set to a given
285 value. If no value is given zero is assumed.
286
287 :param depth: depth below which the density is fixed. If not set, no
288 constraint at depth is applied.
289 :type depth: ``float``
290 """
291 self.__fix_density_below=depth
292
293 def fixSusceptibilityBelow(self, depth=None):
294 """
295 Defines the depth below which the susceptibility anomaly is set to a
296 given value. If no value is given zero is assumed.
297
298 :param depth: depth below which the susceptibility is fixed. If not
299 set, no constraint at depth is applied.
300 :type depth: ``float``
301 """
302 self.__fix_susceptibility_below=depth
303
304 def fixVelocityBelow(self, depth=None):
305 """
306 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 Returns a list of `Data` objects for all surveys of type `datatype`
319 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
322 :return: List of surveys which are tuples (anomaly,error).
323 :rtype: ``list``
324 """
325 surveys=[]
326 for src in self.__sources:
327 if src.getDataType()==datatype:
328 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 return surveys
331
332 def setBackgroundMagneticFluxDensity(self, B):
333 """
334 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
335 """
336 self.__background_magnetic_field=B
337
338 def getBackgroundMagneticFluxDensity(self):
339 """
340 Returns the background magnetic flux density.
341 """
342 B = self.__background_magnetic_field
343 if B is None:
344 raise ValueError("No background magnetic flux density set!")
345
346 if self.__dim < 3 :
347 return np.array([B[0], B[2]])
348 else:
349 return np.array(B)
350
351 def getSetDensityMask(self):
352 """
353 Returns the density mask data object which is non-zero for cells
354 whose density value is fixed, zero otherwise.
355 """
356 z=self.getDomain().getX()[self.__dim-1]
357 m = whereNonNegative(z)
358 if self.__fix_density_below:
359 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
360 return m
361
362 def getSetSusceptibilityMask(self):
363 """
364 Returns the susceptibility mask data object which is non-zero for
365 cells whose susceptibility value is fixed, zero otherwise.
366 """
367 z=self.getDomain().getX()[self.__dim-1]
368 m = whereNonNegative(z)
369 if self.__fix_susceptibility_below:
370 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
371 return m
372
373 def getDomain(self):
374 """
375 Returns a domain that spans the data area plus padding.
376
377 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 :rtype: `esys.escript.Domain`
383 """
384 if self.__domain is None:
385 self.__domain=self.__createDomain()
386 return self.__domain
387
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 :param depth: Depth of the domain (in meters)
394 :type depth: ``float``
395 :param air_layer: Depth of the layer above sea level (in meters)
396 :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 if self.__domain is not None:
402 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
403 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 Helper method that computes origin and number of data elements
410 after adding padding to the bounding box of all available survey data.
411 """
412 X0, NX, DX = self.__getTotalExtents()
413 DATA_DIM=len(X0)
414 frac=[]
415 # padding is applied to each side so multiply by 2 to get the total
416 # amount of padding per dimension
417 pad, pt = self._padding
418 for i in range(DATA_DIM):
419 if pad[i] is None:
420 frac.append(0.)
421 continue
422 if pt == 'f' : # fraction of side length
423 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
430 # calculate new number of elements
431 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 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 if len(self.__sources)==0:
442 raise ValueError("No data")
443 X0, NE, DX = self.__sources[0].getDataExtents()
444 # 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
448 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 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 origin = X0 + [-self._v_depth*self.__v_scale]
476
477 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
484 # cell size / point spacing
485 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
486 #self._spacing = [float(np.floor(si)) for si in spacing]
487 self._spacing = spacing
488
489 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
490
491 if self.__dim==3:
492 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
493 else:
494 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
495
496 # ripley may internally adjust NE and length, so recompute
497 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
498 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
499
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