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

  ViewVC Help
Powered by ViewVC 1.1.26