/[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 4992 - (hide 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 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 caltinay 4060 from esys.ripley import Brick, Rectangle
33     from .datasources import DataSource
34 gross 4373 from .coordinates import ReferenceSystem, CartesianReferenceSystem
35 caltinay 4060
36 gross 4729
37 caltinay 4060 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 caltinay 4435 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 caltinay 4060 """
49    
50 gross 4373 def __init__(self, dim=3, reference_system=None):
51 caltinay 4060 """
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 gross 4373 :param reference_system: reference coordinate system. By default the
59     Cartesian coordinate system is used.
60     :type reference_system: `ReferenceSystem`
61 caltinay 4060 """
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 caltinay 4394 if not reference_system:
66     self.__reference_system=CartesianReferenceSystem()
67     else:
68     self.__reference_system=reference_system
69 caltinay 4435
70 caltinay 4436 if self.__reference_system.isCartesian():
71     self.__v_scale=1.
72     else:
73     self.__v_scale=1./self.getReferenceSystem().getHeightUnit()
74    
75 caltinay 4145 self.__domain=None
76     self.__dim=dim
77     self.__sources=[]
78     self.__background_magnetic_field=None
79 gross 4729 self.__tags=[] # list of all tags being used by all data sources beeing attached:
80 gross 4373 self.setElementPadding()
81 caltinay 4060 self.setVerticalExtents()
82 gross 4120 self.fixDensityBelow()
83     self.fixSusceptibilityBelow()
84 gross 4729 self.fixVelocityBelow()
85 gross 4373
86     def getReferenceSystem(self):
87     """
88     returns the reference coordinate system
89    
90     :rtype: `ReferenceSystem`
91     """
92     return self.__reference_system
93 gross 4729 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 caltinay 4060 def addSource(self, source):
101     """
102     Adds a survey data provider to the domain builder.
103 caltinay 4992 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 caltinay 4060
113 gross 4373 :param source: The data source to be added. Its reference system needs
114     to match the reference system of the DomainBuilder.
115 caltinay 4060 :type source: `DataSource`
116     """
117 caltinay 4150 if self.__domain is not None:
118     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
119 caltinay 4060 if not isinstance(source, DataSource):
120     raise TypeError("source is not a DataSource")
121 caltinay 4394 if not source.getReferenceSystem() == self.getReferenceSystem():
122     raise ValueError("source reference system does not match.")
123 caltinay 4201
124 caltinay 4362 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 caltinay 4992 if len(self.__sources)>0 and self.getReferenceSystem().isCartesian():
128 caltinay 4201 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 caltinay 4145 self.__sources.append(source)
132 gross 4729 if source.getTags(): self.__tags=list(set(self.__tags + source.getTags()))
133 caltinay 4060
134 gross 4373 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
135 caltinay 4060 """
136 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
137     dataset side lengths.
138    
139 caltinay 4150 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 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
142 caltinay 4060
143 caltinay 4108 :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 gross 4373 :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 caltinay 4060 """
153 gross 4373 if not pad_lat == None:
154 caltinay 4394 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 gross 4373 if not pad_lon == None:
159 caltinay 4394 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 caltinay 4150 if self.__domain is not None:
164     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
165 caltinay 4108 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 caltinay 4213
177 gross 4376 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
178 caltinay 4108 """
179     Sets the amount of padding around the dataset in absolute length units.
180 caltinay 4060
181 caltinay 4108 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 gross 4373 :type pad_x: ``float`` in units of length (meter)
187 caltinay 4108 :param pad_y: Padding per side in y direction (default: no padding)
188 gross 4373 :type pad_y: ``float`` in units of length (meter)
189 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
190 gross 4373 :note: this function can only be used if the reference system is Cartesian
191 caltinay 4108 """
192 gross 4373 if not self.getReferenceSystem().isCartesian():
193 caltinay 4394 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
194 caltinay 4150 if self.__domain is not None:
195     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
196 caltinay 4108 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 gross 4373
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 caltinay 4108
208 caltinay 4435 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 gross 4373
212 caltinay 4435 :param pad_lat: Padding per side in latitudinal direction (default: 0)
213 gross 4373 :type pad_lat: ``float`` in units of degree
214 caltinay 4435 :param pad_lon: Padding per side in longitudinal direction (default: 0)
215 gross 4373 :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 caltinay 4108 """
219 gross 4373 if self.getReferenceSystem().isCartesian():
220 caltinay 4394 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
221 gross 4373 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 caltinay 4145 Sets the amount of padding around the dataset in number of elements
234     (cells).
235 caltinay 4108
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 gross 4373 if not pad_lat == None:
247 caltinay 4394 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 gross 4373 if not pad_lon == None:
252 caltinay 4394 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 caltinay 4150 if self.__domain is not None:
258     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
259 caltinay 4108 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 caltinay 4213
271 caltinay 4060 def getGravitySurveys(self):
272 caltinay 4145 """
273     Returns a list of gravity surveys, see `getSurveys` for details.
274     """
275     return self.getSurveys(DataSource.GRAVITY)
276 caltinay 4213
277 gross 4115 def getMagneticSurveys(self):
278 gross 4120 """
279 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
280 gross 4120 """
281 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
282 caltinay 4213
283 gross 4729
284 caltinay 4145 def fixDensityBelow(self, depth=None):
285     """
286 caltinay 4435 Defines the depth below which the density anomaly is set to a given
287     value. If no value is given zero is assumed.
288 gross 4343
289 caltinay 4435 :param depth: depth below which the density is fixed. If not set, no
290     constraint at depth is applied.
291 gross 4343 :type depth: ``float``
292 caltinay 4145 """
293 gross 4120 self.__fix_density_below=depth
294 caltinay 4213
295 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
296 gross 4120 """
297 caltinay 4435 Defines the depth below which the susceptibility anomaly is set to a
298     given value. If no value is given zero is assumed.
299 gross 4343
300 caltinay 4435 :param depth: depth below which the susceptibility is fixed. If not
301     set, no constraint at depth is applied.
302 gross 4343 :type depth: ``float``
303 gross 4120 """
304     self.__fix_susceptibility_below=depth
305    
306 gross 4729 def fixVelocityBelow(self, depth=None):
307 caltinay 4060 """
308 gross 4729 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 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
321 gross 4729 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 caltinay 4060
324 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
325 caltinay 4060 :rtype: ``list``
326     """
327 gross 4124 surveys=[]
328 caltinay 4145 for src in self.__sources:
329     if src.getDataType()==datatype:
330 gross 4729 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 gross 4124 return surveys
333 caltinay 4060
334 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
335 caltinay 4060 """
336 gross 4433 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
337 caltinay 4060 """
338 gross 4106 self.__background_magnetic_field=B
339 caltinay 4108
340 gross 4115 def getBackgroundMagneticFluxDensity(self):
341 gross 4106 """
342 caltinay 4145 Returns the background magnetic flux density.
343 caltinay 4108 """
344 gross 4106 B = self.__background_magnetic_field
345 caltinay 4462 if B is None:
346     raise ValueError("No background magnetic flux density set!")
347    
348 gross 4433 if self.__dim < 3 :
349 caltinay 4287 return np.array([B[0], B[2]])
350 caltinay 4060 else:
351 caltinay 4287 return np.array(B)
352 caltinay 4060
353     def getSetDensityMask(self):
354 caltinay 4150 """
355     Returns the density mask data object which is non-zero for cells
356     whose density value is fixed, zero otherwise.
357     """
358 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
359 gross 4120 m = whereNonNegative(z)
360     if self.__fix_density_below:
361 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
362 gross 4120 return m
363 caltinay 4145
364 caltinay 4060 def getSetSusceptibilityMask(self):
365 caltinay 4150 """
366     Returns the susceptibility mask data object which is non-zero for
367     cells whose susceptibility value is fixed, zero otherwise.
368     """
369 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
370 gross 4120 m = whereNonNegative(z)
371     if self.__fix_susceptibility_below:
372 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
373 gross 4120 return m
374 caltinay 4060
375     def getDomain(self):
376     """
377     Returns a domain that spans the data area plus padding.
378    
379 caltinay 4145 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 caltinay 4060 :rtype: `esys.escript.Domain`
385     """
386 caltinay 4145 if self.__domain is None:
387     self.__domain=self.__createDomain()
388     return self.__domain
389 caltinay 4060
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 caltinay 4145 :param depth: Depth of the domain (in meters)
396 caltinay 4060 :type depth: ``float``
397 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
398 caltinay 4060 :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 caltinay 4150 if self.__domain is not None:
404     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
405 caltinay 4060 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 caltinay 4362 Helper method that computes origin and number of data elements
412 caltinay 4060 after adding padding to the bounding box of all available survey data.
413     """
414     X0, NX, DX = self.__getTotalExtents()
415 caltinay 4362 DATA_DIM=len(X0)
416 caltinay 4060 frac=[]
417     # padding is applied to each side so multiply by 2 to get the total
418     # amount of padding per dimension
419 caltinay 4108 pad, pt = self._padding
420 caltinay 4362 for i in range(DATA_DIM):
421 azadeh 4127 if pad[i] is None:
422 caltinay 4145 frac.append(0.)
423     continue
424 gross 4373 if pt == 'f' : # fraction of side length
425 caltinay 4108 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 caltinay 4060
432     # calculate new number of elements
433 caltinay 4362 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 caltinay 4060 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 caltinay 4145 if len(self.__sources)==0:
444 caltinay 4060 raise ValueError("No data")
445 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
446 caltinay 4613 # 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 caltinay 4060
450 caltinay 4613 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 caltinay 4060 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 caltinay 4435 origin = X0 + [-self._v_depth*self.__v_scale]
478    
479 caltinay 4394 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 caltinay 4435
486 caltinay 4060 # cell size / point spacing
487 caltinay 4435 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
488 gross 4373 #self._spacing = [float(np.floor(si)) for si in spacing]
489     self._spacing = spacing
490 gross 4433
491 caltinay 4435 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
492 caltinay 4613
493 caltinay 4145 if self.__dim==3:
494 caltinay 4060 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
495     else:
496 caltinay 4434 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
497 caltinay 4108
498 caltinay 4060 # ripley may internally adjust NE and length, so recompute
499 caltinay 4145 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
500 caltinay 4435 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
501 caltinay 4060
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