/[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 4462 - (show annotations)
Mon Jun 17 04:12:06 2013 UTC (6 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 19638 byte(s)
Raise exception if background magnetic field has not been set but is required.

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

  ViewVC Help
Powered by ViewVC 1.1.26