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

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

  ViewVC Help
Powered by ViewVC 1.1.26