/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Diff of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4362 by caltinay, Tue Apr 16 04:37:13 2013 UTC revision 4729 by gross, Sat Mar 8 00:56:37 2014 UTC
# Line 1  Line 1 
1    
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2013 by University of Queensland  # Copyright (c) 2003-2014 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 9  Line 9 
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  # Development since 2012 by School of Earth Sciences  # 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"""  """Domain construction from survey data for inversions"""
18    
19  __copyright__="""Copyright (c) 2003-2013 by University of Queensland  __copyright__="""Copyright (c) 2003-2014 by University of Queensland
20  http://www.uq.edu.au  http://www.uq.edu.au
21  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
22  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 30  from esys.escript.util import * Line 31  from esys.escript.util import *
31  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
32  from esys.ripley import Brick, Rectangle  from esys.ripley import Brick, Rectangle
33  from .datasources import DataSource  from .datasources import DataSource
34    from .coordinates import ReferenceSystem, CartesianReferenceSystem
35    
36    
37  class DomainBuilder(object):  class DomainBuilder(object):
38      """      """
39      This class is responsible for constructing an escript Domain object with      This class is responsible for constructing an escript Domain object with
40      suitable extents and resolution for survey data (`DataSource` objects)      suitable extents and resolution for survey data (`DataSource` objects)
41      that is added to it.      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):      def __init__(self, dim=3, reference_system=None):
51          """          """
52          Constructor.          Constructor.
53    
# Line 46  class DomainBuilder(object): Line 55  class DomainBuilder(object):
55                      This has implications for the survey data than can be                      This has implications for the survey data than can be
56                      added. By default a 3D domain is created.                      added. By default a 3D domain is created.
57          :type dim: ``int``          :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__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
63          if dim not in (2,3):          if dim not in (2,3):
64              raise ValueError("Number of dimensions must be 2 or 3")              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          self.__domain=None
76          self.__dim=dim          self.__dim=dim
77          self.__sources=[]          self.__sources=[]
78          self.__background_magnetic_field=None          self.__background_magnetic_field=None
79          self.setPadding()          self.__tags=[]   # list of all tags being used by all data sources beeing attached:
80            self.setElementPadding()
81          self.setVerticalExtents()          self.setVerticalExtents()
82          self.fixDensityBelow()          self.fixDensityBelow()
83          self.fixSusceptibilityBelow()          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):      def addSource(self, source):
101          """          """
102          Adds a survey data provider to the domain builder.          Adds a survey data provider to the domain builder.
# Line 70  class DomainBuilder(object): Line 108  class DomainBuilder(object):
108          dimensionality of the data must be one less than the dimensionality          dimensionality of the data must be one less than the dimensionality
109          of the domain (specified in the constructor).          of the domain (specified in the constructor).
110    
111          :param source: The data source to be added          :param source: The data source to be added. Its reference system needs
112                           to match the reference system of the DomainBuilder.
113          :type source: `DataSource`          :type source: `DataSource`
114          """          """
115          if self.__domain is not None:          if self.__domain is not None:
116              raise RuntimeError("Invalid call to addSource(). Domain is already built.")              raise RuntimeError("Invalid call to addSource(). Domain is already built.")
117          if not isinstance(source, DataSource):          if not isinstance(source, DataSource):
118              raise TypeError("source is not a DataSource")              raise TypeError("source is not a DataSource")
119            if not source.getReferenceSystem() == self.getReferenceSystem():
120               raise ValueError("source reference system does not match.")
121    
122          DATA_DIM = len(source.getDataExtents()[0])          DATA_DIM = len(source.getDataExtents()[0])
123          if DATA_DIM != self.__dim-1:          if DATA_DIM != self.__dim-1:
# Line 86  class DomainBuilder(object): Line 127  class DomainBuilder(object):
127                  raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")                  raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")
128    
129          self.__sources.append(source)          self.__sources.append(source)
130            if source.getTags(): self.__tags=list(set(self.__tags + source.getTags()))
131    
132      def setFractionalPadding(self, pad_x=None, pad_y=None):      def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
133          """          """
134          Sets the amount of padding around the dataset as a fraction of the          Sets the amount of padding around the dataset as a fraction of the
135          dataset side lengths.          dataset side lengths.
# Line 100  class DomainBuilder(object): Line 142  class DomainBuilder(object):
142          :type pad_x: ``float``          :type pad_x: ``float``
143          :param pad_y: Padding per side in y direction (default: no padding)          :param pad_y: Padding per side in y direction (default: no padding)
144          :type pad_y: ``float``          :type pad_y: ``float``
145          :note: `pad_y` is ignored for 2-dimensional domains.          :param pad_lat: Padding per side in latitudinal direction (default: no padding)
146          """          :type pad_lat: ``float``
147            :param pad_lon: Padding per side in longitudinal direction (default: no padding)
148            :type pad_lon: ``float``        
149            :note: `pad_y` is ignored for 2-dimensional domains.
150            """
151            if not pad_lat == None:
152                if not pad_x == None:
153                   raise ValueError("Either pad_lat or pad_x can be set.")
154                else:
155                  pad_x = pad_lat
156            if not pad_lon == None:
157                if not pad_y == None:
158                  raise ValueError("Either pad_lon or pad_y can be set.")
159                else:
160                  pad_y = pad_lan        
161          if self.__domain is not None:          if self.__domain is not None:
162              raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")              raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
163          if pad_x is not None:          if pad_x is not None:
# Line 116  class DomainBuilder(object): Line 172  class DomainBuilder(object):
172                  raise ValueError("setFractionalPadding: Argument too large")                  raise ValueError("setFractionalPadding: Argument too large")
173          self._padding = [pad_x,pad_y], 'f'          self._padding = [pad_x,pad_y], 'f'
174    
175      def setPadding(self, pad_x=None, pad_y=None):      def setPadding(self, pad_x=None, pad_y=None,  pad_lat=None, pad_lon=None):
176          """          """
177          Sets the amount of padding around the dataset in absolute length units.          Sets the amount of padding around the dataset in absolute length units.
178    
# Line 125  class DomainBuilder(object): Line 181  class DomainBuilder(object):
181          non-negative.          non-negative.
182    
183          :param pad_x: Padding per side in x direction (default: no padding)          :param pad_x: Padding per side in x direction (default: no padding)
184          :type pad_x: ``float``          :type pad_x: ``float`` in units of length (meter)
185          :param pad_y: Padding per side in y direction (default: no padding)          :param pad_y: Padding per side in y direction (default: no padding)
186          :type pad_y: ``float``          :type pad_y: ``float`` in units of length (meter)
187          :note: `pad_y` is ignored for 2-dimensional domains.          :note: `pad_y` is ignored for 2-dimensional domains.
188            :note: this function can only be used if the reference system is Cartesian
189          """          """
190            if not self.getReferenceSystem().isCartesian():
191                raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
192          if self.__domain is not None:          if self.__domain is not None:
193              raise RuntimeError("Invalid call to setPadding(). Domain is already built.")              raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
194          if pad_x is not None:          if pad_x is not None:
# Line 139  class DomainBuilder(object): Line 198  class DomainBuilder(object):
198              if pad_y < 0:              if pad_y < 0:
199                  raise ValueError("setPadding: Arguments must be non-negative")                  raise ValueError("setPadding: Arguments must be non-negative")
200          self._padding = [pad_x,pad_y], 'l'          self._padding = [pad_x,pad_y], 'l'
201            
202        def setGeoPadding(self, pad_lat=None, pad_lon=None):
203            """
204            Sets the amount of padding around the dataset in longitude and latitude.
205    
206      def setElementPadding(self, pad_x=None, pad_y=None):          The final domain size will be the extent in the latitudinal (in
207            longitudinal) direction of the dataset plus twice the value of
208            `pad_lat` (`pad_lon`). The arguments must be non-negative.
209    
210            :param pad_lat: Padding per side in latitudinal direction (default: 0)
211            :type pad_lat: ``float`` in units of degree
212            :param pad_lon: Padding per side in longitudinal direction (default: 0)
213            :type pad_lon: ``float``  in units of degree  
214            :note: `pad_lon` is ignored for 2-dimensional domains.
215            :note: this function can only be used if the reference system is not Cartesian
216            """
217            if self.getReferenceSystem().isCartesian():
218                raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
219            if self.__domain is not None:
220                raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
221            if pad_lat is not None:
222                if pad_lat < 0:
223                    raise ValueError("setPadding: Arguments must be non-negative")
224            if pad_lon is not None:
225                if pad_lon < 0:
226                    raise ValueError("setPadding: Arguments must be non-negative")
227            self._padding = [pad_lat,pad_lon], 'd'
228            
229        def setElementPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
230          """          """
231          Sets the amount of padding around the dataset in number of elements          Sets the amount of padding around the dataset in number of elements
232          (cells).          (cells).
# Line 155  class DomainBuilder(object): Line 241  class DomainBuilder(object):
241          :type pad_y: ``int``          :type pad_y: ``int``
242          :note: `pad_y` is ignored for 2-dimensional datasets.          :note: `pad_y` is ignored for 2-dimensional datasets.
243          """          """
244            if not pad_lat == None:
245                if not pad_x == None:
246                  raise ValueError("Either pad_lat or pad_x can be set.")
247                else:
248                  pad_x = pad_lat
249            if not pad_lon == None:
250                if not pad_y == None:
251                  raise ValueError("Either pad_lon or pad_y can be set.")
252                else:
253                  pad_y = pad_lan
254                  
255          if self.__domain is not None:          if self.__domain is not None:
256              raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")              raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
257          if pad_x is not None:          if pad_x is not None:
# Line 181  class DomainBuilder(object): Line 278  class DomainBuilder(object):
278          """          """
279          return self.getSurveys(DataSource.MAGNETIC)          return self.getSurveys(DataSource.MAGNETIC)
280    
281            
282      def fixDensityBelow(self, depth=None):      def fixDensityBelow(self, depth=None):
283          """          """
284          Defines the depth below which the density anomaly is set to a given value.          Defines the depth below which the density anomaly is set to a given
285          if no value is given zero is assumed.          value. If no value is given zero is assumed.
286                    
287          :param depth: depth below which the density is fixed. If not set , no constraint          :param depth: depth below which the density is fixed. If not set, no
288                        at depth is applied.                        constraint at depth is applied.
289          :type depth: ``float``          :type depth: ``float``
290          """          """
291          self.__fix_density_below=depth          self.__fix_density_below=depth
292    
293      def fixSusceptibilityBelow(self, depth=None):      def fixSusceptibilityBelow(self, depth=None):
294          """          """
295          Defines the depth below which the susceptibility anomaly is set to a given value.          Defines the depth below which the susceptibility anomaly is set to a
296          if no value is given zero is assumed.          given value. If no value is given zero is assumed.
297                    
298          :param depth: depth below which the susceptibility is fixed. If not set , no constraint          :param depth: depth below which the susceptibility is fixed. If not
299                        at depth is applied.                        set, no constraint at depth is applied.
300          :type depth: ``float``          :type depth: ``float``
301          """          """
302          self.__fix_susceptibility_below=depth          self.__fix_susceptibility_below=depth
303    
304      def getSurveys(self, datatype):      def fixVelocityBelow(self, depth=None):
305            """
306            Defines the depth below which the velocity and Q index is set to a
307            given value. If no value is given zero is assumed.
308            
309            :param depth: depth below which the velocity is fixed. If not
310                          set, no constraint at depth is applied.
311            :type depth: ``float``
312            """
313            self.__fix_velocity_below=depth
314    
315    
316        def getSurveys(self, datatype, tags=None):
317          """          """
318          Returns a list of `Data` objects for all surveys of type `datatype`          Returns a list of `Data` objects for all surveys of type `datatype`
319          available to this domain builder.          available to this domain builder. If a list of `tags` is given
320            only data sources whose tag matching the tag list are returned
321    
322          :return: List of surveys which are tuples (anomaly,error).          :return: List of surveys which are tuples (anomaly,error).
323          :rtype: ``list``          :rtype: ``list``
# Line 214  class DomainBuilder(object): Line 325  class DomainBuilder(object):
325          surveys=[]          surveys=[]
326          for src in self.__sources:          for src in self.__sources:
327              if src.getDataType()==datatype:              if src.getDataType()==datatype:
328                  surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))                  if tags is None or ( src.getTags() is not None and all( [ t in tags for t in src.getTags() ] )  ) :
329                        surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
330          return surveys          return surveys
331    
332      def setBackgroundMagneticFluxDensity(self, B):      def setBackgroundMagneticFluxDensity(self, B):
333          """          """
334          Sets the background magnetic flux density B=(B_North, B_East, B_Vertical)          Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
335          """          """
336          self.__background_magnetic_field=B          self.__background_magnetic_field=B
337    
# Line 228  class DomainBuilder(object): Line 340  class DomainBuilder(object):
340          Returns the background magnetic flux density.          Returns the background magnetic flux density.
341          """          """
342          B = self.__background_magnetic_field          B = self.__background_magnetic_field
343          # this is for Cartesian (FIXME ?)          if B is None:
344          if self.__dim<3:              raise ValueError("No background magnetic flux density set!")
345    
346            if self.__dim < 3 :
347              return np.array([B[0], B[2]])              return np.array([B[0], B[2]])
348          else:          else:
349              return np.array(B)              return np.array(B)
# Line 242  class DomainBuilder(object): Line 356  class DomainBuilder(object):
356          z=self.getDomain().getX()[self.__dim-1]          z=self.getDomain().getX()[self.__dim-1]
357          m = whereNonNegative(z)          m = whereNonNegative(z)
358          if self.__fix_density_below:          if self.__fix_density_below:
359              m += whereNonPositive(z+self.__fix_density_below)              m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
360          return m          return m
361    
362      def getSetSusceptibilityMask(self):      def getSetSusceptibilityMask(self):
# Line 253  class DomainBuilder(object): Line 367  class DomainBuilder(object):
367          z=self.getDomain().getX()[self.__dim-1]          z=self.getDomain().getX()[self.__dim-1]
368          m = whereNonNegative(z)          m = whereNonNegative(z)
369          if self.__fix_susceptibility_below:          if self.__fix_susceptibility_below:
370              m += whereNonPositive(z+self.__fix_susceptibility_below)              m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
371          return m          return m
372    
373      def getDomain(self):      def getDomain(self):
# Line 305  class DomainBuilder(object): Line 419  class DomainBuilder(object):
419              if pad[i] is None:              if pad[i] is None:
420                  frac.append(0.)                  frac.append(0.)
421                  continue                  continue
422              if pt == 'f': # fraction of side length              if pt == 'f' : # fraction of side length
423                  frac.append(2.*pad[i])                  frac.append(2.*pad[i])
424              elif pt == 'e': # number of elements              elif pt == 'e': # number of elements
425                  frac.append(2.*pad[i]/float(NX[i]))                  frac.append(2.*pad[i]/float(NX[i]))
# Line 327  class DomainBuilder(object): Line 441  class DomainBuilder(object):
441          if len(self.__sources)==0:          if len(self.__sources)==0:
442              raise ValueError("No data")              raise ValueError("No data")
443          X0, NE, DX = self.__sources[0].getDataExtents()          X0, NE, DX = self.__sources[0].getDataExtents()
444          XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]          # do not mess with the values if only one source used
445            if len(self.__sources)>1:
446          for src in self.__sources[1:]:              XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
447              d_x0, d_ne, d_dx = src.getDataExtents()  
448              for i in range(len(d_x0)):              for src in self.__sources[1:]:
449                  X0[i]=min(X0[i], d_x0[i])                  d_x0, d_ne, d_dx = src.getDataExtents()
450              for i in range(len(d_dx)):                  for i in range(len(d_x0)):
451                  DX[i]=min(DX[i], d_dx[i])                      X0[i]=min(X0[i], d_x0[i])
452              for i in range(len(d_ne)):                  for i in range(len(d_dx)):
453                  XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])                      DX[i]=min(DX[i], d_dx[i])
454          NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]                  for i in range(len(d_ne)):
455                        XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
456                # FIXME: should this be rounded up instead?
457                NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
458          return X0, NE, DX          return X0, NE, DX
459    
460      def __createDomain(self):      def __createDomain(self):
# Line 355  class DomainBuilder(object): Line 472  class DomainBuilder(object):
472          NE = NX + [self._v_num_cells]          NE = NX + [self._v_num_cells]
473    
474          # origin of domain          # origin of domain
475          origin = X0 + [-self._v_depth]          origin = X0 + [-self._v_depth*self.__v_scale]
476          self._dom_origin = [np.floor(oi) for oi in origin]  
477            if self.getReferenceSystem().isCartesian():
478                # rounding will give us about meter-accuracy with UTM coordinates
479                self._dom_origin = [np.floor(oi) for oi in origin]
480            else:
481                # this should give us about meter-accuracy with lat/lon coords
482                self._dom_origin = [1e-5*np.floor(oi*1e5) for oi in origin]
483    
484          # cell size / point spacing          # cell size / point spacing
485          spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]          spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
486          self._spacing = [float(np.floor(si)) for si in spacing]          #self._spacing = [float(np.floor(si)) for si in spacing]
487            self._spacing = spacing
488    
489          lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]          lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
490    
491          if self.__dim==3:          if self.__dim==3:
492              dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])              dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
493          else:          else:

Legend:
Removed from v.4362  
changed lines
  Added in v.4729

  ViewVC Help
Powered by ViewVC 1.1.26