/[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 4435 - (hide annotations)
Mon Jun 3 04:27:57 2013 UTC (6 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 19443 byte(s)
Fixed a few scaling problems for non-cartesian coordinates in domainbuilder.

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

  ViewVC Help
Powered by ViewVC 1.1.26