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

  ViewVC Help
Powered by ViewVC 1.1.26