/[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 4992 - (show annotations)
Thu Jun 5 01:38:33 2014 UTC (5 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 21046 byte(s)
Don't obtain/report UTM zone for non-cartesian coordinate systems.

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

  ViewVC Help
Powered by ViewVC 1.1.26