/[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 4394 - (show annotations)
Tue May 7 04:56:59 2013 UTC (6 years, 3 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
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by 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 since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 """Domain construction from survey data for inversions"""
17
18 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 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 from esys.escript import unitsSI as U
31 from esys.ripley import Brick, Rectangle
32 from .datasources import DataSource
33 from .coordinates import ReferenceSystem, CartesianReferenceSystem
34
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
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 """
47
48 def __init__(self, dim=3, reference_system=None):
49 """
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 :param reference_system: reference coordinate system. By default the
57 Cartesian coordinate system is used.
58 :type reference_system: `ReferenceSystem`
59 """
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 if not reference_system:
64 self.__reference_system=CartesianReferenceSystem()
65 else:
66 self.__reference_system=reference_system
67 self.__domain=None
68 self.__dim=dim
69 self.__sources=[]
70 self.__background_magnetic_field=None
71 self.setElementPadding()
72 self.setVerticalExtents()
73 self.fixDensityBelow()
74 self.fixSusceptibilityBelow()
75
76 def getReferenceSystem(self):
77 """
78 returns the reference coordinate system
79
80 :rtype: `ReferenceSystem`
81 """
82 return self.__reference_system
83
84 def addSource(self, source):
85 """
86 Adds a survey data provider to the domain builder.
87 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 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
95 :param source: The data source to be added. Its reference system needs
96 to match the reference system of the DomainBuilder.
97 :type source: `DataSource`
98 """
99 if self.__domain is not None:
100 raise RuntimeError("Invalid call to addSource(). Domain is already built.")
101 if not isinstance(source, DataSource):
102 raise TypeError("source is not a DataSource")
103 if not source.getReferenceSystem() == self.getReferenceSystem():
104 raise ValueError("source reference system does not match.")
105
106 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 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 self.__sources.append(source)
114
115 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
116 """
117 Sets the amount of padding around the dataset as a fraction of the
118 dataset side lengths.
119
120 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 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
123
124 :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 :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 """
134 if not pad_lat == None:
135 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 if not pad_lon == None:
140 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 if self.__domain is not None:
145 raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
146 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
158 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
159 """
160 Sets the amount of padding around the dataset in absolute length units.
161
162 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 :type pad_x: ``float`` in units of length (meter)
168 :param pad_y: Padding per side in y direction (default: no padding)
169 :type pad_y: ``float`` in units of length (meter)
170 :note: `pad_y` is ignored for 2-dimensional domains.
171 :note: this function can only be used if the reference system is Cartesian
172 """
173 if not self.getReferenceSystem().isCartesian():
174 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
175 if self.__domain is not None:
176 raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
177 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
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
189 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 """
200 if self.getReferenceSystem().isCartesian():
201 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
202 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 Sets the amount of padding around the dataset in number of elements
215 (cells).
216
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 if not pad_lat == None:
228 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 if not pad_lon == None:
233 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 if self.__domain is not None:
239 raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
240 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
252 def getGravitySurveys(self):
253 """
254 Returns a list of gravity surveys, see `getSurveys` for details.
255 """
256 return self.getSurveys(DataSource.GRAVITY)
257
258 def getMagneticSurveys(self):
259 """
260 Returns a list of magnetic surveys, see `getSurveys` for details.
261 """
262 return self.getSurveys(DataSource.MAGNETIC)
263
264 def fixDensityBelow(self, depth=None):
265 """
266 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 """
273 self.__fix_density_below=depth
274
275 def fixSusceptibilityBelow(self, depth=None):
276 """
277 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 """
284 self.__fix_susceptibility_below=depth
285
286 def getSurveys(self, datatype):
287 """
288 Returns a list of `Data` objects for all surveys of type `datatype`
289 available to this domain builder.
290
291 :return: List of surveys which are tuples (anomaly,error).
292 :rtype: ``list``
293 """
294 surveys=[]
295 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 return surveys
299
300 def setBackgroundMagneticFluxDensity(self, B):
301 """
302 Sets the background magnetic flux density B=(B_North, B_East, B_Vertical)
303 """
304 self.__background_magnetic_field=B
305
306 def getBackgroundMagneticFluxDensity(self):
307 """
308 Returns the background magnetic flux density.
309 """
310 B = self.__background_magnetic_field
311 if self.__dim<3:
312 return np.array([B[0], B[2]])
313 else:
314 return np.array(B)
315
316 def getSetDensityMask(self):
317 """
318 Returns the density mask data object which is non-zero for cells
319 whose density value is fixed, zero otherwise.
320 """
321 z=self.getDomain().getX()[self.__dim-1]
322 m = whereNonNegative(z)
323 if self.__fix_density_below:
324 m += whereNonPositive(z+self.__fix_density_below)
325 return m
326
327 def getSetSusceptibilityMask(self):
328 """
329 Returns the susceptibility mask data object which is non-zero for
330 cells whose susceptibility value is fixed, zero otherwise.
331 """
332 z=self.getDomain().getX()[self.__dim-1]
333 m = whereNonNegative(z)
334 if self.__fix_susceptibility_below:
335 m += whereNonPositive(z+self.__fix_susceptibility_below)
336 return m
337
338 def getDomain(self):
339 """
340 Returns a domain that spans the data area plus padding.
341
342 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 :rtype: `esys.escript.Domain`
348 """
349 if self.__domain is None:
350 self.__domain=self.__createDomain()
351 return self.__domain
352
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 :param depth: Depth of the domain (in meters)
359 :type depth: ``float``
360 :param air_layer: Depth of the layer above sea level (in meters)
361 :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 if self.__domain is not None:
367 raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
368 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 Helper method that computes origin and number of data elements
375 after adding padding to the bounding box of all available survey data.
376 """
377 X0, NX, DX = self.__getTotalExtents()
378 DATA_DIM=len(X0)
379 frac=[]
380 # padding is applied to each side so multiply by 2 to get the total
381 # amount of padding per dimension
382 pad, pt = self._padding
383 for i in range(DATA_DIM):
384 if pad[i] is None:
385 frac.append(0.)
386 continue
387 if pt == 'f' : # fraction of side length
388 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
395 # calculate new number of elements
396 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 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 if len(self.__sources)==0:
407 raise ValueError("No data")
408 X0, NE, DX = self.__sources[0].getDataExtents()
409 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
410
411 for src in self.__sources[1:]:
412 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 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
445 # cell size / point spacing
446 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 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 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
456 # ripley may internally adjust NE and length, so recompute
457 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
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