/[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 4061 by caltinay, Fri Nov 9 04:05:28 2012 UTC revision 4115 by gross, Fri Dec 14 04:48:48 2012 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__ = ['DataSource','ErMapperData','SyntheticData','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 simpleGeoMagneticFluxDensity(latitude, longitude=0.):
82            theta = (90-latitude)/180.*np.pi
83            B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)
84            B_theta= B_0 * sin(theta)
85            B_r= 2 * B_0 * cos(theta)
86            return B_r, B_theta, 0.
87    
88  class DataSource(object):  class DataSource(object):
89      """      """
# Line 85  class DataSource(object): Line 91  class DataSource(object):
91      This is an abstract base class that implements common functionality.      This is an abstract base class that implements common functionality.
92      Methods to be overwritten by subclasses are marked as such.      Methods to be overwritten by subclasses are marked as such.
93      This class assumes 2D data which is mapped to a slice of a 3D domain.      This class assumes 2D data which is mapped to a slice of a 3D domain.
94      For other setups override the `_createDomain` method.      For other setups override the methods as required.
95      """      """
96    
97      GRAVITY, MAGNETIC = list(range(2))      GRAVITY, MAGNETIC = list(range(2))
# Line 96  class DataSource(object): Line 102  class DataSource(object):
102          """          """
103          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
104          self.__subsampling_factor=1          self.__subsampling_factor=1
105            self.__background_magnetic_field=None
106    
107      def getDataExtents(self):      def getDataExtents(self):
108          """          """
# Line 317  class NetCdfData(DataSource): Line 324  class NetCdfData(DataSource):
324          :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`          :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`
325          :type datatype: ``int``          :type datatype: ``int``
326          :param altitude: altitude of measurements in meters          :param altitude: altitude of measurements in meters
327          :type datatype: ``float``          :type altitude: ``float``
328          """          """
329          super(NetCdfData,self).__init__()          super(NetCdfData,self).__init__()
330          self.__filename=filename          self.__filename=filename
# Line 410  class NetCdfData(DataSource): Line 417  class NetCdfData(DataSource):
417          #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]          #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
418          lon_range=longitude.data.min(),longitude.data.max()          lon_range=longitude.data.min(),longitude.data.max()
419          lat_range=latitude.data.min(),latitude.data.max()          lat_range=latitude.data.min(),latitude.data.max()
420          lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)          if lon_range[1]<180:
421                lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
422          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
423          f.close()          f.close()
424    
# Line 507  class SmoothAnomaly(SourceFeature): Line 515  class SmoothAnomaly(SourceFeature):
515          return m          return m
516    
517  ##############################################################################  ##############################################################################
518  class SyntheticData(DataSource):  class SyntheticDataBase(DataSource):
519      def __init__(self, datatype, DIM, NE, l, features):    """
520          super(SyntheticData,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        psi_ref=pde.getSolution()
620        del pde
621        if self.getDataType()==DataSource.GRAVITY:
622            data = -grad(psi_ref, ReducedFunction(domain))
623        else:
624            data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
625          
626        x=ReducedFunction(domain).getX()    
627        if self.__full_knowledge:
628            sigma = whereNegative(x[DIM-1])
629        else:
630          
631            sigma=1.
632            # limit mask to non-padding in horizontal area        
633            for i in range(DIM-1):
634            x_i=x[i]
635            sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
636            # limit mask to one cell thickness at z=0
637            z=x[DIM-1]
638            oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
639            sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
640        return data,sigma
641        
642      def getReferenceProperty(self, domain=None):
643        """
644        Returns the reference density Data object that was used to generate
645        the gravity/susceptibility anomaly data.
646        
647        :return: the density or susceptibility anomaly used to create the survey data.
648        :note: it can be assumed that in the first call the ``domain`` argument is present so the
649        actual anomaly data can be created. In subsequent calls this may not be true.
650        :note: method needs to be overwritten
651        """
652        raise NotImplementedError()      
653          
654    class SyntheticFeatureData(SyntheticDataBase):
655        """
656        uses a list of ``SourceFeature`` to define synthetic anomaly data
657        """
658        def __init__(self, datatype,
659                           features,
660                           DIM=2,
661                           number_of_elements=10,
662                           length=1*U.km,
663                           B_b=None,
664                           data_offset=0,
665                           full_knowledge=False,
666                           spherical=False):
667            """
668            :param datatype: data type indicator
669            :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
670            :param features: list of features. It is recommended that the features do entirely lay below surface.
671            :type features: ``list`` of ``SourceFeature``
672            :param DIM: spatial dimension
673            :type DIM: 2 or 3
674            :param number_of_elements: lateral number of elements in the region where data are collected
675            :type number_of_elements: ``int``
676            :param length: lateral extend of the region where data are collected
677            :type length: ``float``
678            :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
679            :type B_b: ``list`` of ``Scalar``
680            :param data_offset: offset of the data collection region from the surface
681            :type data_offset: ``float``
682            :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
683            :type full_knowledge: ``Bool``
684            :param spherical: if ``True`` sperical coordinates are used (ignored)
685            :type spherical: ``Bool``
686            """
687            super(SyntheticFeatureData,self).__init__(
688                                     datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,
689                                     length=length, B_b=B_b,
690                                     data_offset=data_offset,
691                                     full_knowledge=full_knowledge,
692                                     spherical=spherical)
693          self._features = features          self._features = features
         self.__origin = [0]*(DIM-1)  
         self.__nPts = [NE]*(DIM-1)  
         self.__delta = [float(l)/NE]*(DIM-1)  
         self.DIM=DIM  
         self.NE=NE  
         self.l=l  
694    
     def getDataExtents(self):  
         return (list(self.__origin), list(self.__nPts), list(self.__delta))  
695    
696      def getDataType(self):      def getReferenceProperty(self, domain=None):
         return self.__datatype  
   
     def getReferenceDensity(self):  
697          """          """
698          Returns the reference density Data object that was used to generate      Returns the reference density Data object that was used to generate
699          the gravity anomaly data.      the gravity/susceptibility anomaly data.
700          """          """
701          return self._rho          if self._reference_data == None:
702                DIM=domain.getDim()
703                x=domain.getX()
704                k=0.
705                for f in self._features:
706                    m=f.getMask(x)
707                    k = k * (1-m) + f.getValue(x) * m
708                self._reference_data= k
709            return self._reference_data
710            
711    class SyntheticData(SyntheticDataBase):
712        """
713        defines synthetic  gravity/magnetic data based on harmonic property anomaly
714        
715            rho = mean + amplitude * sin(n_depth * pi /depth * z) * sin(n_length * pi /length * x - shift )
716            
717        for all x and z<=0. for z>0 rho = 0.        
718        """
719        def __init__(self, datatype,
720                           n_length=1,
721                           n_depth=1,
722                           shift=0.,
723                           amplitude=None,
724                           DIM=2,
725                           number_of_elements=10,
726                           length=1*U.km,
727                           B_b=None,
728                           data_offset=0,
729                           full_knowledge=False,
730                           spherical=False):
731            """
732            :param datatype: data type indicator
733            :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
734            :param n_length: number of oscillation in the anomaly data within the observation region.
735            :type n_length: ``int``
736            :param n_depth: number of oscillation in the anomaly data below surface
737            :param shift: lateral phase shift in the  anomaly data
738            :param amplitude: data amplitude. Default value is 200 *U.kg/U.m**3 for gravity and 0.1 for magnetic data.
739            :param features: list of features. It is recommended that the features do entirely lay below surface.
740            :type features: ``list`` of ``SourceFeature``
741            :param DIM: spatial dimension
742            :type DIM: 2 or 3
743            :param number_of_elements: lateral number of elements in the region where data are collected
744            :type number_of_elements: ``int``
745            :param length: lateral extend of the region where data are collected
746            :type length: ``float``
747            :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
748            :type B_b: ``list`` of ``Scalar``
749            :param data_offset: offset of the data collection region from the surface
750            :type data_offset: ``float``
751            :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
752            :type full_knowledge: ``Bool``
753            :param spherical: if ``True`` sperical coordinates are used (ignored)
754            :type spherical: ``Bool``
755            """      
756            super(SyntheticData,self).__init__(
757                                     datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,
758                                     length=length, B_b=B_b,
759                                     data_offset=data_offset,
760                                     full_knowledge=full_knowledge,
761                                     spherical=spherical)
762            self.__n_length = n_length
763            self.__n_depth = n_depth
764            self.__shift = shift
765            if amplitude == None:
766            if datatype == DataSource.GRAVITY:
767                amplitude = 200 *U.kg/U.m**3
768            else:
769                 amplitude =0.1
770            self.__amplitude = amplitude
771    
     def getReferenceSusceptibility(self):  
         """  
         Returns the reference magnetic susceptibility Data objects that was  
         used to generate the magnetic field data.  
         """  
         return self._k  
772    
     def getSurveyData(self, domain, origin, NE, spacing):  
         pde=LinearSinglePDE(domain)  
         G=U.Gravitational_Constant  
         m_psi_ref=0.  
         x=domain.getX()  
         for i in range(self.DIM):  
             m_psi_ref=m_psi_ref + whereZero(x[i]-inf(x[i])) \  
                     + whereZero(x[i]-sup(x[i]))  
773    
774          if self.getDataType()==DataSource.GRAVITY:      def getReferenceProperty(self, domain=None):
775              rho_ref=0.          """
776              for f in self._features:          Returns the reference density Data object that was used to generate
777                  m=f.getMask(x)          the gravity anomaly data.
778                  rho_ref = rho_ref * (1-m) + f.getValue(x) * m          """
779              self._rho=rho_ref          if self._reference_data == None:
780              pde.setValue(A=kronecker(domain), Y=-4*np.pi*G*rho_ref, q=m_psi_ref)              DIM=domain.getDim()
781          else:              x=domain.getX()
782              k_ref=0.              # set the reference data
783              for f in self._features:              z=x[DIM-1]
784                  m=f.getMask(x)              k=sin(self.__n_depth * np.pi/inf(z) * z) * whereNegative(z) * self.__amplitude
785                  k_ref = k_ref * (1-m) + f.getValue(x) * m              for i in xrange(DIM-1):
786              self._k=k_ref             k*= sin(self.__n_length * np.pi /self.__length * x[i] - shift )
787              B_b = self.getBackgroundMagneticField()              self._reference_data= k
788              pde.setValue(A=kronecker(domain), X=(1+k_ref)*B_b, q=m_psi_ref)          return self._reference_data
         pde.setSymmetryOn()  
         psi_ref=pde.getSolution()  
         del pde  
         if self.getDataType()==DataSource.GRAVITY:  
             data = -grad(psi_ref, ReducedFunction(domain))  
         else:  
             data = (1+self._k)*B_b-grad(psi_ref, ReducedFunction(domain))  
789    
         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  
   
     def getBackgroundMagneticField(self):  
         #FIXME:  
         latitude=-28.0  
         theta = (90-latitude)/180.*np.pi  
         B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)  
         B_theta= B_0 * sin(theta)  
         B_r= 2 * B_0 * cos(theta)  
         if self.DIM<3:  
             return np.array([0.,  -B_r])  
         else:  
             return np.array([-B_theta, 0.,  -B_r])  
790    

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

  ViewVC Help
Powered by ViewVC 1.1.26