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

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

  ViewVC Help
Powered by ViewVC 1.1.26