/[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 3965 by caltinay, Fri Sep 14 01:23:17 2012 UTC revision 4108 by caltinay, Thu Dec 13 06:38:11 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','ERSDataSource','SyntheticDataSource','SmoothAnomaly']  __all__ = ['simpleBackgroundMagneticField', 'DataSource','ErMapperData','SyntheticFeatureData','SmoothAnomaly']
26    
27  import logging  import logging
28  import numpy as np  import numpy as np
29  import struct  from esys.escript import ReducedFunction
30  from esys.escript import *  from esys.escript.linearPDEs import LinearSinglePDE
31  from esys.escript.linearPDEs import *  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    
# Line 37  def LatLonToUTM(lon, lat, wkt_string=Non Line 42  def LatLonToUTM(lon, lat, wkt_string=Non
42      """      """
43      Converts one or more longitude,latitude pairs to the corresponding x,y      Converts one or more longitude,latitude pairs to the corresponding x,y
44      coordinates in the Universal Transverse Mercator projection.      coordinates in the Universal Transverse Mercator projection.
45      If wkt_string is not given or invalid or the gdal module is not available  
46      to convert the string, then the input values are assumed to be given in the      :note: If the ``pyproj`` module is not installed a warning is printed and
47      Clarke 1866 projection.             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      # not really optimal: if pyproj is not installed we return the input
60      # values without modification.      # values scaled by a constant.
61      try:      try:
62          import pyproj          import pyproj
63      except:      except:
64          print("Warning, pyproj not available. Domain extents will be wrong")          print("Warning, pyproj not available. Domain extents will be wrong")
65          return lon,lat          return np.array(lon)*1000., np.array(lat)*1000.
66    
67      # determine UTM zone from the input data      # determine UTM zone from the input data
68      zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))      zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
# Line 64  def LatLonToUTM(lon, lat, wkt_string=Non Line 78  def LatLonToUTM(lon, lat, wkt_string=Non
78      x,y=pyproj.transform(p_src, p_dest, lon, lat)      x,y=pyproj.transform(p_src, p_dest, lon, lat)
79      return x,y      return x,y
80    
81    def simpleBackgroundMagneticField(latitude, longitude=0.):
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._constrainBottom=False  
         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):
         """  
         Helper method that computes new number of elements, length and origin  
         after adding padding to the input values.  
         """  
         DIM=len(NE)  
         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):  
         """  
         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., Function(dom)).getNumberOfDataPoints()  
         for i in xrange(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 xrange(num_arrays):  
             d=Scalar(0., Function(dom))  
             d.expand()  
             arrays.append(d)  
   
         for entry in data:  
             index=[int((entry[i]-self._dom_origin[i])/self._spacing[i]) for i in xrange(dim)]  
             index=int(idx_mult.dot(index))  
             for i in xrange(num_arrays):  
                 for p in xrange(DPP):  
                     arrays[i].setValueOfDataPoint(index+p, entry[dim+i])  
   
         return arrays  
   
     def _interpolateOnDomain_old(self, data):  
         """  
         Old interpolation method. Works only on ContinuousFunction and thus  
         produces wrong values once interpolated on Function.  
         """  
         dom=self.getDomain()  
         dim=dom.getDim()  
         # shape = number of data points/nodes in each dimension  
         shape=()  
         for i in xrange(dim):  
             shape=(self._dom_NE[i]+1,)+shape  
         # separate data arrays and coordinates  
         arrays=np.zeros(((len(data[0])-dim),)+shape)  
         num_arrays=arrays.shape[0]  
         for entry in data:  
             index=()  
             for i in xrange(dim):  
                 index=(int((entry[i]-self._dom_origin[i])/self._spacing[i]),)+index  
             for i in xrange(num_arrays):  
                 arrays[i][index]=entry[dim+i]  
         x=dom.getX()  
         delta=[self._dom_len[i]/(shape[dim-i-1]-1) for i in xrange(dim)]  
         realorigin=[inf(x[i]) for i in xrange(dim)]  
         res=[]  
         for i in xrange(num_arrays):  
             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 setConstrainBottom(self, constrain):  
         """  
         If `constrain` 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.  
         """  
         self._constrainBottom=constrain  
   
     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  
   
     def getDensityMask(self):  
         """  
         Returns the density mask data object, where mask has value 1 on the  
         padding area, 0 elsewhere.  
108          """          """
109          return self._mask          returns a tuple of tuples ``( (x0, y0), (nx, ny), (dx, dy) )``, where
110    
111      def getGravityAndStdDev(self):          - ``x0``, ``y0`` = coordinates of data origin
112          """          - ``nx``, ``ny`` = number of data points in x and y
113          Returns the gravity anomaly and standard deviation data objects as a          - ``dx``, ``dy`` = spacing of data points in x and y
         tuple.  
         """  
         raise NotImplementedError  
114    
115      def getDataExtents(self):          This method must be implemented in subclasses.
         """  
         returns ( (x0, y0), (nx, ny), (dx, dy) ), where  
             x0, y0 = coordinates of data origin  
             nx, ny = number of data points in x and y  
             dx, dy = spacing of data points in x and y  
116          """          """
117          raise NotImplementedError          raise NotImplementedError
118    
119      def getVerticalExtents(self):      def getDataType(self):
120          """          """
121          returns (z0, nz, dz), where          Returns the type of survey data managed by this source.
122              z0 = minimum z coordinate (origin)          Subclasses must return `GRAVITY` or `MAGNETIC` as appropriate.
             nz = number of nodes in z direction  
             dz = spacing of nodes (= cell size in z)  
123          """          """
124          raise NotImplementedError          raise NotImplementedError
125    
126      def getDomainClass(self):      def getSurveyData(self, domain, origin, NE, spacing):
127          """          """
128          returns the domain generator class (e.g. esys.ripley.Brick)          This method is called by the `DomainBuilder` to retrieve the survey
129            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          X0, NX, DX = self.getDataExtents()          self.__subsampling_factor=f
         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 = [np.round(si) for si in self._spacing]  
   
         # length of domain (without padding)  
         l = [NE[i]*self._spacing[i] for i in xrange(len(NE))]  
   
         # now add padding to the values  
         NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \  
                 NE, l, origin)  
155    
156          # number of padding elements per side      def getSubsamplingFactor(self):
         NE_pad=[(NE_new[i]-NE[i])/2 for i in xrange(3)]  
   
         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 xrange(3)]  
         try:  
             dom=self.getDomainClass()(*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 xrange(3)]  
             self._dom_NE=[self._dom_len[i]/self._spacing[i] for i in xrange(3)]  
             x=dom.getX()-[self._dom_origin[i]+NE_pad[i]*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(dom.getX()[2])  
   
         except TypeError:  
             dom=self.getDomainClass()(*self._dom_NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])  
             x=dom.getX()-[NE_pad[i]*self._spacing[i] for i in xrange(3)]  
             mask=wherePositive(x[2]+self._dom_origin[2])  
   
         # prepare density mask (=1 at padding area, 0 else)  
         if self._constrainBottom:  
             M=3 # constrain bottom  
         else:  
             M=2 # do not constrain bottom  
         for i in xrange(M):  
             mask=mask + whereNegative(x[i]) + \  
                     wherePositive(x[i]-l_new[i]+2*NE_pad[i]*self._spacing[i])  
         self._mask=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, domainclass, meshfile, gravfile, topofile=None):  
         super(UBCDataSource,self).__init__()  
         self.__meshfile=meshfile  
         self.__gravfile=gravfile  
         self.__topofile=topofile  
         self.__domainclass=domainclass  
         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)  
157          """          """
158          return (self.__origin[2], self.__nPts[2], self.__delta[2])          Returns the subsampling factor that was set via `setSubsamplingFactor`
159            (see there).
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
160          """          """
161          return self.__domainclass          return self.__subsampling_factor
   
     #def getDensityMask(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]  
162    
     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=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 xrange(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  
163    
164  ##############################################################################  ##############################################################################
165  class NetCDFDataSource(DataSource):  class ErMapperData(DataSource):
     def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):  
         """  
         vertical_extents - (alt_min, alt_max, num_points)  
         alt_of_data - altitude of measurements  
         """  
         super(NetCDFDataSource,self).__init__()  
         self.__topofile=topofile  
         self.__gravfile=gravfile  
         self.__domainclass=domainclass  
         self.__determineExtents(vertical_extents)  
         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")  
   
         # 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 it seems.  
         # Getting min/max from the arrays is obviously not very efficient...  
         #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]  
         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 interpolateOnDomain issues  
         self.__delta=[np.round(lengths[i]/(self.__nPts[i]-1)) for i in xrange(3)]  
         self.__wkt_string=wkt_string  
         self.__lon=lon_name  
         self.__lat=lat_name  
         self.__grv=grav_name  
   
     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)  
         """  
         return (self.__origin[2], self.__nPts[2], self.__delta[2])  
   
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
         """  
         return self.__domainclass  
   
     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=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")  
   
         topodata=np.column_stack((lon,lat,alt))  
         f.close()  
         return topodata  
   
     def _readGravity(self):  
         f=netcdf_file(self.__gravfile, 'r')  
         #lon=f.variables[self.__lon][:]  
         #lat=f.variables[self.__lat][:]  
         NE=[self.__nPts[i]-1 for i in xrange(2)]  
         lon=np.linspace(self.__origin[0], self.__origin[0]+NE[0]*self.__delta[0], NE[0]+1)  
         lat=np.linspace(self.__origin[1], self.__origin[1]+NE[1]*self.__delta[1], NE[1]+1)  
         lon,lat=np.meshgrid(lon,lat)  
         grav=f.variables[self.__grv][:]  
         f.close()  
         lon=lon.flatten()  
         lat=lat.flatten()  
         grav=grav.flatten()  
         alt=self.__altOfData*np.ones(grav.shape)  
         # error value is an assumption  
         try:  
             missing=f.variables[self.__grv].missing_value  
             err=np.where(grav==missing, 0.0, 20.0)  
         except:  
             err=20.0*np.ones(lon.shape)  
         # convert units  
         err=1e-6*err  
         grav=1e-6*grav  
         gravdata=np.column_stack((lon,lat,alt,grav,err))  
         return gravdata  
   
 ##############################################################################  
 class ERSDataSource(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, domainclass, headerfile, datafile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):      def __init__(self, datatype, headerfile, datafile=None, altitude=0.):
172          """          """
173          headerfile - usually ends in .ers          :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`
174          datafile - usually has the same name as the headerfile without '.ers'          :type datatype: ``int``
175            :param headerfile: ER Mapper header file (usually ends in .ers)
176            :type headerfile: ``str``
177            :param datafile: ER Mapper binary data file name. If not supplied the
178                             name of the header file without '.ers' is assumed
179            :type datafile: ``str``
180            :param altitude: altitude of measurements above ground in meters
181            :type altitude: ``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.__domainclass=domainclass          self.__altitude=altitude
190          self.__readHeader(vertical_extents)          self.__datatype=datatype
191          self.__altOfData=alt_of_data          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
198          for i in xrange(len(metadata)):          for i in range(len(metadata)):
199              if metadata[i].strip() == 'DatasetHeader Begin':              if metadata[i].strip() == 'DatasetHeader Begin':
200                  start=i+1                  start=i+1
201          if start==-1:          if start==-1:
# Line 575  class ERSDataSource(DataSource): Line 203  class ERSDataSource(DataSource):
203    
204          md_dict={}          md_dict={}
205          section=[]          section=[]
206          for i in xrange(start, len(metadata)):          for i in range(start, len(metadata)):
207              line=metadata[i]              line=metadata[i]
208              if line[-6:].strip() == 'Begin':              if line[-6:].strip() == 'Begin':
209                  section.append(line[:-6].strip())                  section.append(line[:-6].strip())
# Line 605  class ERSDataSource(DataSource): Line 233  class ERSDataSource(DataSource):
233          try:          try:
234              maskval = float(md_dict['RasterInfo.NullCellValue'])              maskval = float(md_dict['RasterInfo.NullCellValue'])
235          except:          except:
236              raise RuntimeError("Could not determine NullCellValue")              maskval = 99999
237    
238          try:          try:
239              spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])              spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
# Line 635  class ERSDataSource(DataSource): Line 263  class ERSDataSource(DataSource):
263                  wkt=ds.GetProjection()                  wkt=ds.GetProjection()
264              except:              except:
265                  wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'                  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, assuming GDA94 data')                  self.logger.warn('GDAL not available or file read error, assuming GDA94 data')
267              originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)              originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)
268              op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)              op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
269              # we are rounding to avoid interpolateOnDomain issues              # we are rounding to avoid interpolation issues
270              spacingX=np.round(op1X-originX_UTM)              spacingX=np.round(op1X-originX_UTM)
271              spacingY=np.round(op1Y-originY_UTM)              spacingY=np.round(op1Y-originY_UTM)
272              originX=np.round(originX_UTM)              originX=np.round(originX_UTM)
# Line 647  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]-1))          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 getVerticalExtents(self):      def getDataType(self):
296          """          return self.__datatype
297          returns (z0, nz, dz)  
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          return (self.__origin[2], self.__nPts[2], self.__delta[2])          :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      def getDomainClass(self):          # 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 getDataExtents(self):
436          """          """
437          returns the domain generator class (e.g. esys.ripley.Brick)          returns ( (x0, y0), (nx, ny), (dx, dy) )
438          """          """
439          return self.__domainclass          return (list(self.__origin), list(self.__nPts), list(self.__delta))
440    
441      def getGravityAndStdDev(self):      def getDataType(self):
442          gravlist=self._readGravity() # x,y,z,g,s          return self.__datatype
         g_and_sigma=self._interpolateOnDomain(gravlist)  
         return g_and_sigma[0]*[0,0,1], g_and_sigma[1]  
443    
444      def _readGravity(self):      def getSurveyData(self, domain, origin, NE, spacing):
445          f=open(self.__datafile, 'r')          nValues=self.__nPts
446          n=self.__nPts[0]*self.__nPts[1]          # determine base location of this dataset within the domain
447          grav=[]          first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
448          err=[]          if domain.getDim()==3:
449          for i in xrange(n):              first.append(int((self.__altitude-origin[2])/spacing[2]))
450              v=struct.unpack('f',f.read(4))[0]              nValues=nValues+[1]
451              grav.append(v)  
452              if abs((self.__maskval-v)/v) < 1e-6:          data=ripleycpp._readNcGrid(self.__filename, self.__data_name,
453                  err.append(0.)                    ReducedFunction(domain), first, nValues, (), self.__maskval)
454              else:          sigma=whereNonZero(data-self.__maskval)
455                  err.append(1.)          data=data*self.__scalefactor
456          f.close()          sigma=sigma * 2. * self.__scalefactor
457            return data, sigma
458    
         NE=[self.__nPts[i] for i in xrange(2)]  
         x=self.__dataorigin[0]+np.arange(start=0, stop=NE[0]*self.__delta[0], step=self.__delta[0])  
         # data sets have origin in top-left corner so y runs top-down  
         y=self.__dataorigin[1]-np.arange(start=0, stop=NE[1]*self.__delta[1], step=self.__delta[1])  
         x,y=np.meshgrid(x,y)  
         x,y=x.flatten(),y.flatten()  
         alt=self.__altOfData*np.ones((n,))  
   
         # convert units  
         err=2e-5*np.array(err)  
         grav=1e-5*np.array(grav)  
         gravdata=np.column_stack((x,y,alt,grav,err))  
         return gravdata  
459    
460  ##############################################################################  ##############################################################################
461  class SourceFeature(object):  class SourceFeature(object):
# Line 711  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 726  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, rho_outer):      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      def getDensity(self):  
494          return self.rho      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
             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 m      
516    
517  ##############################################################################  ##############################################################################
518  class SyntheticDataSource(DataSource):  class SyntheticFeatureData(DataSource):
519      def __init__(self, DIM, NE, l, h, features):      def __init__(self, datatype, DIM, NE, l, features, B_b=None):
520          super(SyntheticDataSource,self).__init__()          super(SyntheticFeatureData,self).__init__()
521            if not datatype in [self.GRAVITY,self.MAGNETIC]:
522                raise ValueError("Invalid value for datatype parameter")
523            self.__datatype = datatype
524          self._features = features          self._features = features
525            self.__origin = [0]*(DIM-1)
526            self.__nPts = [NE]*(DIM-1)
527            self.__delta = [float(l)/NE]*(DIM-1)
528            self.__B_b =None
529          self.DIM=DIM          self.DIM=DIM
530          self.NE=NE          self.NE=NE
531          self.l=l          self.l=l
532          self.h=h          # this is for Cartesian (FIXME ?)
533            if datatype  ==  self.MAGNETIC:
534      def _createDomain(self, padding_l, padding_h):              if self.DIM<3:
535          NE_H=self.NE                 self.__B_b =  np.array([-B_b[2],  -B_b[0]])
536          NE_L=int((self.l/self.h)*NE_H+0.5)              else:
537          l=[self.l]*(self.DIM-1)+[self.h]                 self.__B_b = ([-B_b[1],  -B_b[2],  -B_b[0]])
         NE=[NE_L]*(self.DIM-1)+[NE_H]  
         origin=[0.]*self.DIM  
         NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \  
                 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:  
             from esys.finley import Rectangle  
             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:    
             from esys.finley import Brick  
             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._mask=whereNegative(self._x[self.DIM-1]) + \  
                 wherePositive(self._x[self.DIM-1]-l[self.DIM-1])  
         for i in xrange(self.DIM-1):  
             self._mask=self._mask + whereNegative(self._x[i]) + \  
                     wherePositive(self._x[i]-l[i])  
         self._mask=wherePositive(self._mask)  
   
         rho_ref=0.  
         for f in self._features:  
             m=f.getMask(self._x)  
             rho_ref = rho_ref * (1-m) + f.getDensity() * m  
         self._rho=rho_ref  
538    
539          return dom      def getDataExtents(self):
540            return (list(self.__origin), list(self.__nPts), list(self.__delta))
541    
542      def getDensityMask(self):      def getDataType(self):
543          return self._mask          return self.__datatype
544    
545      def getReferenceDensity(self):      def getReferenceDensity(self):
546            """
547            Returns the reference density Data object that was used to generate
548            the gravity anomaly data.
549            """
550          return self._rho          return self._rho
551    
552      def getGravityAndStdDev(self):      def getReferenceSusceptibility(self):
553          pde=LinearSinglePDE(self.getDomain())          """
554          G=6.6742e-11*U.m**3/(U.kg*U.sec**2)          Returns the reference magnetic susceptibility Data objects that was
555          m_psi_ref=0.          used to generate the magnetic field data.
556          for i in xrange(self.DIM):          """
557              m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \          return self._k
                     + whereZero(self._x[i]-sup(self._x[i]))  
558    
559          pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)      def getSurveyData(self, domain, origin, NE, spacing):
560            pde=LinearSinglePDE(domain)
561            G=U.Gravitational_Constant
562            m_psi_ref=0.
563            x=domain.getX()
564            DIM=domain.getDim()
565            m_psi_ref=whereZero(x[DIM-1]-sup(x[DIM-1])) # + whereZero(x[DIM-1]-inf(x[DIM-1]))
566            if self.getDataType()==DataSource.GRAVITY:
567                rho_ref=0.
568                for f in self._features:
569                    m=f.getMask(x)
570                    rho_ref = rho_ref * (1-m) + f.getValue(x) * m
571                self._rho=rho_ref
572                pde.setValue(A=kronecker(domain), Y=-4*np.pi*G*rho_ref, q=m_psi_ref)
573            else:
574                k_ref=0.
575                for f in self._features:
576                    m=f.getMask(x)
577                    k_ref = k_ref * (1-m) + f.getValue(x) * m
578                self._k=k_ref
579                pde.setValue(A=kronecker(domain), X=k_ref*self.__B_b, q=m_psi_ref)
580          pde.setSymmetryOn()          pde.setSymmetryOn()
581          psi_ref=pde.getSolution()          psi_ref=pde.getSolution()
582          del pde          del pde
583          g=-grad(psi_ref)          if self.getDataType()==DataSource.GRAVITY:
584          sigma=self._g_mask              data = -grad(psi_ref, ReducedFunction(domain))
585          return g,sigma          else:
586                data = self._k*self.__B_b-grad(psi_ref, ReducedFunction(domain))
587    
588            sigma=1.
589            # limit mask to non-padding in horizontal area
590            x=ReducedFunction(domain).getX()
591            for i in range(self.DIM-1):
592                sigma=sigma * wherePositive(x[i]) \
593                            * whereNegative(x[i]-(sup(x[i])+inf(x[i])))
594            # limit mask to one cell thickness at z=0
595            sigma = sigma * whereNonNegative(x[self.DIM-1]) \
596                    * whereNonPositive(x[self.DIM-1]-spacing[self.DIM-1])
597            return data,sigma
598    
599    
600    

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

  ViewVC Help
Powered by ViewVC 1.1.26