/[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 4108 by caltinay, Thu Dec 13 06:38:11 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__ = ['simpleBackgroundMagneticField', 'DataSource','ErMapperData','SyntheticFeatureData','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 78  def LatLonToUTM(lon, lat, wkt_string=Non Line 80  def LatLonToUTM(lon, lat, wkt_string=Non
80      x,y=pyproj.transform(p_src, p_dest, lon, lat)      x,y=pyproj.transform(p_src, p_dest, lon, lat)
81      return x,y      return x,y
82    
83  def simpleBackgroundMagneticField(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 515  class SmoothAnomaly(SourceFeature): Line 518  class SmoothAnomaly(SourceFeature):
518          return m          return m
519    
520  ##############################################################################  ##############################################################################
521  class SyntheticFeatureData(DataSource):  class SyntheticDataBase(DataSource):
522      def __init__(self, datatype, DIM, NE, l, features, B_b=None):      """
523          super(SyntheticFeatureData,self).__init__()      Base class to define reference data based on a given property distribution
524        (density or susceptibility). Data are collected from a square region of
525        vertical extent `length` on a grid with `number_of_elements` cells in each
526        direction.
527        
528        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        using the entire subsurface region.
531        """
532        def __init__(self, datatype,
533                     DIM=2,
534                     number_of_elements=10,
535                     length=1*U.km,
536                     B_b=None,
537                     data_offset=0,
538                     full_knowledge=False,
539                     spherical=False):
540            """
541            :param datatype: data type indicator
542            :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
543            :param DIM: number of spatial dimension
544            :type DIM: ``int`` (2 or 3)
545            :param number_of_elements: lateral number of elements in the region
546                                       where data are collected
547            :type number_of_elements: ``int``
548            :param length: lateral extend of the region where data are collected
549            :type length: ``float``
550            :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
551            :type B_b: ``list`` of ``Scalar``
552            :param data_offset: offset of the data collection region from the surface
553            :type data_offset: ``float``
554            :param full_knowledge: if ``True`` data are collected from the entire
555                                   subsurface region. This is mainly for testing.
556            :type full_knowledge: ``Bool``
557            :param spherical: if ``True`` sperical coordinates are used (ignored for now)
558            :type spherical: ``Bool``
559            """
560            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
564            self.number_of_elements=number_of_elements
565            self.length=length
566          self.__datatype = datatype          self.__datatype = datatype
567          self._features = features          
568          self.__origin = [0]*(DIM-1)          self.__spherical = spherical  
569          self.__nPts = [NE]*(DIM-1)          self.__full_knowledge= full_knowledge
570          self.__delta = [float(l)/NE]*(DIM-1)          self.__data_offset=data_offset
571          self.__B_b =None          self.__B_b =None
         self.DIM=DIM  
         self.NE=NE  
         self.l=l  
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)
579            self.__delta = [float(length)/number_of_elements]*(DIM-1)
580            self.__nPts = [number_of_elements]*(DIM-1)
581            self._reference_data=None
582            
583      def getDataExtents(self):      def getDataExtents(self):
584            """
585            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):
         return self.__datatype  
   
     def getReferenceDensity(self):  
590          """          """
591          Returns the reference density Data object that was used to generate          returns the data type
         the gravity anomaly data.  
592          """          """
593          return self._rho          return self.__datatype
594            
595      def getReferenceSusceptibility(self):      def getSurveyData(self, domain, origin, number_of_elements, spacing):
596          """          """
597          Returns the reference magnetic susceptibility Data objects that was          returns the survey data placed on a given domain.
598          used to generate the magnetic field data.          
599            :param domain: domain on which the data are to be placed
600            :type domain: ``Domain``
601            :param origin: origin of the domain
602            :type origin: ``list`` of ``float``
603            :param number_of_elements: number of elements (or cells) in each
604                                       spatial direction used to span the domain
605            :type number_of_elements: ``list`` of ``int``
606            :param spacing: cell size in each spatial direction
607            :type spacing: ``list`` of ``float``
608            :return: observed gravity field or magnetic flux density for each cell
609                     in the domain and for each cell an indicator 1/0 if the data
610                     are valid or not.
611            :rtype: pair of ``Scalar``
612          """          """
         return self._k  
   
     def getSurveyData(self, domain, origin, NE, spacing):  
613          pde=LinearSinglePDE(domain)          pde=LinearSinglePDE(domain)
         G=U.Gravitational_Constant  
         m_psi_ref=0.  
         x=domain.getX()  
614          DIM=domain.getDim()          DIM=domain.getDim()
615          m_psi_ref=whereZero(x[DIM-1]-sup(x[DIM-1])) # + whereZero(x[DIM-1]-inf(x[DIM-1]))          x=domain.getX()
616            # set the reference data
617            
618            k=self.getReferenceProperty(domain)
619            # calculate the corresponding potential
620            z=x[DIM-1]
621            m_psi_ref=whereZero(z-sup(z))
622          if self.getDataType()==DataSource.GRAVITY:          if self.getDataType()==DataSource.GRAVITY:
623              rho_ref=0.              pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
             for f in self._features:  
                 m=f.getMask(x)  
                 rho_ref = rho_ref * (1-m) + f.getValue(x) * m  
             self._rho=rho_ref  
             pde.setValue(A=kronecker(domain), Y=-4*np.pi*G*rho_ref, q=m_psi_ref)  
624          else:          else:
625              k_ref=0.              pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
             for f in self._features:  
                 m=f.getMask(x)  
                 k_ref = k_ref * (1-m) + f.getValue(x) * m  
             self._k=k_ref  
             pde.setValue(A=kronecker(domain), X=k_ref*self.__B_b, q=m_psi_ref)  
626          pde.setSymmetryOn()          pde.setSymmetryOn()
627            #pde.getSolverOptions().setTolerance(1e-13)
628          psi_ref=pde.getSolution()          psi_ref=pde.getSolution()
629          del pde          del pde
630          if self.getDataType()==DataSource.GRAVITY:          if self.getDataType()==DataSource.GRAVITY:
631              data = -grad(psi_ref, ReducedFunction(domain))              data = -grad(psi_ref, ReducedFunction(domain))
632          else:          else:
633              data = self._k*self.__B_b-grad(psi_ref, ReducedFunction(domain))              data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
634    
635          sigma=1.          x=ReducedFunction(domain).getX()    
636          # limit mask to non-padding in horizontal area          if self.__full_knowledge:
637          x=ReducedFunction(domain).getX()              sigma = whereNegative(x[DIM-1])
638          for i in range(self.DIM-1):          else:
639              sigma=sigma * wherePositive(x[i]) \              sigma=1.
640                          * whereNegative(x[i]-(sup(x[i])+inf(x[i])))              # limit mask to non-padding in horizontal area        
641          # limit mask to one cell thickness at z=0              for i in range(DIM-1):
642          sigma = sigma * whereNonNegative(x[self.DIM-1]) \                  x_i=x[i]
643                  * whereNonPositive(x[self.DIM-1]-spacing[self.DIM-1])                  sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
644                # limit mask to one cell thickness at z=0
645                z=x[DIM-1]
646                oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
647                sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
648          return data,sigma          return data,sigma
649            
650        def getReferenceProperty(self, domain=None):
651            """
652            Returns the reference density Data object that was used to generate
653            the gravity/susceptibility anomaly data.
654            
655            :return: the density or susceptibility anomaly used to create the
656                     survey data
657            :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                   In subsequent calls this may not be true.
660            :note: method needs to be overwritten
661            """
662            raise NotImplementedError()      
663          
664    class SyntheticFeatureData(SyntheticDataBase):
665        """
666        uses a list of ``SourceFeature`` to define synthetic anomaly data
667        """
668        def __init__(self, datatype,
669                           features,
670                           DIM=2,
671                           number_of_elements=10,
672                           length=1*U.km,
673                           B_b=None,
674                           data_offset=0,
675                           full_knowledge=False,
676                           spherical=False):
677            """
678            :param datatype: data type indicator
679            :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
680            :param features: list of features. It is recommended that the features do entirely lay below surface.
681            :type features: ``list`` of ``SourceFeature``
682            :param DIM: spatial dimension
683            :type DIM: 2 or 3
684            :param number_of_elements: lateral number of elements in the region where data are collected
685            :type number_of_elements: ``int``
686            :param length: lateral extend of the region where data are collected
687            :type length: ``float``
688            :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
689            :type B_b: ``list`` of ``Scalar``
690            :param data_offset: offset of the data collection region from the surface
691            :type data_offset: ``float``
692            :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
693            :type full_knowledge: ``Bool``
694            :param spherical: if ``True`` sperical coordinates are used (ignored)
695            :type spherical: ``Bool``
696            """
697            super(SyntheticFeatureData,self).__init__(
698                                     datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,
699                                     length=length, B_b=B_b,
700                                     data_offset=data_offset,
701                                     full_knowledge=full_knowledge,
702                                     spherical=spherical)
703            self._features = features
704    
705    
706        def getReferenceProperty(self, domain=None):
707            """
708            Returns the reference density Data object that was used to generate
709            the gravity/susceptibility anomaly data.
710            """
711            if self._reference_data == None:
712                DIM=domain.getDim()
713                x=domain.getX()
714                k=0.
715                for f in self._features:
716                    m=f.getMask(x)
717                    k = k * (1-m) + f.getValue(x) * m
718                self._reference_data= k
719            return self._reference_data
720            
721    class SyntheticData(SyntheticDataBase):
722        """
723        defines synthetic  gravity/magnetic data based on harmonic property anomaly
724        
725            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.        
728        """
729        def __init__(self, datatype,
730                           n_length=1,
731                           n_depth=1,
732                           depth_offset=0.,
733                           depth=None,
734                           amplitude=None,
735                           DIM=2,
736                           number_of_elements=10,
737                           length=1*U.km,
738                           B_b=None,
739                           data_offset=0,
740                           full_knowledge=False,
741                           spherical=False):
742            """
743            :param datatype: data type indicator
744            :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
745            :param n_length: number of oscillations in the anomaly data within the
746                             observation region
747            :type n_length: ``int``
748            :param n_depth: number of oscillations in the anomaly data below surface
749            :param depth: vertical extent in the anomaly data. If not present the
750                          depth of the domain is used.
751            :type depth: ``float``
752            :param amplitude: data amplitude. Default value is 200 U.kg/U.m**3 for
753                              gravity and 0.1 for magnetic data.
754            :param DIM: spatial dimensionality
755            :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``
759            :param length: lateral extent of the region where data are collected
760            :type length: ``float``
761            :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``
764            :param data_offset: offset of the data collection region from the surface
765            :type data_offset: ``float``
766            :param full_knowledge: if ``True`` data are collected from the entire
767                                   subsurface region. This is mainly for testing.
768            :type full_knowledge: ``Bool``
769            :param spherical: if ``True`` sperical coordinates are used (ignored for now)
770            :type spherical: ``Bool``
771            """      
772            super(SyntheticData,self).__init__(
773                                     datatype=datatype, DIM=DIM,
774                                     number_of_elements=number_of_elements,
775                                     length=length, B_b=B_b,
776                                     data_offset=data_offset,
777                                     full_knowledge=full_knowledge,
778                                     spherical=spherical)
779            self.__n_length = n_length
780            self.__n_depth = n_depth
781            self.depth = depth
782            self.depth_offset=depth_offset
783            if amplitude == None:
784                if datatype == DataSource.GRAVITY:
785                    amplitude = 200 *U.kg/U.m**3
786                else:
787                    amplitude = 0.1
788            self.__amplitude = amplitude
789    
790        def getReferenceProperty(self, domain=None):
791            """
792            Returns the reference density Data object that was used to generate
793            the gravity anomaly data.
794            """
795            if self._reference_data == None:
796                DIM=domain.getDim()
797                x=domain.getX()
798                # set the reference data
799                z=x[DIM-1]
800                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):
805                   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
810            return self._reference_data
811    

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

  ViewVC Help
Powered by ViewVC 1.1.26