/[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 3947 by caltinay, Wed Aug 22 23:19:10 2012 UTC revision 4121 by gross, Wed Dec 19 00:24:50 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16    """Data readers/providers for inversions"""
17    
18  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
20  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
21  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
22  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
23  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
24    
25  __all__ = ['DataSource','UBCDataSource','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 *  from esys.escript import ReducedFunction
30  from esys.escript.linearPDEs import *  from esys.escript.linearPDEs import LinearSinglePDE
31    from esys.escript.util import *
32  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
33    from esys.ripley import Brick, Rectangle, ripleycpp
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    
41    def LatLonToUTM(lon, lat, wkt_string=None):
42        """
43        Converts one or more longitude,latitude pairs to the corresponding x,y
44        coordinates in the Universal Transverse Mercator projection.
45    
46        :note: If the ``pyproj`` module is not installed a warning is printed and
47               the input values are scaled by a constant and returned.
48        :note: If `wkt_string` is not given or invalid or the ``gdal`` module is
49               not available to convert the string, then the input values are
50               assumed to be given in the Clarke 1866 projection.
51    
52        :param lon: longitude value(s)
53        :type lon: `float`, `list`, `tuple`, or ``numpy.array``
54        :param lat: latitude value(s)
55        :type lat: `float`, `list`, `tuple`, or ``numpy.array``
56        :rtype: ``numpy.array``
57        """
58    
59        # not really optimal: if pyproj is not installed we return the input
60        # values scaled by a constant.
61        try:
62            import pyproj
63        except:
64            print("Warning, pyproj not available. Domain extents will be wrong")
65            return np.array(lon)*1000., np.array(lat)*1000.
66    
67        # determine UTM zone from the input data
68        zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
69        try:
70            import osgeo.osr
71            srs = osgeo.osr.SpatialReference()
72            srs.ImportFromWkt(wkt_string)
73            p_src = pyproj.Proj(srs.ExportToProj4())
74        except:
75            p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
76        # we assume southern hemisphere here
77        p_dest = pyproj.Proj('+proj=utm +zone=%d +south +units=m'%zone)
78        x,y=pyproj.transform(p_src, p_dest, lon, lat)
79        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.
92        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.
94        For other setups override the methods as required.
95      """      """
96      # this is currently specific to gravity inversion and should be generalised  
97        GRAVITY, MAGNETIC = list(range(2))
98    
99      def __init__(self):      def __init__(self):
100          """          """
101            Constructor. Sets some defaults and initializes logger.
102          """          """
         self._domain=None  
         self._pad_l=0.1  
         self._pad_h=0.1  
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            self.__background_magnetic_field=None
106    
107      def _addPadding(self, pad_l, pad_h, NE, l, origin):      def getDataExtents(self):
108          """          """
109          Helper method that computes new number of elements, length and origin          returns a tuple of tuples ``( (x0, y0), (nx, ny), (dx, dy) )``, where
110          after adding padding to the input values.  
111            - ``x0``, ``y0`` = coordinates of data origin
112            - ``nx``, ``ny`` = number of data points in x and y
113            - ``dx``, ``dy`` = spacing of data points in x and y
114    
115            This method must be implemented in subclasses.
116          """          """
117          DIM=len(NE)          raise NotImplementedError
         frac=[0.]*DIM  
         # padding is applied to each side so multiply by 2 to get the total  
         # amount of padding per dimension  
         if pad_l>0 and pad_l<1:  
             for i in xrange(DIM-1):  
                 frac[i]=2*pad_l  
         elif pad_l>=1:  
             for i in xrange(DIM-1):  
                 frac[i]=2*pad_l/float(NE[i])  
         if pad_h>0 and pad_h<1:  
             frac[DIM-1]=2*pad_h  
         elif pad_h>=1:  
             frac[DIM-1]=2*pad_h/(float(NE[DIM-1]))  
         # calculate new number of elements  
         NE_new=[int(NE[i]*(1+frac[i])) for i in xrange(DIM)]  
         NEdiff=[NE_new[i]-NE[i] for i in xrange(DIM)]  
         spacing=[l[i]/NE[i] for i in xrange(DIM)]  
         l_new=[NE_new[i]*spacing[i] for i in xrange(DIM)]  
         origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in xrange(DIM)]  
         return NE_new, l_new, origin_new  
   
     def _interpolateOnDomain(self, data, shape, origin, spacing, length):  
         """  
         Helper method that interpolates data arrays onto the domain.  
         """  
         dim=len(shape)  
         arrays=np.zeros(((len(data[0])-dim),)+tuple(shape))  
         for entry in data:  
             index=()  
             for i in range(dim):  
                 index=((entry[i]-origin[i])/spacing[i],)+index  
             for i in range(arrays.shape[0]):  
                 arrays[i][index]=entry[dim+i]  
         dom=self.getDomain()  
         x=dom.getX()  
         delta=[length[i]/(shape[dim-i-1]-1) for i in xrange(dim)]  
         realorigin=[inf(x[i]) for i in xrange(dim)]  
         res=[]  
         for i in range(arrays.shape[0]):  
             res.append(interpolateTable(arrays[i], x[:dim], realorigin, delta, 1e9))  
         return res  
   
     def setPadding(self, pad_l=0.1, pad_h=0.1):  
         """  
         Sets the amount of padding around the dataset. If pad_l/pad_h is >=1  
         they are treated as number of elements to be added to the domain.  
         If 0 < pad_l;pad_h < 1, the padding amount is relative.  
         """  
         self._pad_l=pad_l  
         self._pad_h=pad_h  
   
     def getDomain(self):  
         """  
         Returns a domain that spans the data area plus padding.  
         """  
         if self._domain is None:  
             self._domain=self._createDomain(self._pad_l, self._pad_h)  
         return self._domain  
118    
119      def getDensityMask(self):      def getDataType(self):
120          """          """
121          Returns the density mask data object, where mask has value 1 on the          Returns the type of survey data managed by this source.
122          padding area, 0 elsewhere.          Subclasses must return `GRAVITY` or `MAGNETIC` as appropriate.
123          """          """
124          raise NotImplementedError          raise NotImplementedError
125    
126      def getGravityAndStdDev(self):      def getSurveyData(self, domain, origin, NE, spacing):
127          """          """
128          Returns the gravity anomaly and standard deviation data objects as a          This method is called by the `DomainBuilder` to retrieve the survey
129          tuple.          data as `Data` objects on the given domain.
130            Subclasses should return one or more `Data` objects with survey data
131            interpolated on the given ripley domain. The exact return type
132            depends on the type of data.
133    
134            :param domain: the escript domain to use
135            :type domain: `esys.escript.Domain`
136            :param origin: the origin coordinates of the domain
137            :type origin: ``tuple`` or ``list``
138            :param NE: the number of domain elements in each dimension
139            :type NE: ``tuple`` or ``list``
140            :param spacing: the cell sizes (node spacing) in the domain
141            :type spacing: ``tuple`` or ``list``
142          """          """
143          raise NotImplementedError          raise NotImplementedError
144    
145      def _createDomain(self, padding_l, padding_h):      def setSubsamplingFactor(self, f):
146          """          """
147          creates and returns an escript domain that spans the entire area of          Sets the data subsampling factor (default=1).
148          available data plus a buffer zone.          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          raise NotImplementedError          self.__subsampling_factor=f
155    
156        def getSubsamplingFactor(self):
157            """
158            Returns the subsampling factor that was set via `setSubsamplingFactor`
159            (see there).
160            """
161            return self.__subsampling_factor
162    
163    
164  ##############################################################################  ##############################################################################
165  class UBCDataSource(DataSource):  class ErMapperData(DataSource):
166      def __init__(self, domainclass, meshfile, gravfile, topofile=None):      """
167          super(UBCDataSource,self).__init__()      Data Source for ER Mapper raster data.
168          self._meshfile=meshfile      Note that this class only accepts a very specific type of ER Mapper data
169          self._gravfile=gravfile      input and will raise an exception if other data is found.
170          self._topofile=topofile      """
171          self._domainclass=domainclass      def __init__(self, datatype, headerfile, datafile=None, altitude=0.):
172          self._readMesh()          """
173            :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`
174      def getDensityMask(self):          :type datatype: ``int``
175          #topodata=self._readTopography()          :param headerfile: ER Mapper header file (usually ends in .ers)
176          #shape=[self.NE[1]+1, self.NE[0]+1]          :type headerfile: ``str``
177          #mask=self._interpolateOnDomain(topodata, shape, self._origin, self._spacing, self._meshlen)          :param datafile: ER Mapper binary data file name. If not supplied the
178          #mask=wherePositive(self.getDomain().getX()[2]-mask[0])                           name of the header file without '.ers' is assumed
179          return self._mask          :type datafile: ``str``
180            :param altitude: altitude of measurements above ground in meters
181      def getGravityAndStdDev(self):          :type altitude: ``float``
182          gravlist=self._readGravity() # x,y,z,g,s          """
183          shape=[self.NE[2]+1, self.NE[1]+1, self.NE[0]+1]          super(ErMapperData,self).__init__()
184          g_and_sigma=self._interpolateOnDomain(gravlist, shape, self._origin, self._spacing, self._meshlen)          self.__headerfile=headerfile
185          return g_and_sigma[0]*[0,0,1], g_and_sigma[1]          if datafile is None:
186                self.__datafile=headerfile[:-4]
187      def _readMesh(self):          else:
188          meshdata=open(self._meshfile).readlines()              self.__datafile=datafile
189          numDataPoints=meshdata[0].split()          self.__altitude=altitude
190          origin=meshdata[1].split()          self.__datatype=datatype
191          self._numDataPoints=[int(x) for x in numDataPoints]          self.__readHeader()
192          self._origin=[float(x) for x in origin]  
193          self._spacing=[float(X.split('*')[1]) for X in meshdata[2:]]      def __readHeader(self):
194          # vertical data is upside down          self.logger.debug("Checking Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))
195          self._origin[2]-=(self._numDataPoints[2]-1)*self._spacing[2]          metadata=open(self.__headerfile, 'r').readlines()
196            # parse metadata
197      def _createDomain(self, padding_l, padding_h):          start=-1
198          NE=[self._numDataPoints[i]-1 for i in xrange(3)]          for i in range(len(metadata)):
199          l=[NE[i]*self._spacing[i] for i in xrange(3)]              if metadata[i].strip() == 'DatasetHeader Begin':
200          NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \                  start=i+1
201                  NE, l, self._origin)          if start==-1:
202                raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
203          self._meshlen=l_new  
204          self.NE=NE_new          md_dict={}
205          self._origin=origin_new          section=[]
206          lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]          for i in range(start, len(metadata)):
207          NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]              line=metadata[i]
208                if line[-6:].strip() == 'Begin':
209                    section.append(line[:-6].strip())
210                elif line[-4:].strip() == 'End':
211                    if len(section)>0:
212                        section.pop()
213                else:
214                    vals=line.split('=')
215                    if len(vals)==2:
216                        key = vals[0].strip()
217                        value = vals[1].strip()
218                        fullkey='.'.join(section+[key])
219                        md_dict[fullkey]=value
220    
221          try:          try:
222              dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])              if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
223              x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]                  raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
224              mask=wherePositive(dom.getX()[2])          except KeyError:
225                self.logger.warn("Cell type not specified. Assuming IEEE4ByteReal.")
         except TypeError:  
             dom=self._domainclass(*self.NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])  
             x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(x[2]+self._origin[2])  
   
         M=2 # do not constrain bottom  
         #M=3 # constrain bottom  
         for i in xrange(M):  
             mask=mask + whereNegative(x[i]) + \  
                     wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])  
         self._mask=wherePositive(mask)  
   
         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]))  
   
         return dom  
   
     def _readTopography(self):  
         f=open(self._topofile)  
         n=int(f.readline())  
         topodata=np.zeros((n,3))  
         for i in xrange(n):  
             x=f.readline().split()  
             x[0]=float(x[0])  
             x[1]=float(x[1])  
             x[2]=float(x[2])  
             topodata[i]=x  
         f.close()  
         return topodata  
226    
227      def _readGravity(self):          try:
228          f=open(self._gravfile)              NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
229          n=int(f.readline())              NY = int(md_dict['RasterInfo.NrOfLines'])
230          gravdata=np.zeros((n,5))          except:
231          for i in xrange(n):              raise RuntimeError("Could not determine extents of data")
             x=f.readline().split()  
             x[0]=float(x[0]) # x  
             x[1]=float(x[1]) # y  
             x[2]=float(x[2]) # z  
             x[3]=float(x[3]) # gravity anomaly in mGal  
             x[4]=float(x[4]) # 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  
232    
233  ##############################################################################          try:
234  class NetCDFDataSource(DataSource):              maskval = float(md_dict['RasterInfo.NullCellValue'])
235      def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):          except:
236                maskval = 99999
237    
238            try:
239                spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
240                spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
241            except:
242                raise RuntimeError("Could not determine cell dimensions")
243    
244            try:
245                if md_dict['CoordinateSpace.CoordinateType']=='EN':
246                    originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
247                    originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
248                elif md_dict['CoordinateSpace.CoordinateType']=='LL':
249                    originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
250                    originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
251                else:
252                    raise RuntimeError("Unknown CoordinateType")
253            except:
254                self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
255                originX,originY = 0.0, 0.0
256    
257            if 'GEODETIC' in md_dict['CoordinateSpace.Projection']:
258                # it appears we have lat/lon coordinates so need to convert
259                # origin and spacing. Try using gdal to get the wkt if available:
260                try:
261                    from osgeo import gdal
262                    ds=gdal.Open(self.__headerfile)
263                    wkt=ds.GetProjection()
264                except:
265                    wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
266                    self.logger.warn('GDAL not available or file read error, assuming GDA94 data')
267                originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)
268                op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
269                # we are rounding to avoid interpolation issues
270                spacingX=np.round(op1X-originX_UTM)
271                spacingY=np.round(op1Y-originY_UTM)
272                originX=np.round(originX_UTM)
273                originY=np.round(originY_UTM)
274    
275            # data sets have origin in top-left corner so y runs top-down
276            self.__dataorigin=[originX, originY]
277            originY-=(NY-1)*spacingY
278            self.__delta = [spacingX, spacingY]
279            self.__maskval = maskval
280            self.__nPts = [NX, NY]
281            self.__origin = [originX, originY]
282            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):
290          """          """
291          vertical_extents - (alt_min, alt_max, num_elements)          returns ( (x0, y0), (nx, ny), (dx, dy) )
         alt_of_data - altitude of measurements  
292          """          """
293          super(NetCDFDataSource,self).__init__()          return (list(self.__origin), list(self.__nPts), list(self.__delta))
         self._topofile=topofile  
         self._gravfile=gravfile  
         self._domainclass=domainclass  
         self._determineExtents(vertical_extents)  
         self._altOfData=alt_of_data  
   
     def getDensityMask(self):  
         #topodata=self._readTopography()  
         #shape=self._numDataPoints[1::-1]  
         #mask=self._interpolateOnDomain(topodata, shape, self._origin, self._spacing, self._meshlen)  
         #mask=wherePositive(self.getDomain().getX()[2]-mask[0])  
         #rho=fill*(1.-mask) + RHO_AIR*mask  
         return self._mask  
   
     def getGravityAndStdDev(self):  
         gravlist=self._readGravity() # x,y,z,g,s  
         shape=[self.NE[2]+1, self.NE[1]+1, self.NE[0]+1]  
         g_and_sigma=self._interpolateOnDomain(gravlist, shape, self._origin, self._spacing, self._meshlen)  
         return g_and_sigma[0]*[0,0,1], g_and_sigma[1]  
294    
295      def _determineExtents(self, ve):      def getDataType(self):
296          f=netcdf_file(self._gravfile, 'r')          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          NX=0
341          for n in ['lon','longitude','x']:          for n in ['lon','longitude','x']:
342              if n in f.dimensions:              if n in f.dimensions:
# Line 278  class NetCDFDataSource(DataSource): Line 353  class NetCDFDataSource(DataSource):
353              raise RuntimeError("Could not determine extents of data")              raise RuntimeError("Could not determine extents of data")
354    
355          # find longitude and latitude variables          # find longitude and latitude variables
         variables=f.variables  
356          lon_name=None          lon_name=None
357          for n in ['lon','longitude']:          for n in ['lon','longitude']:
358              if n in variables:              if n in f.variables:
359                  lon_name=n                  lon_name=n
360                  variables.pop(n)                  longitude=f.variables.pop(n)
361                  break                  break
362          if lon_name is None:          if lon_name is None:
363              raise RuntimeError("Could not determine longitude variable")              raise RuntimeError("Could not determine longitude variable")
364          lat_name=None          lat_name=None
365          for n in ['lat','latitude']:          for n in ['lat','latitude']:
366              if n in variables:              if n in f.variables:
367                  lat_name=n                  lat_name=n
368                  variables.pop(n)                  latitude=f.variables.pop(n)
369                  break                  break
370          if lat_name is None:          if lat_name is None:
371              raise RuntimeError("Could not determine latitude variable")              raise RuntimeError("Could not determine latitude variable")
372    
373          # try to figure out gravity variable name          # try to figure out data variable name
374          grav_name=None          data_name=None
375          if len(variables)==1:          if len(f.variables)==1:
376              grav_name=variables.keys()[0]              data_name=f.variables.keys()[0]
377          else:          else:
378              for n in variables.keys():              for n in f.variables.keys():
379                  dims=variables[n].dimensions                  dims=f.variables[n].dimensions
380                  if (lat_name in dims) and (lon_name in dims):                  if (lat_name in dims) and (lon_name in dims):
381                      grav_name=n                      data_name=n
382                      break                      break
383          if grav_name is None:          if data_name is None:
384              raise RuntimeError("Could not determine gravity variable")              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          origin=[0.,0.,ve[0]]          # see if there is a wkt string to convert coordinates
         lengths=[1.,1.,ve[1]-ve[0]]  
406          try:          try:
407              lon_range=f.variables[lon_name].actual_range              wkt_string=f.variables[data_name].esri_pe_string
             lat_range=f.variables[lat_name].actual_range  
             #origin[:2]=[lon_range[0],lat_range[0]]  
             lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]  
408          except:          except:
409              try:              wkt_string=None
                 #origin[:2]=[f.geospatial_lon_min, f.geospatial_lat_min]  
                 lengths[:2]=[f.geospatial_lon_max-origin[0], f.geospatial_lat_max-origin[1]]  
             except:  
                 pass  
410    
411          # hack          # we don't trust actual_range & geospatial_lon_min/max since subset
412          # grid spacing ~800m          # data does not seem to have these fields updated.
413          SPACING_X=SPACING_Y=800.          # Getting min/max from the arrays is obviously not very efficient but..
414          lengths[0]=(NX-1)*SPACING_X          #lon_range=longitude.actual_range
415          lengths[1]=(NY-1)*SPACING_Y          #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()          f.close()
         self._numDataPoints=[NX, NY, ve[2]]  
         self._origin=origin  
         self._meshlen=lengths  
         self._spacing=[lengths[i]/(self._numDataPoints[i]-1) for i in xrange(3)]  
         self._lon=lon_name  
         self._lat=lat_name  
         self._grv=grav_name  
   
     def _createDomain(self, padding_l, padding_h):  
         NE=[self._numDataPoints[i]-1 for i in xrange(3)]  
         l=self._meshlen  
         NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \  
                 NE, l, self._origin)  
   
         self._meshlen=l_new  
         self.NE=NE_new  
         self._origin=origin_new  
         lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]  
         NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]  
         try:  
             dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])  
             # ripley may adjust NE and length  
             self._meshlen=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in xrange(3)]  
             self.NE=[self._meshlen[i]/(self._spacing[i]) for i in xrange(3)]  
             x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(dom.getX()[2])  
   
         except TypeError:  
             dom=self._domainclass(*self.NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])  
             x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(x[2]+self._origin[2])  
   
         M=2 # do not constrain bottom  
         #M=3 # constrain bottom  
         for i in xrange(M):  
             mask=mask + whereNegative(x[i]) + \  
                     wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])  
         self._mask=wherePositive(mask)  
   
         self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))  
         self.logger.debug("     length: %s x %s x %s"%(self._meshlen[0],self._meshlen[1],self._meshlen[2]))  
         self.logger.debug("     origin: %s x %s x %s"%(self._origin[0],self._origin[1],self._origin[2]))  
   
         return dom  
   
     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")  
424    
425          topodata=np.column_stack((lon,lat,alt))          self.__nPts=[NX, NY]
426          f.close()          self.__origin=[lon_range[0],lat_range[0]]
427          return topodata          # 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 getDataExtents(self):
436            """
437            returns ( (x0, y0), (nx, ny), (dx, dy) )
438            """
439            return (list(self.__origin), list(self.__nPts), list(self.__delta))
440    
441        def getDataType(self):
442            return self.__datatype
443    
444        def getSurveyData(self, domain, origin, NE, spacing):
445            nValues=self.__nPts
446            # determine base location of this dataset within the domain
447            first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
448            if domain.getDim()==3:
449                first.append(int((self.__altitude-origin[2])/spacing[2]))
450                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    
     def _readGravity(self):  
         f=netcdf_file(self._gravfile, 'r')  
         #lon=f.variables[self._lon][:]  
         #lat=f.variables[self._lat][:]  
         grav=f.variables[self._grv][:]  
         f.close()  
         lon=np.linspace(0, (grav.shape[1]-1)*self._spacing[0], num=grav.shape[1])  
         lat=np.linspace(0, (grav.shape[0]-1)*self._spacing[1], num=grav.shape[0])  
         lon,lat=np.meshgrid(lon,lat)  
         lon=lon.flatten()  
         lat=lat.flatten()  
         grav=grav.flatten()  
         alt=self._altOfData*np.ones(grav.shape)  
         try:  
             missing=grav.missing_value  
             err=np.where(grav==missing, 1.0, 0.0)  
         except:  
             err=np.ones(lon.shape)  
         err=1e-6*err  
         grav=1e-6*grav  
         gravdata=np.column_stack((lon,lat,alt,grav,err))  
         return gravdata  
459    
460  class SmoothAnomaly(object):  ##############################################################################
461      def __init__(self, lx, ly, lz, x, y, depth, rho_inner, rho_outer):  class SourceFeature(object):
462        """
463        A feature adds a density distribution to (parts of) a domain of a synthetic
464        data source, for example a layer of a specific rock type or a simulated
465        ore body.
466        """
467        def getValue(self):
468            """
469            Returns the value for the area covered by mask. It can be constant
470            or a data object with spatial dependency.
471            """
472            raise NotImplementedError
473        def getMask(self, x):
474            """
475            Returns the mask of the area of interest for this feature. That is,
476            mask is non-zero where the density returned by getDensity() should be
477            applied, zero elsewhere.
478            """
479            raise NotImplementedError
480    
481    class SmoothAnomaly(SourceFeature):
482        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.rho=None          self.value=None
492            self.mask=None
493    
494        def getValue(self,x):
495            if self.value is None:
496                if self.v_outer is None or self.v_inner is None:
497                    self.value=0
498                else:
499                    DIM=x.getDomain().getDim()
500                    alpha=-log(abs(self.v_outer/self.v_inner))*4
501                    value=exp(-alpha*((x[0]-self.x)/self.lx)**2)
502                    value=value*exp(-alpha*((x[DIM-1]+self.depth)/self.lz)**2)
503                    self.value=maximum(abs(self.v_outer), abs(self.v_inner*value))
504                    if self.v_inner<0: self.value=-self.value
505    
506            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))
514          if self.rho is None:          self.mask = m
515              alpha=-log(abs(self.rho_outer/self.rho_inner))*4          return m
516              rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)  
517              rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)  ##############################################################################
518              self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))  class SyntheticDataBase(DataSource):
519              if self.rho_inner<0: self.rho=-self.rho    """
520          return m        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  class SyntheticDataSource(DataSource):    grid with ``number_of_elements`` cells in each direction.
523      def __init__(self, DIM, NE, l, h, features):    
524          super(SyntheticDataSource,self).__init__()    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
695          self.DIM=DIM  
696          self.NE=NE  
697          self.l=l      def getReferenceProperty(self, domain=None):
698          self.h=h          """
699        Returns the reference density Data object that was used to generate
700      def _createDomain(self, padding_l, padding_h):      the gravity/susceptibility anomaly data.
701          NE_H=self.NE          """
702          NE_L=int((self.l/self.h)*NE_H+0.5)          if self._reference_data == None:
703          l=[self.l]*(self.DIM-1)+[self.h]              DIM=domain.getDim()
704          NE=[NE_L]*(self.DIM-1)+[NE_H]              x=domain.getX()
705          origin=[0.]*self.DIM              k=0.
706          NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \              for f in self._features:
707                  NE, l, origin)                  m=f.getMask(x)
708                    k = k * (1-m) + f.getValue(x) * m
709          self.NE=NE_new              self._reference_data= k
710          self.l=l_new[0]          return self._reference_data
711          self.h=l_new[self.DIM-1]          
712    class SyntheticData(SyntheticDataBase):
713          if self.DIM==2:      """
714              from esys.finley import Rectangle      defines synthetic  gravity/magnetic data based on harmonic property anomaly
715              dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])      
716              self._x = dom.getX() + origin_new          rho = mean + amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * x - shift )
717              self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))          
718              self.logger.debug("     length: %g x %g"%(l_new[0],l_new[1]))      for all x and z<=0. for z>0 rho = 0.        
719              self.logger.debug("     origin: %g x %g"%(origin_new[0],origin_new[1]))      """
720          else:        def __init__(self, datatype,
721              from esys.finley import Brick                         n_length=1,
722              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])                         n_depth=1,
723              self._x = dom.getX() + origin_new                         depth_offset=0.,
724              self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))                         depth=None,
725              self.logger.debug("     length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))                         amplitude=None,
726              self.logger.debug("     origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))                         DIM=2,
727                           number_of_elements=10,
728          dz=l_new[self.DIM-1]/NE_new[self.DIM-1]                         length=1*U.km,
729          self._g_mask=wherePositive(dom.getX()[0]-origin_new[0]) \                         B_b=None,
730                  * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \                         data_offset=0,
731                  * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \                         full_knowledge=False,
732                  * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))                         spherical=False):
733          self._mask=whereNegative(self._x[self.DIM-1]) + \          """
734                  wherePositive(self._x[self.DIM-1]-l[self.DIM-1])          :param datatype: data type indicator
735          for i in xrange(self.DIM-1):          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
736              self._mask=self._mask + whereNegative(self._x[i]) + \          :param n_length: number of oscillation in the anomaly data within the observation region.
737                      wherePositive(self._x[i]-l[i])          :type n_length: ``int``
738          self._mask=wherePositive(self._mask)          :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          rho_ref=0.          :type depth: ``float``
741          for f in self._features:          :param amplitude: data amplitude. Default value is 200 *U.kg/U.m**3 for gravity and 0.1 for magnetic data.
742              m=f.getMask(self._x)          :param features: list of features. It is recommended that the features do entirely lay below surface.
743              rho_ref = rho_ref * (1-m) + f.rho * m          :type features: ``list`` of ``SourceFeature``
744          self._rho=rho_ref          :param DIM: spatial dimension
745            :type DIM: 2 or 3
746          return dom          :param number_of_elements: lateral number of elements in the region where data are collected
747            :type number_of_elements: ``int``
748      def getDensityMask(self):          :param length: lateral extend of the region where data are collected
749          return self._mask          :type length: ``float``
750            :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
751      def getReferenceDensity(self):          :type B_b: ``list`` of ``Scalar``
752          return self._rho          :param data_offset: offset of the data collection region from the surface
753            :type data_offset: ``float``
754      def getGravityAndStdDev(self):          :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
755          pde=LinearSinglePDE(self.getDomain())          :type full_knowledge: ``Bool``
756          G=6.6742e-11*U.m**3/(U.kg*U.sec**2)          :param spherical: if ``True`` sperical coordinates are used (ignored)
757          m_psi_ref=0.          :type spherical: ``Bool``
758          for i in range(self.DIM):          """      
759              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \          super(SyntheticData,self).__init__(
760                      + whereZero(self._x[i]-sup(self._x[i]))                                   datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,
761                                     length=length, B_b=B_b,
762          pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)                                   data_offset=data_offset,
763          pde.setSymmetryOn()                                   full_knowledge=full_knowledge,
764          psi_ref=pde.getSolution()                                   spherical=spherical)
765          del pde          self.__n_length = n_length
766          g=-grad(psi_ref)          self.__n_depth = n_depth
767          sigma=self._g_mask          self.depth = depth
768          return g,sigma          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.3947  
changed lines
  Added in v.4121

  ViewVC Help
Powered by ViewVC 1.1.26