/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26