/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 19890 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26