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

  ViewVC Help
Powered by ViewVC 1.1.26