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

1 caltinay 4060
2     ##############################################################################
3     #
4 caltinay 4145 # Copyright (c) 2003-2013 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     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15    
16     """Domain construction from survey data for inversions"""
17    
18 caltinay 4145 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 caltinay 4060 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 caltinay 4097 from esys.escript import unitsSI as U
31 caltinay 4060 from esys.ripley import Brick, Rectangle
32     from .datasources import DataSource
33 gross 4373 from .coordinates import ReferenceSystem, CartesianReferenceSystem
34 caltinay 4060
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 caltinay 4435 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 caltinay 4060 """
47    
48 gross 4373 def __init__(self, dim=3, reference_system=None):
49 caltinay 4060 """
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 gross 4373 :param reference_system: reference coordinate system. By default the
57     Cartesian coordinate system is used.
58     :type reference_system: `ReferenceSystem`
59 caltinay 4060 """
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 caltinay 4394 if not reference_system:
64     self.__reference_system=CartesianReferenceSystem()
65     else:
66     self.__reference_system=reference_system
67 caltinay 4435
68 caltinay 4436 if self.__reference_system.isCartesian():
69     self.__v_scale=1.
70     else:
71     self.__v_scale=1./self.getReferenceSystem().getHeightUnit()
72    
73 caltinay 4145 self.__domain=None
74     self.__dim=dim
75     self.__sources=[]
76     self.__background_magnetic_field=None
77 gross 4373 self.setElementPadding()
78 caltinay 4060 self.setVerticalExtents()
79 gross 4120 self.fixDensityBelow()
80     self.fixSusceptibilityBelow()
81 gross 4373
82     def getReferenceSystem(self):
83     """
84     returns the reference coordinate system
85    
86     :rtype: `ReferenceSystem`
87     """
88     return self.__reference_system
89    
90 caltinay 4060 def addSource(self, source):
91     """
92     Adds a survey data provider to the domain builder.
93 caltinay 4201 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 caltinay 4362 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 caltinay 4060
101 gross 4373 :param source: The data source to be added. Its reference system needs
102     to match the reference system of the DomainBuilder.
103 caltinay 4060 :type source: `DataSource`
104     """
105 caltinay 4150 if self.__domain is not None:
106     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
107 caltinay 4060 if not isinstance(source, DataSource):
108     raise TypeError("source is not a DataSource")
109 caltinay 4394 if not source.getReferenceSystem() == self.getReferenceSystem():
110     raise ValueError("source reference system does not match.")
111 caltinay 4201
112 caltinay 4362 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 caltinay 4201 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 caltinay 4145 self.__sources.append(source)
120 caltinay 4060
121 gross 4373 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
122 caltinay 4060 """
123 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
124     dataset side lengths.
125    
126 caltinay 4150 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 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
129 caltinay 4060
130 caltinay 4108 :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 gross 4373 :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 caltinay 4060 """
140 gross 4373 if not pad_lat == None:
141 caltinay 4394 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 gross 4373 if not pad_lon == None:
146 caltinay 4394 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 caltinay 4150 if self.__domain is not None:
151     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
152 caltinay 4108 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 caltinay 4213
164 gross 4376 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
165 caltinay 4108 """
166     Sets the amount of padding around the dataset in absolute length units.
167 caltinay 4060
168 caltinay 4108 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 gross 4373 :type pad_x: ``float`` in units of length (meter)
174 caltinay 4108 :param pad_y: Padding per side in y direction (default: no padding)
175 gross 4373 :type pad_y: ``float`` in units of length (meter)
176 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
177 gross 4373 :note: this function can only be used if the reference system is Cartesian
178 caltinay 4108 """
179 gross 4373 if not self.getReferenceSystem().isCartesian():
180 caltinay 4394 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
181 caltinay 4150 if self.__domain is not None:
182     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
183 caltinay 4108 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 gross 4373
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 caltinay 4108
195 caltinay 4435 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 gross 4373
199 caltinay 4435 :param pad_lat: Padding per side in latitudinal direction (default: 0)
200 gross 4373 :type pad_lat: ``float`` in units of degree
201 caltinay 4435 :param pad_lon: Padding per side in longitudinal direction (default: 0)
202 gross 4373 :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 caltinay 4108 """
206 gross 4373 if self.getReferenceSystem().isCartesian():
207 caltinay 4394 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
208 gross 4373 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 caltinay 4145 Sets the amount of padding around the dataset in number of elements
221     (cells).
222 caltinay 4108
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 gross 4373 if not pad_lat == None:
234 caltinay 4394 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 gross 4373 if not pad_lon == None:
239 caltinay 4394 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 caltinay 4150 if self.__domain is not None:
245     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
246 caltinay 4108 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 caltinay 4213
258 caltinay 4060 def getGravitySurveys(self):
259 caltinay 4145 """
260     Returns a list of gravity surveys, see `getSurveys` for details.
261     """
262     return self.getSurveys(DataSource.GRAVITY)
263 caltinay 4213
264 gross 4115 def getMagneticSurveys(self):
265 gross 4120 """
266 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
267 gross 4120 """
268 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
269 caltinay 4213
270 caltinay 4145 def fixDensityBelow(self, depth=None):
271     """
272 caltinay 4435 Defines the depth below which the density anomaly is set to a given
273     value. If no value is given zero is assumed.
274 gross 4343
275 caltinay 4435 :param depth: depth below which the density is fixed. If not set, no
276     constraint at depth is applied.
277 gross 4343 :type depth: ``float``
278 caltinay 4145 """
279 gross 4120 self.__fix_density_below=depth
280 caltinay 4213
281 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
282 gross 4120 """
283 caltinay 4435 Defines the depth below which the susceptibility anomaly is set to a
284     given value. If no value is given zero is assumed.
285 gross 4343
286 caltinay 4435 :param depth: depth below which the susceptibility is fixed. If not
287     set, no constraint at depth is applied.
288 gross 4343 :type depth: ``float``
289 gross 4120 """
290     self.__fix_susceptibility_below=depth
291    
292 gross 4115 def getSurveys(self, datatype):
293 caltinay 4060 """
294 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
295     available to this domain builder.
296 caltinay 4060
297 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
298 caltinay 4060 :rtype: ``list``
299     """
300 gross 4124 surveys=[]
301 caltinay 4145 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 gross 4124 return surveys
305 caltinay 4060
306 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
307 caltinay 4060 """
308 gross 4433 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
309 caltinay 4060 """
310 gross 4106 self.__background_magnetic_field=B
311 caltinay 4108
312 gross 4115 def getBackgroundMagneticFluxDensity(self):
313 gross 4106 """
314 caltinay 4145 Returns the background magnetic flux density.
315 caltinay 4108 """
316 gross 4106 B = self.__background_magnetic_field
317 gross 4433 if self.__dim < 3 :
318 caltinay 4287 return np.array([B[0], B[2]])
319 caltinay 4060 else:
320 caltinay 4287 return np.array(B)
321 caltinay 4060
322     def getSetDensityMask(self):
323 caltinay 4150 """
324     Returns the density mask data object which is non-zero for cells
325     whose density value is fixed, zero otherwise.
326     """
327 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
328 gross 4120 m = whereNonNegative(z)
329     if self.__fix_density_below:
330 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
331 gross 4120 return m
332 caltinay 4145
333 caltinay 4060 def getSetSusceptibilityMask(self):
334 caltinay 4150 """
335     Returns the susceptibility mask data object which is non-zero for
336     cells whose susceptibility value is fixed, zero otherwise.
337     """
338 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
339 gross 4120 m = whereNonNegative(z)
340     if self.__fix_susceptibility_below:
341 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
342 gross 4120 return m
343 caltinay 4060
344     def getDomain(self):
345     """
346     Returns a domain that spans the data area plus padding.
347    
348 caltinay 4145 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 caltinay 4060 :rtype: `esys.escript.Domain`
354     """
355 caltinay 4145 if self.__domain is None:
356     self.__domain=self.__createDomain()
357     return self.__domain
358 caltinay 4060
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 caltinay 4145 :param depth: Depth of the domain (in meters)
365 caltinay 4060 :type depth: ``float``
366 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
367 caltinay 4060 :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 caltinay 4150 if self.__domain is not None:
373     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
374 caltinay 4060 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 caltinay 4362 Helper method that computes origin and number of data elements
381 caltinay 4060 after adding padding to the bounding box of all available survey data.
382     """
383     X0, NX, DX = self.__getTotalExtents()
384 caltinay 4362 DATA_DIM=len(X0)
385 caltinay 4060 frac=[]
386     # padding is applied to each side so multiply by 2 to get the total
387     # amount of padding per dimension
388 caltinay 4108 pad, pt = self._padding
389 caltinay 4362 for i in range(DATA_DIM):
390 azadeh 4127 if pad[i] is None:
391 caltinay 4145 frac.append(0.)
392     continue
393 gross 4373 if pt == 'f' : # fraction of side length
394 caltinay 4108 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 caltinay 4060
401     # calculate new number of elements
402 caltinay 4362 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 caltinay 4060 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 caltinay 4145 if len(self.__sources)==0:
413 caltinay 4060 raise ValueError("No data")
414 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
415 caltinay 4060 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
416    
417 caltinay 4145 for src in self.__sources[1:]:
418 caltinay 4060 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 caltinay 4435 origin = X0 + [-self._v_depth*self.__v_scale]
444    
445 caltinay 4394 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 caltinay 4435
452 caltinay 4060 # cell size / point spacing
453 caltinay 4435 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
454 gross 4373 #self._spacing = [float(np.floor(si)) for si in spacing]
455     self._spacing = spacing
456 gross 4433
457 caltinay 4435 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
458 gross 4433
459 caltinay 4145 if self.__dim==3:
460 caltinay 4060 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
461     else:
462 caltinay 4434 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
463 caltinay 4108
464 caltinay 4060 # ripley may internally adjust NE and length, so recompute
465 caltinay 4145 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
466 caltinay 4435 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
467 caltinay 4060
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