/[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 4053 by caltinay, Wed Oct 31 07:27:35 2012 UTC revision 4121 by gross, Wed Dec 19 00:24:50 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','UBCDataSource','ERSDataSource','SyntheticDataSource','SmoothAnomaly']  __all__ = ['simpleGeoMagneticFluxDensity', 'DataSource','ErMapperData', 'SyntheticDataBase' , 'SyntheticFeatureData', 'SyntheticData','SmoothAnomaly']
26    
27  import logging  import logging
28  import numpy as np  import numpy as np
29  from esys.escript import ReducedFunction, Scalar  from esys.escript import ReducedFunction
30  from esys.escript.linearPDEs import LinearSinglePDE  from esys.escript.linearPDEs import LinearSinglePDE
31  from esys.escript.util import *  from esys.escript.util import *
32  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
# Line 34  from esys.ripley import Brick, Rectangle Line 34  from esys.ripley import Brick, Rectangle
34    
35  try:  try:
36      from scipy.io.netcdf import netcdf_file      from scipy.io.netcdf import netcdf_file
37      __all__ += ['NetCDFDataSource']      __all__ += ['NetCdfData']
38  except:  except:
39      pass      pass
40    
# 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      """      """
90      A class that provides survey data for the inversion process.      A class that provides survey data for the inversion process.
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))
98    
99      def __init__(self):      def __init__(self):
100          """          """
101          Constructor. Sets some defaults and initializes logger.          Constructor. Sets some defaults and initializes logger.
102          """          """
         self._constrainBottom=False  
         self._constrainSides=True  
         self._domain=None  
         self.__set_density_mask=None  
         self.__set_susceptibility_mask=None  
         self.setPadding()  
103          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
104            self.__subsampling_factor=1
105      def setPadding(self, pad_x=10, pad_y=10):          self.__background_magnetic_field=None
         """  
         Sets the amount of padding around the dataset. If `pad_x`/`pad_y`  
         is >=1 the value is treated as number of elements to be added to the  
         domain (per side).  
         If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the  
         dataset size. For example, calling ``setPadding(3, 0.1)`` to a data  
         source with size 10x20 will result in the padded data set size  
         16x24 (10+2*3, 20*(1+2*0.1))  
   
         :param pad_x: Padding per side in x direction (default: 10 elements)  
         :type pad_x: ``int`` or ``float``  
         :param pad_y: Padding per side in y direction (default: 10 elements).  
                       This value is only used for 3-dimensional datasets  
         :type pad_y: ``int`` or ``float``  
         """  
         self._pad_x=pad_x  
         self._pad_y=pad_y  
   
     def setConstraints(self, bottom=False, sides=True):  
         """  
         If `bottom` is True, then the density mask will be set to 1 in the  
         padding area at the bottom of the domain. By default this area is  
         unconstrained. Similarly, if `sides` is True (default) then the  
         horizontal padding area is constrained, otherwise not.  
   
         :param bottom: Whether to constrain the density at the bottom of the  
                        domain  
         :type bottom: ``bool``  
         :param sides: Whether to constrain the density in the padding area  
                       surrounding the data  
         :type sides: ``bool``  
         """  
         self._constrainBottom=bottom  
         self._constrainSides=sides  
   
     def getDomain(self):  
         """  
         Returns a domain that spans the data area plus padding.  
         The domain is created the first time this method is called, subsequent  
         calls return the same domain so anything that affects the domain  
         (such as padding) needs to be set beforehand.  
   
         :return: The escript domain for this data source.  
         :rtype: `esys.escript.Domain`  
         """  
         if self._domain is None:  
             self._domain=self._createDomain(self._pad_x, self._pad_y)  
         return self._domain  
   
     def getSetDensityMask(self):  
         """  
         Returns the density mask data object, where mask has value 1 in the  
         padding area, 0 elsewhere.  
   
         :return: The mask for the density.  
         :rtype: `esys.escript.Data`  
         """  
         return self.__set_density_mask  
   
     def setSetDensityMask(self, mask):  
         """  
         Sets the density mask to use.  
         """  
         self.__set_density_mask=mask  
   
     def getSetSusceptibilityMask(self):  
         """  
         Returns the susceptibility mask data object, where mask has value 1  
         in the padding area, 0 elsewhere.  
   
         :return: The mask for the susceptibility.  
         :rtype: `esys.escript.Data`  
         """  
         return self.__set_susceptibility_mask  
   
     def setSetSusceptibilityMask(self, mask):  
         """  
         Sets the susceptibility mask to use.  
         """  
         self.__set_susceptibility_mask=mask  
   
     def getGravityAndStdDev(self):  
         """  
         Returns the gravity anomaly and standard deviation data objects as a  
         tuple. This method must be implemented in subclasses that supply  
         gravity data.  
         """  
         raise NotImplementedError  
   
     def getMagneticFieldAndStdDev(self):  
         """  
         Returns the magnetic field and standard deviation data objects as a  
         tuple. This method must be implemented in subclasses that supply  
         magnetic data.  
         """  
         raise NotImplementedError  
   
     def getBackgroundMagneticField(self):  
         """  
         Returns the background magnetic field. This method must be  
         implemented in subclasses that supply magnetic data.  
         """  
         raise NotImplementedError  
106    
107      def getDataExtents(self):      def getDataExtents(self):
108          """          """
# Line 216  class DataSource(object): Line 116  class DataSource(object):
116          """          """
117          raise NotImplementedError          raise NotImplementedError
118    
119      def getVerticalExtents(self):      def getDataType(self):
120          """          """
121          returns a tuple ``(z0, nz, dz)``, where          Returns the type of survey data managed by this source.
122            Subclasses must return `GRAVITY` or `MAGNETIC` as appropriate.
         - ``z0`` = minimum z coordinate (origin)  
         - ``nz`` = number of elements in z direction  
         - ``dz`` = spacing of elements (= cell size in z)  
   
         This method must be implemented in subclasses.  
123          """          """
124          raise NotImplementedError          raise NotImplementedError
125    
126      def _addPadding(self, pad_x, pad_y, NE, l, origin):      def getSurveyData(self, domain, origin, NE, spacing):
127          """          """
128          Helper method that computes new number of elements, length and origin          This method is called by the `DomainBuilder` to retrieve the survey
129          after adding padding to the input values.          data as `Data` objects on the given domain.
130            Subclasses should return one or more `Data` objects with survey data
131          :param pad_x: Number of elements or fraction of padding in x direction          interpolated on the given ripley domain. The exact return type
132          :type pad_x: ``int`` or ``float``          depends on the type of data.
133          :param pad_y: Number of elements or fraction of padding in y direction  
134          :type pad_y: ``int`` or ``float``          :param domain: the escript domain to use
135          :param NE: Initial number of elements          :type domain: `esys.escript.Domain`
136          :type NE: ``tuple`` or ``list``          :param origin: the origin coordinates of the domain
         :param l: Initial side lengths  
         :type l: ``tuple`` or ``list``  
         :param origin: Initial origin  
137          :type origin: ``tuple`` or ``list``          :type origin: ``tuple`` or ``list``
138          :return: tuple with three elements ``(NE_padded, l_padded, origin_padded)``,          :param NE: the number of domain elements in each dimension
139                   which are lists of the updated input parameters          :type NE: ``tuple`` or ``list``
140          """          :param spacing: the cell sizes (node spacing) in the domain
141          DIM=len(NE)          :type spacing: ``tuple`` or ``list``
         frac=[0.]*(DIM-1)+[0]  
         # padding is applied to each side so multiply by 2 to get the total  
         # amount of padding per dimension  
         if pad_x>0 and pad_y<1:  
             frac[0]=2.*pad_x  
         elif pad_x>=1:  
             frac[0]=2.*pad_x/float(NE[0])  
         if DIM>2:  
             if pad_y>0 and pad_y<1:  
                 frac[1]=2.*pad_y  
             elif pad_y>=1:  
                 frac[1]=2.*pad_y/(float(NE[1]))  
   
         # calculate new number of elements  
         NE_new=[int(NE[i]*(1+frac[i])) for i in range(DIM)]  
         NEdiff=[NE_new[i]-NE[i] for i in range(DIM)]  
         spacing=[l[i]/NE[i] for i in range(DIM)]  
         l_new=[NE_new[i]*spacing[i] for i in range(DIM)]  
         origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in range(DIM)]  
         return NE_new, l_new, origin_new  
   
     def _interpolateOnDomain(self, data):  
         """  
         Helper method that interpolates data arrays onto the domain.  
         Currently this works like a nearest neighbour mapping, i.e. values  
         are directly inserted into data objects at closest location.  
         """  
         dom=self.getDomain()  
         dim=dom.getDim()  
         # determine number of values required per element  
         DPP=Scalar(0., ReducedFunction(dom)).getNumberOfDataPoints()  
         for i in range(dim):  
             DPP=DPP/self._dom_NE[i]  
         DPP=int(DPP)  
   
         # idx_mult.dot([x,y,z]) = flat index into data object  
         idx_mult=np.array([DPP]+self._dom_NE[:dim-1]).cumprod()  
   
         # separate data arrays and coordinates  
         num_arrays=len(data[0])-dim  
         arrays=[]  
         for i in range(num_arrays):  
             d=Scalar(0., ReducedFunction(dom))  
             d.expand()  
             arrays.append(d)  
   
         for entry in data:  
             index=[int((entry[i]-self._dom_origin[i])/self._spacing[i]) for i in range(dim)]  
             index=int(idx_mult.dot(index))  
             for i in range(num_arrays):  
                 for p in range(DPP):  
                     arrays[i].setValueOfDataPoint(index+p, entry[dim+i])  
   
         return arrays  
   
     def _createDomain(self, padding_x, padding_y):  
         """  
         Creates and returns an escript domain that spans the entire area of  
         available data plus a buffer zone. This method is called only once  
         the first time `getDomain()` is invoked and may be overwritten if  
         required.  
   
         :return: The escript domain for this data source.  
         :rtype: `esys.escript.Domain`  
         """  
         X0, NX, DX = self.getDataExtents()  
         z0, nz, dz = self.getVerticalExtents()  
   
         # number of elements (without padding)  
         NE = [NX[0], NX[1], nz]  
   
         # origin of domain (without padding)  
         origin = [X0[0], X0[1], z0]  
         origin = [np.round(oi) for oi in origin]  
   
         # cell size / point spacing  
         self._spacing = DX+[dz]  
         self._spacing = [float(np.round(si)) for si in self._spacing]  
   
         # length of domain (without padding)  
         l = [NE[i]*self._spacing[i] for i in range(len(NE))]  
   
         # now add padding to the values  
         NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \  
                 NE, l, origin)  
   
         # number of padding elements per side  
         NE_pad=[(NE_new[i]-NE[i])//2 for i in range(3)]  
   
         self._dom_NE_pad = NE_pad  
         self._dom_len = l_new  
         self._dom_NE = NE_new  
         self._dom_origin = origin_new  
         lo=[(origin_new[i], origin_new[i]+l_new[i]) for i in range(3)]  
         dom=Brick(*self._dom_NE, l0=lo[0], l1=lo[1], l2=lo[2])  
         # ripley may internally adjust NE and length, so recompute  
         self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(3)]  
         self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(3)]  
         x=dom.getX()-[self._dom_origin[i]+NE_pad[i]*self._spacing[i] for i in range(3)]  
         mask=wherePositive(dom.getX()[2])  
   
         # prepare density mask (=1 at padding area, 0 else)  
         if self._constrainSides:  
             for i in range(2):  
                 mask=mask + whereNegative(x[i]) + \  
                         wherePositive(x[i]-l_new[i]+2*NE_pad[i]*self._spacing[i])  
   
         if self._constrainBottom:  
             mask = mask + whereNonPositive(x[2])  
         self.setSetDensityMask(wherePositive(mask))  
         self.setSetSusceptibilityMask(wherePositive(mask))  
   
         self.logger.debug("Domain size: %d x %d x %d elements"%(self._dom_NE[0],self._dom_NE[1],self._dom_NE[2]))  
         self.logger.debug("     length: %g x %g x %g"%(self._dom_len[0],self._dom_len[1],self._dom_len[2]))  
         self.logger.debug("     origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))  
   
         return dom  
   
   
 ##############################################################################  
 class UBCDataSource(DataSource):  
     def __init__(self, meshfile, gravfile, topofile=None):  
         super(UBCDataSource,self).__init__()  
         self.__meshfile=meshfile  
         self.__gravfile=gravfile  
         self.__topofile=topofile  
         self.__readMesh()  
   
     def __readMesh(self):  
         meshdata=open(self.__meshfile).readlines()  
         numDataPoints=meshdata[0].split()  
         origin=meshdata[1].split()  
         self.__nPts=map(int, numDataPoints)  
         self.__origin=map(float, origin)  
         self.__delta=[float(X.split('*')[1]) for X in meshdata[2:]]  
         # vertical data is upside down  
         self.__origin[2]-=(self.__nPts[2]-1)*self.__delta[2]  
         self.logger.debug("Data Source: %s (mesh file: %s)"%(self.__gravfile, self.__meshfile))  
   
     def getDataExtents(self):  
         """  
         returns ( (x0, y0), (nx, ny), (dx, dy) )  
         """  
         return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])  
   
     def getVerticalExtents(self):  
         """  
         returns (z0, nz, dz)  
142          """          """
143          return (self.__origin[2], self.__nPts[2], self.__delta[2])          raise NotImplementedError
   
     #def getSetDensityMask(self):  
     #    topodata=self.__readTopography()  
     #    mask=self._interpolateOnDomain(topodata)  
     #    mask=wherePositive(self.getDomain().getX()[2]-mask[0])  
     #    return mask  
   
     def getGravityAndStdDev(self):  
         gravlist=self.__readGravity() # x,y,z,g,s  
         g_and_sigma=self._interpolateOnDomain(gravlist)  
         return g_and_sigma[0]*[0,0,1], g_and_sigma[1]  
   
     def __readTopography(self):  
         f=open(self.__topofile)  
         n=int(f.readline())  
         topodata=np.zeros((n,3))  
         for i in range(n):  
             x=f.readline().split()  
             x=map(float, x)  
             topodata[i]=x  
         f.close()  
         return topodata  
   
     def __readGravity(self):  
         f=open(self.__gravfile)  
         n=int(f.readline())  
         gravdata=np.zeros((n,5))  
         for i in range(n):  
             x=f.readline().split()  
             x=map(float, x) # x, y, z, anomaly in mGal, stddev  
             # convert gravity anomaly units to m/s^2 and rescale error  
             x[3]*=-1e-5  
             x[4]*=1e-5  
             gravdata[i]=x  
         f.close()  
         return gravdata  
   
 ##############################################################################  
 class NetCDFDataSource(DataSource):  
     """  
     Data Source for gridded netCDF data that use CF/COARDS conventions.  
     """  
     def __init__(self, gravfile=None, magfile=None, topofile=None,  
                        depth=40000.,  
                        air_layer=10000.,  
                        vertical_cells=25,  
                        alt_of_data=0.):  
         """  
         :param gravfile: file with gravity data in netcdf format  
         :type gravfile: ``str``  
         :param magfile: file with magnetic datafile in netcdf format  
         :type magfile: ``str``  
         :param topofile: file with topography data in netcdf format  
         :type topofile: ``str``  
         :param depth: depth below surface filled with rock  
         :type depth: ``float``  
         :param air_layer: height above surface filled with air. It is  
                           assumed that that density correction is zero  
                           in this region.  
         :type air_layer: ``float``  
         :param vertical_cells: number of vertical cells  
         :type vertical_cells: ``int``  
         :param alt_of_data: altitude of measurements above ground in meters  
         :type alt_of_data: ``float``  
         """  
         super(NetCDFDataSource,self).__init__()  
         self.__topofile=topofile  
         self.__gravfile=gravfile  
         self.__magfile=magfile  
         self.__determineExtents((-depth,air_layer,vertical_cells))  
         self.__altOfData=alt_of_data  
   
     def __determineExtents(self, ve):  
         self.logger.debug("Data Source: %s"%self.__gravfile)  
         f=netcdf_file(self.__gravfile, 'r')  
         NX=0  
         for n in ['lon','longitude','x']:  
             if n in f.dimensions:  
                 NX=f.dimensions[n]  
                 break  
         if NX==0:  
             raise RuntimeError("Could not determine extents of data")  
         NY=0  
         for n in ['lat','latitude','y']:  
             if n in f.dimensions:  
                 NY=f.dimensions[n]  
                 break  
         if NY==0:  
             raise RuntimeError("Could not determine extents of data")  
   
         # find longitude and latitude variables  
         lon_name=None  
         for n in ['lon','longitude']:  
             if n in f.variables:  
                 lon_name=n  
                 longitude=f.variables.pop(n)  
                 break  
         if lon_name is None:  
             raise RuntimeError("Could not determine longitude variable")  
         lat_name=None  
         for n in ['lat','latitude']:  
             if n in f.variables:  
                 lat_name=n  
                 latitude=f.variables.pop(n)  
                 break  
         if lat_name is None:  
             raise RuntimeError("Could not determine latitude variable")  
   
         # try to figure out gravity variable name  
         grav_name=None  
         if len(f.variables)==1:  
             grav_name=f.variables.keys()[0]  
         else:  
             for n in f.variables.keys():  
                 dims=f.variables[n].dimensions  
                 if (lat_name in dims) and (lon_name in dims):  
                     grav_name=n  
                     break  
         if grav_name is None:  
             raise RuntimeError("Could not determine gravity variable")  
   
         # try to determine value for unused data  
         if hasattr(f.variables[grav_name], 'missing_value'):  
             maskval = float(f.variables[grav_name].missing_value)  
         elif hasattr(f.variables[grav_name], '_FillValue'):  
             maskval = float(f.variables[grav_name]._FillValue)  
         else:  
             self.logger.debug("missing_value attribute not found, using default.")  
             maskval = 99999  
   
         # see if there is a wkt string to convert coordinates  
         try:  
             wkt_string=f.variables[grav_name].esri_pe_string  
         except:  
             wkt_string=None  
   
         # we don't trust actual_range & geospatial_lon_min/max since subset  
         # data does not seem to have these fields updated.  
         # Getting min/max from the arrays is obviously not very efficient but..  
         #lon_range=longitude.actual_range  
         #lat_range=latitude.actual_range  
         #lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]  
         #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]  
         self.latitude=latitude.data.min()  
         lon_range=longitude.data.min(),longitude.data.max()  
         lat_range=latitude.data.min(),latitude.data.max()  
         lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)  
         origin=[lon_range[0],lat_range[0],ve[0]]  
         lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0],ve[1]-ve[0]]  
   
         f.close()  
   
         self.__nPts=[NX, NY, ve[2]]  
         self.__origin=origin  
         # we are rounding to avoid interpolation issues  
         self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(3)]  
         self.__wkt_string=wkt_string  
         self.__lon=lon_name  
         self.__lat=lat_name  
         self.__grv=grav_name  
         self.__maskval=maskval  
144    
145      def getDataExtents(self):      def setSubsamplingFactor(self, f):
146          """          """
147          returns ( (x0, y0), (nx, ny), (dx, dy) )          Sets the data subsampling factor (default=1).
148            The factor is applied in all dimensions. For example a 2D dataset
149            with 300 x 150 data points will be reduced to 150 x 75 when a
150            subsampling factor of 2 is used.
151            This becomes important when adding data of varying resolution to
152            a `DomainBuilder`.
153          """          """
154          return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])          self.__subsampling_factor=f
155    
156      def getVerticalExtents(self):      def getSubsamplingFactor(self):
157          """          """
158          returns (z0, nz, dz)          Returns the subsampling factor that was set via `setSubsamplingFactor`
159            (see there).
160          """          """
161          return (self.__origin[2], self.__nPts[2], self.__delta[2])          return self.__subsampling_factor
   
     def getBackgroundMagneticField(self):  
         theta = (90-self.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)  
         return np.array([-B_theta, 0., -B_r])  
   
     def getMagneticFieldAndStdDev(self):  
         B,sigma=self.getGravityAndStdDev()  
         return 1e-3*B,1e-3*sigma  
   
     def getGravityAndStdDev(self):  
         nValues=self.__nPts[:2]+[1]  
         first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]  
         g=ripleycpp._readNcGrid(self.__gravfile, self.__grv,  
                 ReducedFunction(self.getDomain()),  
                 first, nValues, (), self.__maskval)  
         sigma=whereNonZero(g-self.__maskval)  
         g=g*1e-6  
         sigma=sigma*2e-6  
         return g*[0,0,1], sigma  
   
     def _readTopography(self):  
         f=netcdf_file(self.__topofile, 'r')  
         lon=None  
         for n in ['lon','longitude']:  
             if n in f.variables:  
                 lon=f.variables[n][:]  
                 break  
         if lon is None:  
             raise RuntimeError("Could not determine longitude variable")  
         lat=None  
         for n in ['lat','latitude']:  
             if n in f.variables:  
                 lat=f.variables[n][:]  
                 break  
         if lat is None:  
             raise RuntimeError("Could not determine latitude variable")  
         alt=None  
         for n in ['altitude','alt']:  
             if n in f.variables:  
                 alt=f.variables[n][:]  
                 break  
         if alt is None:  
             raise RuntimeError("Could not determine altitude variable")  
162    
         topodata=np.column_stack((lon,lat,alt))  
         f.close()  
         return topodata  
163    
164  ##############################################################################  ##############################################################################
165  class ERSDataSource(DataSource):  class ErMapperData(DataSource):
166      """      """
167      Data Source for ER Mapper raster data.      Data Source for ER Mapper raster data.
168      Note that this class only accepts a very specific type of ER Mapper data      Note that this class only accepts a very specific type of ER Mapper data
169      input and will raise an exception if other data is found.      input and will raise an exception if other data is found.
170      """      """
171      def __init__(self, headerfile, datafile=None, depth=40000.,      def __init__(self, datatype, headerfile, datafile=None, altitude=0.):
                  air_layer=10000., vertical_cells=25, alt_of_data=0.):  
172          """          """
173            :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`
174            :type datatype: ``int``
175          :param headerfile: ER Mapper header file (usually ends in .ers)          :param headerfile: ER Mapper header file (usually ends in .ers)
176          :type headerfile: ``str``          :type headerfile: ``str``
177          :param datafile: ER Mapper binary data file name. If not supplied the          :param datafile: ER Mapper binary data file name. If not supplied the
178                           name of the header file without '.ers' is assumed                           name of the header file without '.ers' is assumed
179          :type datafile: ``str``          :type datafile: ``str``
180          :param depth: depth below surface filled with rock          :param altitude: altitude of measurements above ground in meters
181          :type depth: ``float``          :type altitude: ``float``
         :param air_layer: height above surface filled with air. It is  
                           assumed that that density correction is zero  
                           in this region.  
         :type air_layer: ``float``  
         :param vertical_cells: number of vertical cells  
         :type vertical_cells: ``int``  
         :param alt_of_data: altitude of measurements above ground in meters  
         :type alt_of_data: ``float``  
182          """          """
183          super(ERSDataSource,self).__init__()          super(ErMapperData,self).__init__()
184          self.__headerfile=headerfile          self.__headerfile=headerfile
185          if datafile is None:          if datafile is None:
186              self.__datafile=headerfile[:-4]              self.__datafile=headerfile[:-4]
187          else:          else:
188              self.__datafile=datafile              self.__datafile=datafile
189          self.__readHeader((-depth, air_layer, vertical_cells))          self.__altitude=altitude
190          self.__altOfData=alt_of_data          self.__datatype=datatype
191            self.__readHeader()
192    
193      def __readHeader(self, ve):      def __readHeader(self):
194          self.logger.debug("Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))          self.logger.debug("Checking Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))
195          metadata=open(self.__headerfile, 'r').readlines()          metadata=open(self.__headerfile, 'r').readlines()
196          # parse metadata          # parse metadata
197          start=-1          start=-1
# Line 741  class ERSDataSource(DataSource): Line 275  class ERSDataSource(DataSource):
275          # data sets have origin in top-left corner so y runs top-down          # data sets have origin in top-left corner so y runs top-down
276          self.__dataorigin=[originX, originY]          self.__dataorigin=[originX, originY]
277          originY-=(NY-1)*spacingY          originY-=(NY-1)*spacingY
278          self.__maskval=maskval          self.__delta = [spacingX, spacingY]
279          spacingZ=np.round(float(ve[1]-ve[0])/ve[2])          self.__maskval = maskval
280          self.__delta = [spacingX, spacingY, spacingZ]          self.__nPts = [NX, NY]
281          self.__nPts = [NX, NY, ve[2]]          self.__origin = [originX, originY]
282          self.__origin = [originX, originY, ve[0]]          if self.__datatype == self.GRAVITY:
283                self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
284                self.__scalefactor = 1e-6
285            else:
286                self.logger.info("Assuming magnetic data units are 'nT'.")
287                self.__scalefactor = 1e-9
288    
289      def getDataExtents(self):      def getDataExtents(self):
290          """          """
291          returns ( (x0, y0), (nx, ny), (dx, dy) )          returns ( (x0, y0), (nx, ny), (dx, dy) )
292          """          """
293          return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])          return (list(self.__origin), list(self.__nPts), list(self.__delta))
294    
295        def getDataType(self):
296            return self.__datatype
297    
298        def getSurveyData(self, domain, origin, NE, spacing):
299            nValues=self.__nPts
300            # determine base location of this dataset within the domain
301            first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
302            if domain.getDim()==3:
303                first.append(int((self.__altitude-origin[2])/spacing[2]))
304                nValues=nValues+[1]
305    
306            data=ripleycpp._readBinaryGrid(self.__datafile,
307                    ReducedFunction(domain),
308                    first, nValues, (), self.__maskval)
309            sigma = whereNonZero(data-self.__maskval)
310            data = data*self.__scalefactor
311            sigma = sigma * 2. * self.__scalefactor
312            return data, sigma
313    
314    
315    ##############################################################################
316    class NetCdfData(DataSource):
317        """
318        Data Source for gridded netCDF data that use CF/COARDS conventions.
319        """
320        def __init__(self, datatype, filename, altitude=0.):
321            """
322            :param filename: file name for survey data in netCDF format
323            :type filename: ``str``
324            :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`
325            :type datatype: ``int``
326            :param altitude: altitude of measurements in meters
327            :type altitude: ``float``
328            """
329            super(NetCdfData,self).__init__()
330            self.__filename=filename
331            if not datatype in [self.GRAVITY,self.MAGNETIC]:
332                raise ValueError("Invalid value for datatype parameter")
333            self.__datatype=datatype
334            self.__altitude=altitude
335            self.__readMetadata()
336    
337        def __readMetadata(self):
338            self.logger.debug("Checking Data Source: %s"%self.__filename)
339            f=netcdf_file(self.__filename, 'r')
340            NX=0
341            for n in ['lon','longitude','x']:
342                if n in f.dimensions:
343                    NX=f.dimensions[n]
344                    break
345            if NX==0:
346                raise RuntimeError("Could not determine extents of data")
347            NY=0
348            for n in ['lat','latitude','y']:
349                if n in f.dimensions:
350                    NY=f.dimensions[n]
351                    break
352            if NY==0:
353                raise RuntimeError("Could not determine extents of data")
354    
355            # find longitude and latitude variables
356            lon_name=None
357            for n in ['lon','longitude']:
358                if n in f.variables:
359                    lon_name=n
360                    longitude=f.variables.pop(n)
361                    break
362            if lon_name is None:
363                raise RuntimeError("Could not determine longitude variable")
364            lat_name=None
365            for n in ['lat','latitude']:
366                if n in f.variables:
367                    lat_name=n
368                    latitude=f.variables.pop(n)
369                    break
370            if lat_name is None:
371                raise RuntimeError("Could not determine latitude variable")
372    
373            # try to figure out data variable name
374            data_name=None
375            if len(f.variables)==1:
376                data_name=f.variables.keys()[0]
377            else:
378                for n in f.variables.keys():
379                    dims=f.variables[n].dimensions
380                    if (lat_name in dims) and (lon_name in dims):
381                        data_name=n
382                        break
383            if data_name is None:
384                raise RuntimeError("Could not determine data variable")
385    
386            # try to determine value for unused data
387            if hasattr(f.variables[data_name], 'missing_value'):
388                maskval = float(f.variables[data_name].missing_value)
389            elif hasattr(f.variables[data_name], '_FillValue'):
390                maskval = float(f.variables[data_name]._FillValue)
391            else:
392                self.logger.debug("missing_value attribute not found, using default.")
393                maskval = 99999
394    
395            # try to determine units of data - this is disabled for now
396            #if hasattr(f.variables[data_name], 'units'):
397            #   units=f.variables[data_name].units
398            if self.__datatype == self.GRAVITY:
399                self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
400                self.__scalefactor = 1e-6
401            else:
402                self.logger.info("Assuming magnetic data units are 'nT'.")
403                self.__scalefactor = 1e-9
404    
405            # see if there is a wkt string to convert coordinates
406            try:
407                wkt_string=f.variables[data_name].esri_pe_string
408            except:
409                wkt_string=None
410    
411            # we don't trust actual_range & geospatial_lon_min/max since subset
412            # data does not seem to have these fields updated.
413            # Getting min/max from the arrays is obviously not very efficient but..
414            #lon_range=longitude.actual_range
415            #lat_range=latitude.actual_range
416            #lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
417            #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
418            lon_range=longitude.data.min(),longitude.data.max()
419            lat_range=latitude.data.min(),latitude.data.max()
420            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]]
423            f.close()
424    
425            self.__nPts=[NX, NY]
426            self.__origin=[lon_range[0],lat_range[0]]
427            # we are rounding to avoid interpolation issues
428            self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(2)]
429            #self.__wkt_string=wkt_string
430            #self.__lon_name=lon_name
431            #self.__lat_name=lat_name
432            self.__data_name=data_name
433            self.__maskval=maskval
434    
435      def getVerticalExtents(self):      def getDataExtents(self):
436          """          """
437          returns (z0, nz, dz)          returns ( (x0, y0), (nx, ny), (dx, dy) )
438          """          """
439          return (self.__origin[2], self.__nPts[2], self.__delta[2])          return (list(self.__origin), list(self.__nPts), list(self.__delta))
440    
441      def getGravityAndStdDev(self):      def getDataType(self):
442          nValues=self.__nPts[:2]+[1]          return self.__datatype
443          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]  
444          g=ripleycpp._readBinaryGrid(self.__datafile,      def getSurveyData(self, domain, origin, NE, spacing):
445                  ReducedFunction(self.getDomain()),          nValues=self.__nPts
446                  first, nValues, (), self.__maskval)          # determine base location of this dataset within the domain
447          sigma=whereNonZero(g-self.__maskval)          first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
448          g=g*1e-6          if domain.getDim()==3:
449          sigma=sigma*2e-6              first.append(int((self.__altitude-origin[2])/spacing[2]))
450          return g*[0,0,1], sigma              nValues=nValues+[1]
451    
452            data=ripleycpp._readNcGrid(self.__filename, self.__data_name,
453                      ReducedFunction(domain), first, nValues, (), self.__maskval)
454            sigma=whereNonZero(data-self.__maskval)
455            data=data*self.__scalefactor
456            sigma=sigma * 2. * self.__scalefactor
457            return data, sigma
458    
459    
460  ##############################################################################  ##############################################################################
# Line 778  class SourceFeature(object): Line 464  class SourceFeature(object):
464      data source, for example a layer of a specific rock type or a simulated      data source, for example a layer of a specific rock type or a simulated
465      ore body.      ore body.
466      """      """
467      def getDensity(self):      def getValue(self):
468          """          """
469          Returns the density for the area covered by mask. It can be constant          Returns the value for the area covered by mask. It can be constant
470          or a data object with spatial dependency.          or a data object with spatial dependency.
471          """          """
472          raise NotImplementedError          raise NotImplementedError
# Line 793  class SourceFeature(object): Line 479  class SourceFeature(object):
479          raise NotImplementedError          raise NotImplementedError
480    
481  class SmoothAnomaly(SourceFeature):  class SmoothAnomaly(SourceFeature):
482      def __init__(self, lx, ly, lz, x, y, depth, rho_inner=None, rho_outer=None, k_inner=None, k_outer=None):      def __init__(self, lx, ly, lz, x, y, depth, v_inner=None, v_outer=None):
483          self.x=x          self.x=x
484          self.y=y          self.y=y
485          self.lx=lx          self.lx=lx
486          self.ly=ly          self.ly=ly
487          self.lz=lz          self.lz=lz
488          self.depth=depth          self.depth=depth
489          self.rho_inner=rho_inner          self.v_inner=v_inner
490          self.rho_outer=rho_outer          self.v_outer=v_outer
491          self.k_inner=k_inner          self.value=None
         self.k_outer=k_outer  
         self.rho=None  
         self.k=None  
492          self.mask=None          self.mask=None
493    
494      def getDensity(self,x):      def getValue(self,x):
495          if self.rho is None:          if self.value is None:
496              if self.rho_outer is None or self.rho_inner is None:              if self.v_outer is None or self.v_inner is None:
497                  self.rho=0                  self.value=0
             else:  
                 DIM=x.getDomain().getDim()  
                 alpha=-log(abs(self.rho_outer/self.rho_inner))*4  
                 rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)  
                 rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)  
                 self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))  
                 if self.rho_inner<0: self.rho=-self.rho  
   
         return self.rho  
   
     def getSusceptibility(self,x):  
          if self.k is None:  
             if self.k_outer is None or self.k_inner is None:  
                 self.k=0  
498              else:              else:
499                  DIM=x.getDomain().getDim()                  DIM=x.getDomain().getDim()
500                  alpha=-log(abs(self.k_outer/self.k_inner))*4                  alpha=-log(abs(self.v_outer/self.v_inner))*4
501                  k=exp(-alpha*((x[0]-self.x)/self.lx)**2)                  value=exp(-alpha*((x[0]-self.x)/self.lx)**2)
502                  k=k*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)                  value=value*exp(-alpha*((x[DIM-1]+self.depth)/self.lz)**2)
503                  self.k=maximum(abs(self.k_outer), abs(self.k_inner*k))                  self.value=maximum(abs(self.v_outer), abs(self.v_inner*value))
504                  if self.k_inner<0: self.k=-self.k                  if self.v_inner<0: self.value=-self.value
505    
506           return self.k          return self.value
507    
508      def getMask(self, x):      def getMask(self, x):
509          DIM=x.getDomain().getDim()          DIM=x.getDomain().getDim()
510          m=whereNonNegative(x[DIM-1]-(sup(x[DIM-1])-self.depth-self.lz/2)) * whereNonPositive(x[DIM-1]-(sup(x[DIM-1])-self.depth+self.lz/2)) \          m=whereNonNegative(x[DIM-1]+self.depth+self.lz/2) * whereNonPositive(x[DIM-1]+self.depth-self.lz/2) \
511              *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))              *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
512          if DIM>2:          if DIM>2:
513              m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))              m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
# Line 846  class SmoothAnomaly(SourceFeature): Line 515  class SmoothAnomaly(SourceFeature):
515          return m          return m
516    
517  ##############################################################################  ##############################################################################
518  class SyntheticDataSource(DataSource):  class SyntheticDataBase(DataSource):
519      def __init__(self, DIM, NE, l, h, features, latitude=-32.):    """
520          super(SyntheticDataSource,self).__init__()    Base class to define reference data based on a given property distribution (density or
521      susceptibility). Data are collected from a square region of vertical extend `length` on a
522      grid with ``number_of_elements`` cells in each direction.
523      
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))+whereZero(z-inf(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.DIM=DIM  
         self.NE=NE  
         self.l=l  
         self.h=h  
         self.latitude=latitude  
   
     def _createDomain(self, padding_x, padding_y):  
         NE_H=self.NE  
         NE_L=int((self.l/self.h)*NE_H+0.5)  
         l=[self.l]*(self.DIM-1)+[self.h]  
         NE=[NE_L]*(self.DIM-1)+[NE_H]  
         origin=[0.]*self.DIM  
         NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \  
                 NE, l, origin)  
   
         self.NE=NE_new  
         self.l=l_new[0]  
         self.h=l_new[self.DIM-1]  
   
         self.logger.debug("Data Source: synthetic with %d features"%len(self._features))  
         if self.DIM==2:  
             dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])  
             self._x = dom.getX() + origin_new  
             self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))  
             self.logger.debug("     length: %g x %g"%(l_new[0],l_new[1]))  
             self.logger.debug("     origin: %g x %g"%(origin_new[0],origin_new[1]))  
         else:  
             dom = Brick(n0=NE_new[0], n1=NE_new[1], n2=NE_new[2], l0=l_new[0], l1=l_new[1], l2=l_new[2])  
             self._x = dom.getX() + origin_new  
             self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))  
             self.logger.debug("     length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))  
             self.logger.debug("     origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))  
   
         dz=l_new[self.DIM-1]/NE_new[self.DIM-1]  
         self._g_mask=wherePositive(dom.getX()[0]-origin_new[0]) \  
                 * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \  
                 * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \  
                 * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))  
   
         self._B_mask=self._g_mask  
   
         mask=whereNegative(self._x[self.DIM-1]) + \  
                 wherePositive(self._x[self.DIM-1]-l[self.DIM-1])  
         for i in range(self.DIM-1):  
             mask+= whereNegative(self._x[i]) +  wherePositive(self._x[i]-l[i])  
         self.setSetDensityMask(wherePositive(mask))  
         self.setSetSusceptibilityMask(wherePositive(mask))  
   
         rho_ref=0.  
         k_ref=0  
         for f in self._features:  
             m=f.getMask(self._x)  
             rho_ref = rho_ref * (1-m) + f.getDensity(self._x) * m  
             k_ref = k_ref * (1-m) + f.getSusceptibility(self._x) * m  
         self._rho=rho_ref  
         self._k=k_ref  
   
         return dom  
   
   
     def getReferenceDensity(self):  
         return self._rho  
     def getReferenceSusceptibility(self):  
         return self._k  
   
     def getGravityAndStdDev(self):  
         pde=LinearSinglePDE(self.getDomain())  
         G=U.Gravitational_Constant  
         m_psi_ref=0.  
         for i in range(self.DIM):  
             m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \  
                     + whereZero(self._x[i]-sup(self._x[i]))  
   
         pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)  
         pde.setSymmetryOn()  
         psi_ref=pde.getSolution()  
         del pde  
         g=-grad(psi_ref)  
         sigma=self._g_mask  
         return g,sigma  
   
     def getMagneticFieldAndStdDev(self):  
         pde=LinearSinglePDE(self.getDomain())  
         B_b=self.getBackgroundMagneticField()  
         DIM=self.getDomain().getDim()  
         m_psi_ref=0.  
         for i in range(self.DIM):  
             m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \  
                     + whereZero(self._x[i]-sup(self._x[i]))  
         pde.setValue(A=kronecker(self.getDomain()), X=(1+self._k)*B_b, q=m_psi_ref)  
         pde.setSymmetryOn()  
         psi_ref=pde.getSolution()  
         del pde  
         B= (1+self._k) * B_b -grad(psi_ref)  
         sigma=self._B_mask  
         return B,sigma  
695    
696      def getBackgroundMagneticField(self):  
697          theta = (90-self.latitude)/180.*np.pi      def getReferenceProperty(self, domain=None):
698          B_0=U.Mu_0  * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi *  U.R_Earth**3)          """
699          B_theta= B_0 * sin(theta)      Returns the reference density Data object that was used to generate
700          B_r= 2 * B_0 * cos(theta)      the gravity/susceptibility anomaly data.
701          DIM=self.getDomain().getDim()          """
702          if DIM<3:          if self._reference_data == None:
703              return np.array([0.,  -B_r])              DIM=domain.getDim()
704          else:              x=domain.getX()
705              return np.array([-B_theta, 0.,  -B_r])              k=0.
706                for f in self._features:
707                    m=f.getMask(x)
708                    k = k * (1-m) + f.getValue(x) * m
709                self._reference_data= k
710            return self._reference_data
711            
712    class SyntheticData(SyntheticDataBase):
713        """
714        defines synthetic  gravity/magnetic data based on harmonic property anomaly
715        
716            rho = mean + amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * x - shift )
717            
718        for all x and z<=0. for z>0 rho = 0.        
719        """
720        def __init__(self, datatype,
721                           n_length=1,
722                           n_depth=1,
723                           depth_offset=0.,
724                           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    
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               m=whereNonNegative(x_i)*whereNonNegative(self.length-x_i)
795               k*= sin(self.__n_length*np.pi /self.length * x_i)*m
796                self._reference_data= k
797            return self._reference_data
798    
799    

Legend:
Removed from v.4053  
changed lines
  Added in v.4121

  ViewVC Help
Powered by ViewVC 1.1.26