/[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 4394 - (hide annotations)
Tue May 7 04:56:59 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 19325 byte(s)
Fixed rounding with lat/lon coordinates and removed stray print statements.

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

  ViewVC Help
Powered by ViewVC 1.1.26