/[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 3966 by caltinay, Mon Sep 17 00:45:39 2012 UTC revision 3971 by caltinay, Wed Sep 19 02:55:35 2012 UTC
# Line 27  import struct Line 27  import struct
27  from esys.escript import *  from esys.escript import *
28  from esys.escript.linearPDEs import *  from esys.escript.linearPDEs import *
29  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
30    from esys.ripley import *
31    
32  try:  try:
33      from scipy.io.netcdf import netcdf_file      from scipy.io.netcdf import netcdf_file
34      __all__ += ['NetCDFDataSource']      __all__ += ['NetCDFDataSource']
# Line 261  class DataSource(object): Line 263  class DataSource(object):
263          # number of padding elements per side          # number of padding elements per side
264          NE_pad=[(NE_new[i]-NE[i])/2 for i in xrange(3)]          NE_pad=[(NE_new[i]-NE[i])/2 for i in xrange(3)]
265    
266            self._dom_NE_pad = NE_pad
267          self._dom_len = l_new          self._dom_len = l_new
268          self._dom_NE = NE_new          self._dom_NE = NE_new
269          self._dom_origin = origin_new          self._dom_origin = origin_new
# Line 547  class ERSDataSource(DataSource): Line 550  class ERSDataSource(DataSource):
550      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
551      input and will raise an exception if other data is found.      input and will raise an exception if other data is found.
552      """      """
553      def __init__(self, domainclass, headerfile, datafile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):      def __init__(self, headerfile, datafile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):
554          """          """
555          headerfile - usually ends in .ers          headerfile - usually ends in .ers
556          datafile - usually has the same name as the headerfile without '.ers'          datafile - usually has the same name as the headerfile without '.ers'
# Line 558  class ERSDataSource(DataSource): Line 561  class ERSDataSource(DataSource):
561              self.__datafile=headerfile[:-4]              self.__datafile=headerfile[:-4]
562          else:          else:
563              self.__datafile=datafile              self.__datafile=datafile
         self.__domainclass=domainclass  
564          self.__readHeader(vertical_extents)          self.__readHeader(vertical_extents)
565          self.__altOfData=alt_of_data          self.__altOfData=alt_of_data
566    
# Line 669  class ERSDataSource(DataSource): Line 671  class ERSDataSource(DataSource):
671          """          """
672          returns the domain generator class (e.g. esys.ripley.Brick)          returns the domain generator class (e.g. esys.ripley.Brick)
673          """          """
674          return self.__domainclass          return Brick
675    
676      def getGravityAndStdDev(self):      def getGravityAndStdDev(self):
677          gravlist=self.__readGravity() # x,y,z,g,s          nValues=self.__nPts[:2]+[1]
678          g_and_sigma=self._interpolateOnDomain(gravlist)          first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
679          return g_and_sigma[0]*[0,0,1], g_and_sigma[1]          g=ripleycpp._readBinaryGrid(self.__datafile, Function(self.getDomain()),
680                    first, nValues, (), -1e8)
681            sigma=whereNonNegative(g+1e7)
682            g=g*1e-6
683            sigma=sigma*2e-6
684            return g*[0,0,1], sigma
685    
     def __readGravity(self):  
         f=open(self.__datafile, 'r')  
         n=self.__nPts[0]*self.__nPts[1]  
         grav=[]  
         err=[]  
         for i in xrange(n):  
             v=struct.unpack('f',f.read(4))[0]  
             grav.append(v)  
             if abs((self.__maskval-v)/v) < 1e-6:  
                 err.append(0.)  
             else:  
                 err.append(1.)  
         f.close()  
   
         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-6*np.array(err)  
         grav=1e-6*np.array(grav)  
         gravdata=np.column_stack((x,y,alt,grav,err))  
         return gravdata  
686    
687  ##############################################################################  ##############################################################################
688  class SourceFeature(object):  class SourceFeature(object):

Legend:
Removed from v.3966  
changed lines
  Added in v.3971

  ViewVC Help
Powered by ViewVC 1.1.26