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

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

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

revision 4115 by gross, Fri Dec 14 04:48:48 2012 UTC revision 4145 by caltinay, Fri Jan 18 00:51:49 2013 UTC
# Line 1  Line 1 
1    
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2013 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 15  Line 15 
15    
16  """Data readers/providers for inversions"""  """Data readers/providers for inversions"""
17    
18  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19  http://www.uq.edu.au  http://www.uq.edu.au
20  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
21  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
22  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
23  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
24    
25  __all__ = ['simpleGeoMagneticFluxDensity', 'DataSource','ErMapperData', 'SyntheticDataBase' , 'SyntheticFeatureData', 'SyntheticData','SmoothAnomaly']  __all__ = ['simpleGeoMagneticFluxDensity', 'DataSource', 'ErMapperData', \
26            'SyntheticDataBase', 'SyntheticFeatureData', 'SyntheticData',
27            'SmoothAnomaly']
28    
29  import logging  import logging
30  import numpy as np  import numpy as np
31  from esys.escript import ReducedFunction  from esys.escript import ReducedFunction
32    from esys.escript import unitsSI as U
33  from esys.escript.linearPDEs import LinearSinglePDE  from esys.escript.linearPDEs import LinearSinglePDE
34  from esys.escript.util import *  from esys.escript.util import *
 import esys.escript.unitsSI as U  
35  from esys.ripley import Brick, Rectangle, ripleycpp  from esys.ripley import Brick, Rectangle, ripleycpp
36    
37  try:  try:
# Line 79  def LatLonToUTM(lon, lat, wkt_string=Non Line 81  def LatLonToUTM(lon, lat, wkt_string=Non
81      return x,y      return x,y
82    
83  def simpleGeoMagneticFluxDensity(latitude, longitude=0.):  def simpleGeoMagneticFluxDensity(latitude, longitude=0.):
84          theta = (90-latitude)/180.*np.pi      theta = (90-latitude)/180.*np.pi
85          B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)      B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)
86          B_theta= B_0 * sin(theta)      B_theta= B_0 * sin(theta)
87          B_r= 2 * B_0 * cos(theta)      B_r= 2 * B_0 * cos(theta)
88          return B_r, B_theta, 0.      return B_r, B_theta, 0.
89    
90  class DataSource(object):  class DataSource(object):
91      """      """
# Line 102  class DataSource(object): Line 104  class DataSource(object):
104          """          """
105          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
106          self.__subsampling_factor=1          self.__subsampling_factor=1
         self.__background_magnetic_field=None  
107    
108      def getDataExtents(self):      def getDataExtents(self):
109          """          """
# Line 127  class DataSource(object): Line 128  class DataSource(object):
128          """          """
129          This method is called by the `DomainBuilder` to retrieve the survey          This method is called by the `DomainBuilder` to retrieve the survey
130          data as `Data` objects on the given domain.          data as `Data` objects on the given domain.
131    
132          Subclasses should return one or more `Data` objects with survey data          Subclasses should return one or more `Data` objects with survey data
133          interpolated on the given ripley domain. The exact return type          interpolated on the given ripley domain. The exact return type
134          depends on the type of data.          depends on the type of data.
# Line 145  class DataSource(object): Line 147  class DataSource(object):
147      def setSubsamplingFactor(self, f):      def setSubsamplingFactor(self, f):
148          """          """
149          Sets the data subsampling factor (default=1).          Sets the data subsampling factor (default=1).
150    
151          The factor is applied in all dimensions. For example a 2D dataset          The factor is applied in all dimensions. For example a 2D dataset
152          with 300 x 150 data points will be reduced to 150 x 75 when a          with 300 x 150 data points will be reduced to 150 x 75 when a
153          subsampling factor of 2 is used.          subsampling factor of 2 is used.
# Line 516  class SmoothAnomaly(SourceFeature): Line 519  class SmoothAnomaly(SourceFeature):
519    
520  ##############################################################################  ##############################################################################
521  class SyntheticDataBase(DataSource):  class SyntheticDataBase(DataSource):
522    """      """
523    Base class to define reference data based on a given property distribution (density or      Base class to define reference data based on a given property distribution
524    susceptibility). Data are collected from a square region of vertical extend `length` on a      (density or susceptibility). Data are collected from a square region of
525    grid with ``number_of_elements`` cells in each direction.      vertical extent `length` on a grid with `number_of_elements` cells in each
526          direction.
527    The synthetic data are constructed by solving the appropriate forward problem. Data can be sampled      
528    with an offset from the surface at z=0 or using the entire subsurface region.      The synthetic data are constructed by solving the appropriate forward
529    """      problem. Data can be sampled with an offset from the surface at z=0 or
530    def __init__(self, datatype,      using the entire subsurface region.
531          DIM=2,      """
532          number_of_elements=10,      def __init__(self, datatype,
533          length=1*U.km,                   DIM=2,
534          B_b=None,                   number_of_elements=10,
535          data_offset=0,                   length=1*U.km,
536          full_knowledge=False,                   B_b=None,
537          spherical=False):                   data_offset=0,
538      """                   full_knowledge=False,
539                     spherical=False):
540            """
541          :param datatype: data type indicator          :param datatype: data type indicator
542          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
543          :param DIM: spatial dimension          :param DIM: number of spatial dimension
544          :type DIM: 2 or 3          :type DIM: ``int`` (2 or 3)
545          :param number_of_elements: lateral number of elements in the region where data are collected          :param number_of_elements: lateral number of elements in the region
546                                       where data are collected
547          :type number_of_elements: ``int``          :type number_of_elements: ``int``
548          :param length: lateral extend of the region where data are collected          :param length: lateral extend of the region where data are collected
549          :type length: ``float``          :type length: ``float``
# Line 545  class SyntheticDataBase(DataSource): Line 551  class SyntheticDataBase(DataSource):
551          :type B_b: ``list`` of ``Scalar``          :type B_b: ``list`` of ``Scalar``
552          :param data_offset: offset of the data collection region from the surface          :param data_offset: offset of the data collection region from the surface
553          :type data_offset: ``float``          :type data_offset: ``float``
554          :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.          :param full_knowledge: if ``True`` data are collected from the entire
555                                   subsurface region. This is mainly for testing.
556          :type full_knowledge: ``Bool``          :type full_knowledge: ``Bool``
557          :param spherical: if ``True`` sperical coordinates are used (ignored)          :param spherical: if ``True`` sperical coordinates are used (ignored for now)
558          :type spherical: ``Bool``          :type spherical: ``Bool``
559      """          """
560      super(SyntheticDataBase,self).__init__()          super(SyntheticDataBase,self).__init__()
561      if not datatype in [self.GRAVITY,self.MAGNETIC]:          if not datatype in [self.GRAVITY,self.MAGNETIC]:
562          raise ValueError("Invalid value for datatype parameter")                    raise ValueError("Invalid value for datatype parameter")      
563      self.DIM=DIM          self.DIM=DIM
564      self.number_of_elements=number_of_elements          self.number_of_elements=number_of_elements
565      self.length=length          self.length=length
566      self.__datatype = datatype          self.__datatype = datatype
567                
568      self.__spherical = spherical            self.__spherical = spherical  
569      self.__full_knowledge= full_knowledge          self.__full_knowledge= full_knowledge
570      self.__data_offset=data_offset          self.__data_offset=data_offset
571      self.__B_b =None          self.__B_b =None
572      # this is for Cartesian (FIXME ?)          # this is for Cartesian (FIXME ?)
573      if datatype  ==  self.MAGNETIC:          if datatype  ==  self.MAGNETIC:
574          if self.DIM<3:              if self.DIM < 3:
575            self.__B_b =  np.array([-B_b[2],  -B_b[0]])                  self.__B_b = np.array([-B_b[2], -B_b[0]])
576          else:              else:
577            self.__B_b = ([-B_b[1],  -B_b[2],  -B_b[0]])                      self.__B_b = ([-B_b[1], -B_b[2], -B_b[0]])    
578      self.__origin = [0]*(DIM-1)          self.__origin = [0]*(DIM-1)
579      self.__delta = [float(length)/number_of_elements]*(DIM-1)          self.__delta = [float(length)/number_of_elements]*(DIM-1)
580      self.__nPts = [number_of_elements]*(DIM-1)          self.__nPts = [number_of_elements]*(DIM-1)
581      self._reference_data=None          self._reference_data=None
582                
583    def getDataExtents(self):      def getDataExtents(self):
584      """          """
585      returns the lateral data extend of the data set          returns the lateral data extend of the data set
586      """          """
587      return (list(self.__origin), list(self.__nPts), list(self.__delta))          return (list(self.__origin), list(self.__nPts), list(self.__delta))
588    
589    def getDataType(self):      def getDataType(self):
590      """          """
591      returns the data type          returns the data type
592      """          """
593      return self.__datatype          return self.__datatype
594                
595    def getSurveyData(self, domain, origin, number_of_elements, spacing):      def getSurveyData(self, domain, origin, number_of_elements, spacing):
596      """          """
597      returns the survey data placed on a given domain.          returns the survey data placed on a given domain.
598                
599      :param domain: domain on which the data are to be placed          :param domain: domain on which the data are to be placed
600      :type param: ``Domain``          :type domain: ``Domain``
601      :param origin: origin of the domain          :param origin: origin of the domain
602      :type origin: ``list`` of ``float``          :type origin: ``list`` of ``float``
603      :param number_of_elements: number of elements (or cells) in each spatial direction used          :param number_of_elements: number of elements (or cells) in each
604                    span the domain                                     spatial direction used to span the domain
605      :type number_of_elements: `list`` of ``int``          :type number_of_elements: ``list`` of ``int``
606      :param spacing: cell size in each spatial direction          :param spacing: cell size in each spatial direction
607      :type spacing: ``list`` of ``float``          :type spacing: ``list`` of ``float``
608      :return: observed gravity field or magnetic flux density for each cell in the domain and          :return: observed gravity field or magnetic flux density for each cell
609      for each cell an indicator 1/0 if the data are valid or not.                   in the domain and for each cell an indicator 1/0 if the data
610      :rtype: pair of ``Scalar``                   are valid or not.
611      """          :rtype: pair of ``Scalar``
612      pde=LinearSinglePDE(domain)          """
613      DIM=domain.getDim()          pde=LinearSinglePDE(domain)
614      x=domain.getX()          DIM=domain.getDim()
615      # set the reference data          x=domain.getX()
616                # set the reference data
617      k=self.getReferenceProperty(domain)          
618      # calculate the corresponding potential          k=self.getReferenceProperty(domain)
619      z=x[DIM-1]          # calculate the corresponding potential
620      m_psi_ref=whereZero(z-sup(z))          z=x[DIM-1]
621      if self.getDataType()==DataSource.GRAVITY:          m_psi_ref=whereZero(z-sup(z))
622          pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)          if self.getDataType()==DataSource.GRAVITY:
623      else:              pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
624          pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)          else:
625      pde.setSymmetryOn()              pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
626      psi_ref=pde.getSolution()          pde.setSymmetryOn()
627      del pde          #pde.getSolverOptions().setTolerance(1e-13)
628      if self.getDataType()==DataSource.GRAVITY:          psi_ref=pde.getSolution()
629          data = -grad(psi_ref, ReducedFunction(domain))          del pde
630      else:          if self.getDataType()==DataSource.GRAVITY:
631          data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))              data = -grad(psi_ref, ReducedFunction(domain))
632                  else:
633      x=ReducedFunction(domain).getX()                  data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
634      if self.__full_knowledge:  
635          sigma = whereNegative(x[DIM-1])          x=ReducedFunction(domain).getX()    
636      else:          if self.__full_knowledge:
637                      sigma = whereNegative(x[DIM-1])
638          sigma=1.          else:
639          # limit mask to non-padding in horizontal area                      sigma=1.
640          for i in range(DIM-1):              # limit mask to non-padding in horizontal area        
641          x_i=x[i]              for i in range(DIM-1):
642          sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))                  x_i=x[i]
643          # limit mask to one cell thickness at z=0                  sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
644          z=x[DIM-1]              # limit mask to one cell thickness at z=0
645          oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]              z=x[DIM-1]
646          sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])              oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
647      return data,sigma              sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
648                return data,sigma
649    def getReferenceProperty(self, domain=None):          
650      """      def getReferenceProperty(self, domain=None):
651      Returns the reference density Data object that was used to generate          """
652      the gravity/susceptibility anomaly data.          Returns the reference density Data object that was used to generate
653                the gravity/susceptibility anomaly data.
654      :return: the density or susceptibility anomaly used to create the survey data.          
655      :note: it can be assumed that in the first call the ``domain`` argument is present so the          :return: the density or susceptibility anomaly used to create the
656      actual anomaly data can be created. In subsequent calls this may not be true.                   survey data
657      :note: method needs to be overwritten          :note: it can be assumed that in the first call the ``domain``
658      """                 argument is present so the actual anomaly data can be created.
659      raise NotImplementedError()                       In subsequent calls this may not be true.
660            :note: method needs to be overwritten
661            """
662            raise NotImplementedError()      
663                
664  class SyntheticFeatureData(SyntheticDataBase):  class SyntheticFeatureData(SyntheticDataBase):
665      """      """
# Line 695  class SyntheticFeatureData(SyntheticData Line 705  class SyntheticFeatureData(SyntheticData
705    
706      def getReferenceProperty(self, domain=None):      def getReferenceProperty(self, domain=None):
707          """          """
708      Returns the reference density Data object that was used to generate          Returns the reference density Data object that was used to generate
709      the gravity/susceptibility anomaly data.          the gravity/susceptibility anomaly data.
710          """          """
711          if self._reference_data == None:          if self._reference_data == None:
712              DIM=domain.getDim()              DIM=domain.getDim()
# Line 712  class SyntheticData(SyntheticDataBase): Line 722  class SyntheticData(SyntheticDataBase):
722      """      """
723      defines synthetic  gravity/magnetic data based on harmonic property anomaly      defines synthetic  gravity/magnetic data based on harmonic property anomaly
724            
725          rho = mean + amplitude * sin(n_depth * pi /depth * z) * sin(n_length * pi /length * x - shift )          rho = mean + amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * x - shift )
726                    
727      for all x and z<=0. for z>0 rho = 0.              for all x and z<=0. for z>0 rho = 0.        
728      """      """
729      def __init__(self, datatype,      def __init__(self, datatype,
730                         n_length=1,                         n_length=1,
731                         n_depth=1,                         n_depth=1,
732                         shift=0.,                         depth_offset=0.,
733                           depth=None,
734                         amplitude=None,                         amplitude=None,
735                         DIM=2,                         DIM=2,
736                         number_of_elements=10,                         number_of_elements=10,
# Line 731  class SyntheticData(SyntheticDataBase): Line 742  class SyntheticData(SyntheticDataBase):
742          """          """
743          :param datatype: data type indicator          :param datatype: data type indicator
744          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
745          :param n_length: number of oscillation in the anomaly data within the observation region.          :param n_length: number of oscillations in the anomaly data within the
746                             observation region
747          :type n_length: ``int``          :type n_length: ``int``
748          :param n_depth: number of oscillation in the anomaly data below surface          :param n_depth: number of oscillations in the anomaly data below surface
749          :param shift: lateral phase shift in the  anomaly data          :param depth: vertical extent in the anomaly data. If not present the
750          :param amplitude: data amplitude. Default value is 200 *U.kg/U.m**3 for gravity and 0.1 for magnetic data.                        depth of the domain is used.
751          :param features: list of features. It is recommended that the features do entirely lay below surface.          :type depth: ``float``
752          :type features: ``list`` of ``SourceFeature``          :param amplitude: data amplitude. Default value is 200 U.kg/U.m**3 for
753          :param DIM: spatial dimension                            gravity and 0.1 for magnetic data.
754          :type DIM: 2 or 3          :param DIM: spatial dimensionality
755          :param number_of_elements: lateral number of elements in the region where data are collected          :type DIM: ``int`` (2 or 3)
756            :param number_of_elements: lateral number of elements in the region
757                                       where data are collected
758          :type number_of_elements: ``int``          :type number_of_elements: ``int``
759          :param length: lateral extend of the region where data are collected          :param length: lateral extent of the region where data are collected
760          :type length: ``float``          :type length: ``float``
761          :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.          :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude].
762                        Only used for magnetic data.
763          :type B_b: ``list`` of ``Scalar``          :type B_b: ``list`` of ``Scalar``
764          :param data_offset: offset of the data collection region from the surface          :param data_offset: offset of the data collection region from the surface
765          :type data_offset: ``float``          :type data_offset: ``float``
766          :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.          :param full_knowledge: if ``True`` data are collected from the entire
767                                   subsurface region. This is mainly for testing.
768          :type full_knowledge: ``Bool``          :type full_knowledge: ``Bool``
769          :param spherical: if ``True`` sperical coordinates are used (ignored)          :param spherical: if ``True`` sperical coordinates are used (ignored for now)
770          :type spherical: ``Bool``          :type spherical: ``Bool``
771          """                """      
772          super(SyntheticData,self).__init__(          super(SyntheticData,self).__init__(
773                                   datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,                                   datatype=datatype, DIM=DIM,
774                                     number_of_elements=number_of_elements,
775                                   length=length, B_b=B_b,                                   length=length, B_b=B_b,
776                                   data_offset=data_offset,                                   data_offset=data_offset,
777                                   full_knowledge=full_knowledge,                                   full_knowledge=full_knowledge,
778                                   spherical=spherical)                                   spherical=spherical)
779          self.__n_length = n_length          self.__n_length = n_length
780          self.__n_depth = n_depth          self.__n_depth = n_depth
781          self.__shift = shift          self.depth = depth
782            self.depth_offset=depth_offset
783          if amplitude == None:          if amplitude == None:
784          if datatype == DataSource.GRAVITY:              if datatype == DataSource.GRAVITY:
785              amplitude = 200 *U.kg/U.m**3                  amplitude = 200 *U.kg/U.m**3
786          else:              else:
787               amplitude =0.1                  amplitude = 0.1
788          self.__amplitude = amplitude          self.__amplitude = amplitude
789    
   
   
790      def getReferenceProperty(self, domain=None):      def getReferenceProperty(self, domain=None):
791          """          """
792          Returns the reference density Data object that was used to generate          Returns the reference density Data object that was used to generate
# Line 781  class SyntheticData(SyntheticDataBase): Line 797  class SyntheticData(SyntheticDataBase):
797              x=domain.getX()              x=domain.getX()
798              # set the reference data              # set the reference data
799              z=x[DIM-1]              z=x[DIM-1]
800              k=sin(self.__n_depth * np.pi/inf(z) * z) * whereNegative(z) * self.__amplitude              dd=self.depth
801                if  dd ==None :  dd=inf(z)
802                z2=(z+self.depth_offset)/(self.depth_offset-dd)
803                k=sin(self.__n_depth * np.pi  * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude
804              for i in xrange(DIM-1):              for i in xrange(DIM-1):
805             k*= sin(self.__n_length * np.pi /self.__length * x[i] - shift )                 x_i=x[i]
806                   min_x=inf(x_i)
807                   max_x=sup(x_i)
808                   k*= sin(self.__n_length*np.pi*(x_i-min_x)/(max_x-min_x))
809              self._reference_data= k              self._reference_data= k
810          return self._reference_data          return self._reference_data
811    
   

Legend:
Removed from v.4115  
changed lines
  Added in v.4145

  ViewVC Help
Powered by ViewVC 1.1.26