/[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 4004 by caltinay, Mon Sep 24 01:08:19 2012 UTC revision 4005 by caltinay, Fri Sep 28 06:09:03 2012 UTC
# Line 24  __all__ = ['DataSource','UBCDataSource', Line 24  __all__ = ['DataSource','UBCDataSource',
24    
25  import logging  import logging
26  import numpy as np  import numpy as np
 import struct  
27  from esys.escript import *  from esys.escript import *
28  from esys.escript.linearPDEs import *  from esys.escript.linearPDEs import LinearSinglePDE
29  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
30  from esys.ripley import *  from esys.ripley import Brick, Rectangle, ripleycpp
31    
32  try:  try:
33      from scipy.io.netcdf import netcdf_file      from scipy.io.netcdf import netcdf_file
# Line 40  def LatLonToUTM(lon, lat, wkt_string=Non Line 39  def LatLonToUTM(lon, lat, wkt_string=Non
39      """      """
40      Converts one or more longitude,latitude pairs to the corresponding x,y      Converts one or more longitude,latitude pairs to the corresponding x,y
41      coordinates in the Universal Transverse Mercator projection.      coordinates in the Universal Transverse Mercator projection.
42      If wkt_string is not given or invalid or the gdal module is not available      If `wkt_string` is not given or invalid or the ``gdal`` module is not
43      to convert the string, then the input values are assumed to be given in the      available to convert the string, then the input values are assumed to be
44      Clarke 1866 projection.      given in the Clarke 1866 projection.
45      """      """
46    
47      # not really optimal: if pyproj is not installed we return the input      # not really optimal: if pyproj is not installed we return the input
48      # values without modification.      # values scaled by a constant.
49      try:      try:
50          import pyproj          import pyproj
51      except:      except:
# Line 70  def LatLonToUTM(lon, lat, wkt_string=Non Line 69  def LatLonToUTM(lon, lat, wkt_string=Non
69  class DataSource(object):  class DataSource(object):
70      """      """
71      A class that provides survey data for the inversion process.      A class that provides survey data for the inversion process.
72        This is an abstract base class that implements common functionality.
73        Methods to be overwritten by subclasses are marked as such.
74      """      """
75      # this is currently specific to gravity inversion and should be generalised      # this is currently specific to gravity inversion and should be generalised
76    
77      def __init__(self):      def __init__(self):
78          """          """
79            Constructor. Sets some defaults and initializes logger.
80          """          """
81          self._constrainBottom=False          self._constrainBottom=False
82          self._constrainSides=True          self._constrainSides=True
83          self._domain=None          self._domain=None
84          self._pad_l=0.1          self.setPadding()
         self._pad_h=0.1  
85          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)          self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
86    
87      def _addPadding(self, pad_l, pad_h, NE, l, origin):      def setPadding(self, pad_x=10, pad_y=10):
88            """
89            Sets the amount of padding around the dataset. If `pad_x`/`pad_y`
90            is >=1 the value is treated as number of elements to be added to the
91            domain (per side).
92            If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the
93            dataset size. For example, calling ``setPadding(3, 0.1)`` to a data
94            source with size 10x20 will result in the padded data set size
95            16x24 (10+2*3, 20*(1+2*0.1))
96    
97            :param pad_x: Padding per side in x direction (default: 10 elements)
98            :type pad_x: ``int`` or ``float``
99            :param pad_y: Padding per side in y direction (default: 10 elements).
100                          This value is only used for 3-dimensional datasets
101            :type pad_y: ``int`` or ``float``
102            """
103            self._pad_x=pad_x
104            self._pad_y=pad_y
105    
106        def setConstraints(self, bottom=False, sides=True):
107            """
108            If `bottom` is True, then the density mask will be set to 1 in the
109            padding area at the bottom of the domain. By default this area is
110            unconstrained. Similarly, if `sides` is True (default) then the
111            horizontal padding area is constrained, otherwise not.
112    
113            :param bottom: Whether to constrain the density at the bottom of the
114                           domain
115            :type bottom: ``bool``
116            :param sides: Whether to constrain the density in the padding area
117                          surrounding the data
118            :type sides: ``bool``
119            """
120            self._constrainBottom=bottom
121            self._constrainSides=sides
122    
123        def getDomain(self):
124            """
125            Returns a domain that spans the data area plus padding.
126            The domain is created the first time this method is called, subsequent
127            calls return the same domain so anything that affects the domain
128            (such as padding) needs to be set beforehand.
129    
130            :return: The escript domain for this data source.
131            :rtype: `esys.escript.Domain`
132            """
133            if self._domain is None:
134                self._domain=self._createDomain(self._pad_x, self._pad_y)
135            return self._domain
136    
137        def getDensityMask(self):
138            """
139            Returns the density mask data object, where mask has value 1 in the
140            padding area, 0 elsewhere.
141    
142            :return: The mask for the density.
143            :rtype: `esys.escript.Data`
144            """
145            return self._mask
146    
147        def getGravityAndStdDev(self):
148            """
149            Returns the gravity anomaly and standard deviation data objects as a
150            tuple. This method must be implemented in subclasses.
151            """
152            raise NotImplementedError
153    
154        def getDataExtents(self):
155            """
156            returns a tuple of tuples ``( (x0, y0), (nx, ny), (dx, dy) )``, where
157    
158            - ``x0``, ``y0`` = coordinates of data origin
159            - ``nx``, ``ny`` = number of data points in x and y
160            - ``dx``, ``dy`` = spacing of data points in x and y
161    
162            This method must be implemented in subclasses.
163            """
164            raise NotImplementedError
165    
166        def getVerticalExtents(self):
167            """
168            returns a tuple ``(z0, nz, dz)``, where
169    
170            - ``z0`` = minimum z coordinate (origin)
171            - ``nz`` = number of nodes in z direction
172            - ``dz`` = spacing of nodes (= cell size in z)
173    
174            This method must be implemented in subclasses.
175            """
176            raise NotImplementedError
177    
178        def getDomainClass(self):
179            """
180            returns the domain generator class (e.g. esys.ripley.Brick).
181            Must be implemented in subclasses.
182            """
183            raise NotImplementedError
184    
185        def _addPadding(self, pad_x, pad_y, NE, l, origin):
186          """          """
187          Helper method that computes new number of elements, length and origin          Helper method that computes new number of elements, length and origin
188          after adding padding to the input values.          after adding padding to the input values.
189    
190            :param pad_x: Number of elements or fraction of padding in x direction
191            :type pad_x: ``int`` or ``float``
192            :param pad_y: Number of elements or fraction of padding in y direction
193            :type pad_y: ``int`` or ``float``
194            :param NE: Initial number of elements
195            :type NE: ``tuple`` or ``list``
196            :param l: Initial side lengths
197            :type l: ``tuple`` or ``list``
198            :param origin: Initial origin
199            :type origin: ``tuple`` or ``list``
200            :return: tuple with three elements ``(NE_padded, l_padded, origin_padded)``,
201                     which are lists of the updated input parameters
202          """          """
203          DIM=len(NE)          DIM=len(NE)
204          frac=[0.]*DIM          frac=[0.]*(DIM-1)+[0]
205          # padding is applied to each side so multiply by 2 to get the total          # padding is applied to each side so multiply by 2 to get the total
206          # amount of padding per dimension          # amount of padding per dimension
207          if pad_l>0 and pad_l<1:          if pad_x>0 and pad_y<1:
208              for i in xrange(DIM-1):              frac[0]=2.*pad_x
209                  frac[i]=2.*pad_l          elif pad_x>=1:
210          elif pad_l>=1:              frac[0]=2.*pad_x/float(NE[0])
211              for i in xrange(DIM-1):          if DIM>2:
212                  frac[i]=2.*pad_l/float(NE[i])              if pad_y>0 and pad_y<1:
213          if pad_h>0 and pad_h<1:                  frac[1]=2.*pad_y
214              frac[DIM-1]=2.*pad_h              elif pad_y>=1:
215          elif pad_h>=1:                  frac[1]=2.*pad_y/(float(NE[1]))
             frac[DIM-1]=2.*pad_h/(float(NE[DIM-1]))  
216    
217          # calculate new number of elements          # calculate new number of elements
218          NE_new=[int(NE[i]*(1+frac[i])) for i in xrange(DIM)]          NE_new=[int(NE[i]*(1+frac[i])) for i in xrange(DIM)]
# Line 145  class DataSource(object): Line 256  class DataSource(object):
256    
257          return arrays          return arrays
258    
     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 setConstraints(self, bottom=False, sides=True):  
         """  
         If `bottom` is True, then the density mask will be set to 1 in the  
         padding area at the bottom of the domain. By default this area is  
         unconstrained. Similarly, if `sides` is True (default) then the  
         horizontal padding area is constrained, otherwise not.  
         """  
         self._constrainBottom=bottom  
         self._constrainSides=sides  
   
     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.  
         """  
         return self._mask  
   
     def getGravityAndStdDev(self):  
         """  
         Returns the gravity anomaly and standard deviation data objects as a  
         tuple.  
         """  
         raise NotImplementedError  
   
     def getDataExtents(self):  
         """  
         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  
         """  
         raise NotImplementedError  
   
     def getVerticalExtents(self):  
         """  
         returns (z0, nz, dz), where  
             z0 = minimum z coordinate (origin)  
             nz = number of nodes in z direction  
             dz = spacing of nodes (= cell size in z)  
         """  
         raise NotImplementedError  
   
     def getDomainClass(self):  
         """  
         returns the domain generator class (e.g. esys.ripley.Brick)  
         """  
         raise NotImplementedError  
   
259      def _createDomain(self, padding_l, padding_h):      def _createDomain(self, padding_l, padding_h):
260          """          """
261          creates and returns an escript domain that spans the entire area of          Creates and returns an escript domain that spans the entire area of
262          available data plus a buffer zone.          available data plus a buffer zone. This method is called only once
263            the first time `getDomain()` is invoked and may be overwritten if
264            required.
265    
266            :return: The escript domain for this data source.
267            :rtype: `esys.escript.Domain`
268          """          """
269          X0, NX, DX = self.getDataExtents()          X0, NX, DX = self.getDataExtents()
270          z0, nz, dz = self.getVerticalExtents()          z0, nz, dz = self.getVerticalExtents()

Legend:
Removed from v.4004  
changed lines
  Added in v.4005

  ViewVC Help
Powered by ViewVC 1.1.26