/[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 4436 - (show annotations)
Mon Jun 3 05:37:31 2013 UTC (6 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 19542 byte(s)
Tweak.

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 self.__dim < 3 :
318 return np.array([B[0], B[2]])
319 else:
320 return np.array(B)
321
322 def getSetDensityMask(self):
323 """
324 Returns the density mask data object which is non-zero for cells
325 whose density value is fixed, zero otherwise.
326 """
327 z=self.getDomain().getX()[self.__dim-1]
328 m = whereNonNegative(z)
329 if self.__fix_density_below:
330 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
331 return m
332
333 def getSetSusceptibilityMask(self):
334 """
335 Returns the susceptibility mask data object which is non-zero for
336 cells whose susceptibility value is fixed, zero otherwise.
337 """
338 z=self.getDomain().getX()[self.__dim-1]
339 m = whereNonNegative(z)
340 if self.__fix_susceptibility_below:
341 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
342 return m
343
344 def getDomain(self):
345 """
346 Returns a domain that spans the data area plus padding.
347
348 The domain is created the first time this method is called,
349 subsequent calls return the same domain so anything that affects
350 the domain (such as padding) needs to be set beforehand.
351
352 :return: The escript domain for this data source
353 :rtype: `esys.escript.Domain`
354 """
355 if self.__domain is None:
356 self.__domain=self.__createDomain()
357 return self.__domain
358
359 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
360 """
361 This method sets the target domain parameters for the vertical
362 dimension.
363
364 :param depth: Depth of the domain (in meters)
365 :type depth: ``float``
366 :param air_layer: Depth of the layer above sea level (in meters)
367 :type air_layer: ``float``
368 :param num_cells: Number of domain elements for the entire vertical
369 dimension
370 :type num_cells: ``int``
371 """
372 if self.__domain is not None:
373 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
374 self._v_depth=depth
375 self._v_air_layer=air_layer
376 self._v_num_cells=num_cells
377
378 def __getTotalExtentsWithPadding(self):
379 """
380 Helper method that computes origin and number of data elements
381 after adding padding to the bounding box of all available survey data.
382 """
383 X0, NX, DX = self.__getTotalExtents()
384 DATA_DIM=len(X0)
385 frac=[]
386 # padding is applied to each side so multiply by 2 to get the total
387 # amount of padding per dimension
388 pad, pt = self._padding
389 for i in range(DATA_DIM):
390 if pad[i] is None:
391 frac.append(0.)
392 continue
393 if pt == 'f' : # fraction of side length
394 frac.append(2.*pad[i])
395 elif pt == 'e': # number of elements
396 frac.append(2.*pad[i]/float(NX[i]))
397 else: # absolute length
398 f=pad[i]/DX[i]
399 frac.append(2.*f/float(NX[i]))
400
401 # calculate new number of elements
402 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DATA_DIM)]
403 NXdiff=[NX_padded[i]-NX[i] for i in range(DATA_DIM)]
404 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DATA_DIM)]
405 return X0_padded, NX_padded, DX
406
407 def __getTotalExtents(self):
408 """
409 Helper method that computes the origin, number of elements and
410 minimal element spacing taking into account all available survey data.
411 """
412 if len(self.__sources)==0:
413 raise ValueError("No data")
414 X0, NE, DX = self.__sources[0].getDataExtents()
415 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
416
417 for src in self.__sources[1:]:
418 d_x0, d_ne, d_dx = src.getDataExtents()
419 for i in range(len(d_x0)):
420 X0[i]=min(X0[i], d_x0[i])
421 for i in range(len(d_dx)):
422 DX[i]=min(DX[i], d_dx[i])
423 for i in range(len(d_ne)):
424 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
425 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
426 return X0, NE, DX
427
428 def __createDomain(self):
429 """
430 Creates and returns an escript domain that spans the entire area of
431 available data plus a padding zone. This method is called only once
432 the first time `getDomain()` is invoked.
433
434 :return: The escript domain
435 :rtype: `esys.escript.Domain`
436 """
437 X0, NX, DX = self.__getTotalExtentsWithPadding()
438
439 # number of domain elements
440 NE = NX + [self._v_num_cells]
441
442 # origin of domain
443 origin = X0 + [-self._v_depth*self.__v_scale]
444
445 if self.getReferenceSystem().isCartesian():
446 # rounding will give us about meter-accuracy with UTM coordinates
447 self._dom_origin = [np.floor(oi) for oi in origin]
448 else:
449 # this should give us about meter-accuracy with lat/lon coords
450 self._dom_origin = [1e-5*np.floor(oi*1e5) for oi in origin]
451
452 # cell size / point spacing
453 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
454 #self._spacing = [float(np.floor(si)) for si in spacing]
455 self._spacing = spacing
456
457 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
458
459 if self.__dim==3:
460 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
461 else:
462 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
463
464 # ripley may internally adjust NE and length, so recompute
465 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
466 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
467
468 self.logger.debug("Domain size: "+str(self._dom_NE))
469 self.logger.debug(" length: "+str(self._dom_len))
470 self.logger.debug(" origin: "+str(self._dom_origin))
471 return dom
472

  ViewVC Help
Powered by ViewVC 1.1.26