/[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 4995 - (show annotations)
Thu Jun 5 03:44:16 2014 UTC (5 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 21039 byte(s)
Doc fixes and ensuring that expensive debug ops will not be processed if debug
mode is disabled.

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 # list of all tags used by all data sources being attached
80 self.__tags=[]
81 self.setElementPadding()
82 self.setVerticalExtents()
83 self.fixDensityBelow()
84 self.fixSusceptibilityBelow()
85 self.fixVelocityBelow()
86
87 def getReferenceSystem(self):
88 """
89 returns the reference coordinate system
90
91 :rtype: `ReferenceSystem`
92 """
93 return self.__reference_system
94
95 def getTags(self):
96 """
97 returns a list of all tags in use by the attached data sources.
98 The list may be empty.
99 """
100 return self.__tags
101
102 def addSource(self, source):
103 """
104 Adds a survey data provider to the domain builder.
105 An exception is raised if the domain has already been built.
106 An exception is also reported if the reference system used is
107 cartesian and the UTM zone of `source` does not match the UTM zone of
108 sources already added to the domain builder (see Inversion Cookbook
109 for more information).
110 The dimensionality of the data source must be compatible with this
111 domain builder. That is, the dimensionality of the data must be one
112 less than the dimensionality of the domain (specified in the
113 constructor).
114
115 :param source: The data source to be added. Its reference system needs
116 to match the reference system of the DomainBuilder.
117 :type source: `DataSource`
118 """
119 if self.__domain is not None:
120 raise RuntimeError("Invalid call to addSource(). Domain is already built.")
121 if not isinstance(source, DataSource):
122 raise TypeError("source is not a DataSource")
123 if not source.getReferenceSystem() == self.getReferenceSystem():
124 raise ValueError("source reference system does not match.")
125
126 DATA_DIM = len(source.getDataExtents()[0])
127 if DATA_DIM != self.__dim-1:
128 raise ValueError("Data must be %d-dimensional."%(self.__dim-1))
129 if len(self.__sources)>0 and self.getReferenceSystem().isCartesian():
130 if self.__sources[0].getUtmZone() != source.getUtmZone():
131 raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")
132
133 self.__sources.append(source)
134 if source.getTags(): self.__tags=list(set(self.__tags + source.getTags()))
135
136 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
137 """
138 Sets the amount of padding around the dataset as a fraction of the
139 dataset side lengths.
140
141 For example, calling ``setFractionalPadding(0.2, 0.1)`` with a data
142 source of size 10x20 will result in the padded data set size
143 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
144
145 :param pad_x: Padding per side in x direction (default: no padding)
146 :type pad_x: ``float``
147 :param pad_y: Padding per side in y direction (default: no padding)
148 :type pad_y: ``float``
149 :param pad_lat: Padding per side in latitudinal direction (default: no padding)
150 :type pad_lat: ``float``
151 :param pad_lon: Padding per side in longitudinal direction (default: no padding)
152 :type pad_lon: ``float``
153 :note: `pad_y` is ignored for 2-dimensional domains.
154 """
155 if not pad_lat == None:
156 if not pad_x == None:
157 raise ValueError("Either pad_lat or pad_x can be set.")
158 else:
159 pad_x = pad_lat
160 if not pad_lon == None:
161 if not pad_y == None:
162 raise ValueError("Either pad_lon or pad_y can be set.")
163 else:
164 pad_y = pad_lan
165 if self.__domain is not None:
166 raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
167 if pad_x is not None:
168 if pad_x < 0:
169 raise ValueError("setFractionalPadding: Arguments must be non-negative")
170 if pad_x > 10:
171 raise ValueError("setFractionalPadding: Argument too large")
172 if pad_y is not None:
173 if pad_y < 0:
174 raise ValueError("setFractionalPadding: Arguments must be non-negative")
175 if pad_y > 10:
176 raise ValueError("setFractionalPadding: Argument too large")
177 self._padding = [pad_x,pad_y], 'f'
178
179 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
180 """
181 Sets the amount of padding around the dataset in absolute length units.
182
183 The final domain size will be the length in x (in y) of the dataset
184 plus twice the value of `pad_x` (`pad_y`). The arguments must be
185 non-negative.
186
187 :param pad_x: Padding per side in x direction (default: no padding)
188 :type pad_x: ``float`` in units of length (meter)
189 :param pad_y: Padding per side in y direction (default: no padding)
190 :type pad_y: ``float`` in units of length (meter)
191 :note: `pad_y` is ignored for 2-dimensional domains.
192 :note: this function can only be used if the reference system is Cartesian
193 """
194 if not self.getReferenceSystem().isCartesian():
195 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
196 if self.__domain is not None:
197 raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
198 if pad_x is not None:
199 if pad_x < 0:
200 raise ValueError("setPadding: Arguments must be non-negative")
201 if pad_y is not None:
202 if pad_y < 0:
203 raise ValueError("setPadding: Arguments must be non-negative")
204 self._padding = [pad_x,pad_y], 'l'
205
206 def setGeoPadding(self, pad_lat=None, pad_lon=None):
207 """
208 Sets the amount of padding around the dataset in longitude and latitude.
209
210 The final domain size will be the extent in the latitudinal (in
211 longitudinal) direction of the dataset plus twice the value of
212 `pad_lat` (`pad_lon`). The arguments must be non-negative.
213
214 :param pad_lat: Padding per side in latitudinal direction (default: 0)
215 :type pad_lat: ``float`` in units of degree
216 :param pad_lon: Padding per side in longitudinal direction (default: 0)
217 :type pad_lon: ``float`` in units of degree
218 :note: `pad_lon` is ignored for 2-dimensional domains.
219 :note: this function can only be used if the reference system is not Cartesian
220 """
221 if self.getReferenceSystem().isCartesian():
222 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
223 if self.__domain is not None:
224 raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
225 if pad_lat is not None:
226 if pad_lat < 0:
227 raise ValueError("setPadding: Arguments must be non-negative")
228 if pad_lon is not None:
229 if pad_lon < 0:
230 raise ValueError("setPadding: Arguments must be non-negative")
231 self._padding = [pad_lat,pad_lon], 'd'
232
233 def setElementPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
234 """
235 Sets the amount of padding around the dataset in number of elements
236 (cells).
237
238 When the domain is constructed `pad_x` (`pad_y`) elements are added
239 on each side of the x- (y-) dimension. The arguments must be
240 non-negative.
241
242 :param pad_x: Padding per side in x direction (default: no padding)
243 :type pad_x: ``int``
244 :param pad_y: Padding per side in y direction (default: no padding)
245 :type pad_y: ``int``
246 :note: `pad_y` is ignored for 2-dimensional datasets.
247 """
248 if not pad_lat == None:
249 if not pad_x == None:
250 raise ValueError("Either pad_lat or pad_x can be set.")
251 else:
252 pad_x = pad_lat
253 if not pad_lon == None:
254 if not pad_y == None:
255 raise ValueError("Either pad_lon or pad_y can be set.")
256 else:
257 pad_y = pad_lan
258
259 if self.__domain is not None:
260 raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
261 if pad_x is not None:
262 if type(pad_x) is not int:
263 raise TypeError("setElementPadding expects integer arguments")
264 if pad_x < 0:
265 raise ValueError("setElementPadding: Arguments must be non-negative")
266 if pad_y is not None:
267 if type(pad_y) is not int:
268 raise TypeError("setElementPadding expects integer arguments")
269 if pad_y < 0:
270 raise ValueError("setElementPadding: Arguments must be non-negative")
271 self._padding = [pad_x,pad_y], 'e'
272
273 def getGravitySurveys(self):
274 """
275 Returns a list of gravity surveys, see `getSurveys` for details.
276 """
277 return self.getSurveys(DataSource.GRAVITY)
278
279 def getMagneticSurveys(self):
280 """
281 Returns a list of magnetic surveys, see `getSurveys` for details.
282 """
283 return self.getSurveys(DataSource.MAGNETIC)
284
285
286 def fixDensityBelow(self, depth=None):
287 """
288 Defines the depth below which the density anomaly is set to a given
289 value. If no value is given zero is assumed.
290
291 :param depth: depth below which the density is fixed. If not set, no
292 constraint at depth is applied.
293 :type depth: ``float``
294 """
295 self.__fix_density_below=depth
296
297 def fixSusceptibilityBelow(self, depth=None):
298 """
299 Defines the depth below which the susceptibility anomaly is set to a
300 given value. If no value is given zero is assumed.
301
302 :param depth: depth below which the susceptibility is fixed. If not
303 set, no constraint at depth is applied.
304 :type depth: ``float``
305 """
306 self.__fix_susceptibility_below=depth
307
308 def fixVelocityBelow(self, depth=None):
309 """
310 Defines the depth below which the velocity and Q index is set to a
311 given value. If no value is given zero is assumed.
312
313 :param depth: depth below which the velocity is fixed. If not
314 set, no constraint at depth is applied.
315 :type depth: ``float``
316 """
317 self.__fix_velocity_below=depth
318
319
320 def getSurveys(self, datatype, tags=None):
321 """
322 Returns a list of `Data` objects for all surveys of type `datatype`
323 available to this domain builder. If a list of `tags` is given
324 only data sources whose tag matches the tag list are returned.
325
326 :return: List of surveys which are tuples (anomaly,error).
327 :rtype: ``list``
328 """
329 surveys=[]
330 for src in self.__sources:
331 if src.getDataType()==datatype:
332 if tags is None or ( src.getTags() is not None and all( [ t in tags for t in src.getTags() ] ) ) :
333 surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
334 return surveys
335
336 def setBackgroundMagneticFluxDensity(self, B):
337 """
338 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
339 """
340 self.__background_magnetic_field=B
341
342 def getBackgroundMagneticFluxDensity(self):
343 """
344 Returns the background magnetic flux density.
345 """
346 B = self.__background_magnetic_field
347 if B is None:
348 raise ValueError("No background magnetic flux density set!")
349
350 if self.__dim < 3 :
351 return np.array([B[0], B[2]])
352 else:
353 return np.array(B)
354
355 def getSetDensityMask(self):
356 """
357 Returns the density mask data object which is non-zero for cells
358 whose density value is fixed, zero otherwise.
359 """
360 z=self.getDomain().getX()[self.__dim-1]
361 m = whereNonNegative(z)
362 if self.__fix_density_below:
363 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
364 return m
365
366 def getSetSusceptibilityMask(self):
367 """
368 Returns the susceptibility mask data object which is non-zero for
369 cells whose susceptibility value is fixed, zero otherwise.
370 """
371 z=self.getDomain().getX()[self.__dim-1]
372 m = whereNonNegative(z)
373 if self.__fix_susceptibility_below:
374 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
375 return m
376
377 def getDomain(self):
378 """
379 Returns a domain that spans the data area plus padding.
380
381 The domain is created the first time this method is called,
382 subsequent calls return the same domain so anything that affects
383 the domain (such as padding) needs to be set beforehand.
384
385 :return: The escript domain for this data source
386 :rtype: `esys.escript.Domain`
387 """
388 if self.__domain is None:
389 self.__domain=self.__createDomain()
390 return self.__domain
391
392 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
393 """
394 This method sets the target domain parameters for the vertical
395 dimension.
396
397 :param depth: Depth of the domain (in meters)
398 :type depth: ``float``
399 :param air_layer: Depth of the layer above sea level (in meters)
400 :type air_layer: ``float``
401 :param num_cells: Number of domain elements for the entire vertical
402 dimension
403 :type num_cells: ``int``
404 """
405 if self.__domain is not None:
406 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
407 self._v_depth=depth
408 self._v_air_layer=air_layer
409 self._v_num_cells=num_cells
410
411 def __getTotalExtentsWithPadding(self):
412 """
413 Helper method that computes origin and number of data elements
414 after adding padding to the bounding box of all available survey data.
415 """
416 X0, NX, DX = self.__getTotalExtents()
417 DATA_DIM=len(X0)
418 frac=[]
419 # padding is applied to each side so multiply by 2 to get the total
420 # amount of padding per dimension
421 pad, pt = self._padding
422 for i in range(DATA_DIM):
423 if pad[i] is None:
424 frac.append(0.)
425 continue
426 if pt == 'f' : # fraction of side length
427 frac.append(2.*pad[i])
428 elif pt == 'e': # number of elements
429 frac.append(2.*pad[i]/float(NX[i]))
430 else: # absolute length
431 f=pad[i]/DX[i]
432 frac.append(2.*f/float(NX[i]))
433
434 # calculate new number of elements
435 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DATA_DIM)]
436 NXdiff=[NX_padded[i]-NX[i] for i in range(DATA_DIM)]
437 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DATA_DIM)]
438 return X0_padded, NX_padded, DX
439
440 def __getTotalExtents(self):
441 """
442 Helper method that computes the origin, number of elements and
443 minimal element spacing taking into account all available survey data.
444 """
445 if len(self.__sources)==0:
446 raise ValueError("No data")
447 X0, NE, DX = self.__sources[0].getDataExtents()
448 # do not mess with the values if only one source used
449 if len(self.__sources)>1:
450 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
451
452 for src in self.__sources[1:]:
453 d_x0, d_ne, d_dx = src.getDataExtents()
454 for i in range(len(d_x0)):
455 X0[i]=min(X0[i], d_x0[i])
456 for i in range(len(d_dx)):
457 DX[i]=min(DX[i], d_dx[i])
458 for i in range(len(d_ne)):
459 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
460 # FIXME: should this be rounded up instead?
461 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
462 return X0, NE, DX
463
464 def __createDomain(self):
465 """
466 Creates and returns an escript domain that spans the entire area of
467 available data plus a padding zone. This method is called only once
468 the first time `getDomain()` is invoked.
469
470 :return: The escript domain
471 :rtype: `esys.escript.Domain`
472 """
473 X0, NX, DX = self.__getTotalExtentsWithPadding()
474
475 # number of domain elements
476 NE = NX + [self._v_num_cells]
477
478 # origin of domain
479 origin = X0 + [-self._v_depth*self.__v_scale]
480
481 if self.getReferenceSystem().isCartesian():
482 # rounding will give us about meter-accuracy with UTM coordinates
483 self._dom_origin = [np.floor(oi) for oi in origin]
484 else:
485 # this should give us about meter-accuracy with lat/lon coords
486 self._dom_origin = [1e-5*np.floor(oi*1e5) for oi in origin]
487
488 # cell size / point spacing
489 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
490 #self._spacing = [float(np.floor(si)) for si in spacing]
491 self._spacing = spacing
492
493 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
494
495 if self.__dim==3:
496 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
497 else:
498 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
499
500 # ripley may internally adjust NE and length, so recompute
501 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
502 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
503
504 self.logger.debug("Domain size: "+str(self._dom_NE))
505 self.logger.debug(" length: "+str(self._dom_len))
506 self.logger.debug(" origin: "+str(self._dom_origin))
507 return dom
508

  ViewVC Help
Powered by ViewVC 1.1.26