/[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 3957 by caltinay, Wed Sep 5 23:49:23 2012 UTC revision 3958 by caltinay, Fri Sep 7 02:56:42 2012 UTC
# Line 19  __license__="""Licensed under the Open S Line 19  __license__="""Licensed under the Open S
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  __all__ = ['DataSource','UBCDataSource','SyntheticDataSource','SmoothAnomaly']  __all__ = ['DataSource','UBCDataSource','ERSDataSource','SyntheticDataSource','SmoothAnomaly']
23    
24  import logging  import logging
25  import numpy as np  import numpy as np
26    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
# Line 71  class DataSource(object): Line 72  class DataSource(object):
72      def __init__(self):      def __init__(self):
73          """          """
74          """          """
75            self._constrainBottom=False
76          self._domain=None          self._domain=None
77          self._pad_l=0.1          self._pad_l=0.1
78          self._pad_h=0.1          self._pad_h=0.1
# Line 111  class DataSource(object): Line 113  class DataSource(object):
113          arrays=np.zeros(((len(data[0])-dim),)+tuple(shape))          arrays=np.zeros(((len(data[0])-dim),)+tuple(shape))
114          for entry in data:          for entry in data:
115              index=()              index=()
116              for i in range(dim):              for i in xrange(dim):
117                  index=((entry[i]-origin[i])/spacing[i],)+index                  index=((entry[i]-origin[i])/spacing[i],)+index
118              for i in range(arrays.shape[0]):              for i in xrange(arrays.shape[0]):
119                  arrays[i][index]=entry[dim+i]                  arrays[i][index]=entry[dim+i]
120          dom=self.getDomain()          dom=self.getDomain()
121          x=dom.getX()          x=dom.getX()
122          delta=[length[i]/(shape[dim-i-1]-1) for i in xrange(dim)]          delta=[length[i]/(shape[dim-i-1]-1) for i in xrange(dim)]
123          realorigin=[inf(x[i]) for i in xrange(dim)]          realorigin=[inf(x[i]) for i in xrange(dim)]
124          res=[]          res=[]
125          for i in range(arrays.shape[0]):          for i in xrange(arrays.shape[0]):
126              res.append(interpolateTable(arrays[i], x[:dim], realorigin, delta, 1e9))              res.append(interpolateTable(arrays[i], x[:dim], realorigin, delta, 1e9))
127          return res          return res
128    
# Line 133  class DataSource(object): Line 135  class DataSource(object):
135          self._pad_l=pad_l          self._pad_l=pad_l
136          self._pad_h=pad_h          self._pad_h=pad_h
137    
138        def setConstrainBottom(self, constrain):
139            """
140            If `constrain` is True, then the density mask will be set to 1 in the
141            padding area at the bottom of the domain. By default this area is
142            unconstrained.
143            """
144            self._constrainBottom=constrain
145    
146      def getDomain(self):      def getDomain(self):
147          """          """
148          Returns a domain that spans the data area plus padding.          Returns a domain that spans the data area plus padding.
# Line 217  class UBCDataSource(DataSource): Line 227  class UBCDataSource(DataSource):
227              x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]              x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
228              mask=wherePositive(x[2]+self._origin[2])              mask=wherePositive(x[2]+self._origin[2])
229    
230          M=2 # do not constrain bottom          if self._constrainBottom:
231          #M=3 # constrain bottom              M=3 # constrain bottom
232            else:
233                M=2 # do not constrain bottom
234          for i in xrange(M):          for i in xrange(M):
235              mask=mask + whereNegative(x[i]) + \              mask=mask + whereNegative(x[i]) + \
236                      wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])                      wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])
# Line 263  class UBCDataSource(DataSource): Line 275  class UBCDataSource(DataSource):
275    
276  ##############################################################################  ##############################################################################
277  class NetCDFDataSource(DataSource):  class NetCDFDataSource(DataSource):
278      def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=1.):      def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):
279          """          """
280          vertical_extents - (alt_min, alt_max, num_elements)          vertical_extents - (alt_min, alt_max, num_elements)
281          alt_of_data - altitude of measurements          alt_of_data - altitude of measurements
# Line 343  class NetCDFDataSource(DataSource): Line 355  class NetCDFDataSource(DataSource):
355          except:          except:
356              wkt_string=None              wkt_string=None
357    
358          origin=[0.,0.,ve[0]]          # we don't trust actual_range & geospatial_lon_min/max since subset
359          lengths=[100000.,100000.,ve[1]-ve[0]]          # data does not seem to have these fields updated it seems.
360            # Getting min/max from the arrays is obviously not very efficient...
361          try:          #lon_range=longitude.actual_range
362              lon_range=longitude.actual_range          #lat_range=latitude.actual_range
363              lat_range=latitude.actual_range          #lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
364              lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)          #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
365              origin[:2]=lon_range[0],lat_range[0]          lon_range=longitude.data.min(),longitude.data.max()
366              lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]          lat_range=latitude.data.min(),latitude.data.max()
367          except:          lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
368              try:          origin=[lon_range[0],lat_range[0],ve[0]]
369                  lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0],ve[1]-ve[0]]
                 lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]  
                 lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)  
                 origin[:2]=lon_range[0],lat_range[0]  
                 lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]  
             except:  
                 pass  
370    
371          f.close()          f.close()
372    
373          self._numDataPoints=[NX, NY, ve[2]]          self._numDataPoints=[NX, NY, ve[2]]
374          self._origin=origin          self._origin=origin
375            # we are rounding to avoid interpolateOnDomain issues
376          self._spacing=[np.round(lengths[i]/(self._numDataPoints[i]-1)) for i in xrange(3)]          self._spacing=[np.round(lengths[i]/(self._numDataPoints[i]-1)) for i in xrange(3)]
377          self._meshlen=[self._numDataPoints[i]*self._spacing[i] for i in xrange(3)]          self._meshlen=[(self._numDataPoints[i]-1)*self._spacing[i] for i in xrange(3)]
378          self._wkt_string=wkt_string          self._wkt_string=wkt_string
379          self._lon=lon_name          self._lon=lon_name
380          self._lat=lat_name          self._lat=lat_name
# Line 381  class NetCDFDataSource(DataSource): Line 388  class NetCDFDataSource(DataSource):
388    
389          self._meshlen=l_new          self._meshlen=l_new
390          self.NE=NE_new          self.NE=NE_new
391            self._dataorigin=self._origin
392          self._origin=origin_new          self._origin=origin_new
393          lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]          lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]
394          NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]          NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]
# Line 397  class NetCDFDataSource(DataSource): Line 405  class NetCDFDataSource(DataSource):
405              x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]              x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
406              mask=wherePositive(x[2]+self._origin[2])              mask=wherePositive(x[2]+self._origin[2])
407    
408          M=2 # do not constrain bottom          if self._constrainBottom:
409          #M=3 # constrain bottom              M=3 # constrain bottom
410            else:
411                M=2 # do not constrain bottom
412          for i in xrange(M):          for i in xrange(M):
413              mask=mask + whereNegative(x[i]) + \              mask=mask + whereNegative(x[i]) + \
414                      wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])                      wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])
# Line 440  class NetCDFDataSource(DataSource): Line 450  class NetCDFDataSource(DataSource):
450    
451      def _readGravity(self):      def _readGravity(self):
452          f=netcdf_file(self._gravfile, 'r')          f=netcdf_file(self._gravfile, 'r')
453          lon=f.variables[self._lon][:]          #lon=f.variables[self._lon][:]
454          lat=f.variables[self._lat][:]          #lat=f.variables[self._lat][:]
455            NE=[self._numDataPoints[i]-1 for i in xrange(2)]
456            lon=np.linspace(self._dataorigin[0], self._dataorigin[0]+NE[0]*self._spacing[0], NE[0]+1)
457            lat=np.linspace(self._dataorigin[1], self._dataorigin[1]+NE[1]*self._spacing[1], NE[1]+1)
458          lon,lat=np.meshgrid(lon,lat)          lon,lat=np.meshgrid(lon,lat)
         lon,lat=LatLonToUTM(lon, lat, self._wkt_string)  
459          grav=f.variables[self._grv][:]          grav=f.variables[self._grv][:]
460          f.close()          f.close()
461          lon=lon.flatten()          lon=lon.flatten()
# Line 452  class NetCDFDataSource(DataSource): Line 464  class NetCDFDataSource(DataSource):
464          alt=self._altOfData*np.ones(grav.shape)          alt=self._altOfData*np.ones(grav.shape)
465          # error value is an assumption          # error value is an assumption
466          try:          try:
467              missing=grav.missing_value              missing=f.variables[self._grv].missing_value
468              err=np.where(grav==missing, 20.0, 0.0)              err=np.where(grav==missing, 0.0, 20.0)
469          except:          except:
470              err=20.0*np.ones(lon.shape)              err=20.0*np.ones(lon.shape)
471          # convert units          # convert units
# Line 462  class NetCDFDataSource(DataSource): Line 474  class NetCDFDataSource(DataSource):
474          gravdata=np.column_stack((lon,lat,alt,grav,err))          gravdata=np.column_stack((lon,lat,alt,grav,err))
475          return gravdata          return gravdata
476    
477  class SmoothAnomaly(object):  ##############################################################################
478    class ERSDataSource(DataSource):
479        """
480        Data Source for ER Mapper raster data.
481        Note that this class only accepts a very specific type of ER Mapper data
482        input and will raise an exception if other data is found.
483        """
484        def __init__(self, domainclass, headerfile, datafile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):
485            """
486            headerfile - usually ends in .ers
487            datafile - usually has the same name as the headerfile without '.ers'
488            """
489            super(ERSDataSource,self).__init__()
490            self._headerfile=headerfile
491            if datafile is None:
492                self._datafile=headerfile[:-4]
493            else:
494                self._datafile=datafile
495            self._domainclass=domainclass
496            self._readHeader(vertical_extents)
497            self._altOfData=alt_of_data
498    
499        def _readHeader(self, ve):
500            metadata=open(self._headerfile, 'r').readlines()
501            # parse metadata
502            start=-1
503            for i in xrange(len(metadata)):
504                if metadata[i].strip() == 'DatasetHeader Begin':
505                    start=i+1
506            if start==-1:
507                raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
508    
509            md_dict={}
510            section=[]
511            for i in xrange(start, len(metadata)):
512                line=metadata[i]
513                if line[-6:].strip() == 'Begin':
514                    section.append(line[:-6].strip())
515                elif line[-4:].strip() == 'End':
516                    if len(section)>0:
517                        section.pop()
518                else:
519                    vals=line.split('=')
520                    if len(vals)==2:
521                        key = vals[0].strip()
522                        value = vals[1].strip()
523                        fullkey='.'.join(section+[key])
524                        md_dict[fullkey]=value
525    
526            try:
527                if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
528                    raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
529            except KeyError:
530                self.logger.warn("Cell type not specified. Assuming IEEE4ByteReal.")
531    
532            try:
533                NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
534                NY = int(md_dict['RasterInfo.NrOfLines'])
535            except:
536                raise RuntimeError("Could not determine extents of data")
537    
538            try:
539                maskval = float(md_dict['RasterInfo.NullCellValue'])
540            except:
541                raise RuntimeError("Could not determine NullCellValue")
542    
543            try:
544                spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
545                spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
546            except:
547                raise RuntimeError("Could not determine cell dimensions")
548    
549            try:
550                originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
551                originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
552            except:
553                self.logger.warn("Could not determine coordinate origin")
554                originX,originY = 0.0, 0.0
555    
556            # test data sets have origin in top-left corner so y runs top-down
557            self._dataorigin=[originX, originY]
558            originY-=(NY-1)*spacingY
559            self._maskval=maskval
560            self._spacing=[spacingX, spacingY, (ve[1]-ve[0])/(ve[2]-1)]
561            self._numDataPoints = [NX, NY, ve[2]]
562            self._origin = [originX, originY, ve[0]]
563            self._meshlen=[self._numDataPoints[i]*self._spacing[i] for i in xrange(3)]
564    
565        def getDensityMask(self):
566            return self._mask
567    
568        def getGravityAndStdDev(self):
569            gravlist=self._readGravity() # x,y,z,g,s
570            shape=[self.NE[2]+1, self.NE[1]+1, self.NE[0]+1]
571            g_and_sigma=self._interpolateOnDomain(gravlist, shape, self._origin, self._spacing, self._meshlen)
572            return g_and_sigma[0]*[0,0,1], g_and_sigma[1]
573    
574        def _createDomain(self, padding_l, padding_h):
575            NE=[self._numDataPoints[i]-1 for i in xrange(3)]
576            l=[NE[i]*self._spacing[i] for i in xrange(3)]
577            NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \
578                    NE, l, self._origin)
579    
580            self._meshlen=l_new
581            self.NE=NE_new
582            self._origin=origin_new
583            lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]
584            NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]
585            try:
586                dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])
587                x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
588                mask=wherePositive(dom.getX()[2])
589    
590            except TypeError:
591                dom=self._domainclass(*self.NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])
592                x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
593                mask=wherePositive(x[2]+self._origin[2])
594    
595            if self._constrainBottom:
596                M=3 # constrain bottom
597            else:
598                M=2 # do not constrain bottom
599            for i in xrange(M):
600                mask=mask + whereNegative(x[i]) + \
601                        wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])
602            self._mask=wherePositive(mask)
603    
604            self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
605            self.logger.debug("     length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))
606            self.logger.debug("     origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
607    
608            return dom
609    
610        def _readGravity(self):
611            f=open(self._datafile, 'r')
612            n=self._numDataPoints[0]*self._numDataPoints[1]
613            grav=[]
614            err=[]
615            for i in xrange(n):
616                v=struct.unpack('f',f.read(4))[0]
617                grav.append(v)
618                if abs((self._maskval-v)/v) < 1e-6:
619                    err.append(0.)
620                else:
621                    err.append(1.)
622            f.close()
623    
624            NE=[self._numDataPoints[i] for i in xrange(2)]
625            x=self._dataorigin[0]+np.arange(start=0, stop=NE[0]*self._spacing[0], step=self._spacing[0])
626            # test data sets have origin in top-left corner so y runs top-down
627            y=self._dataorigin[1]-np.arange(start=0, stop=NE[1]*self._spacing[1], step=self._spacing[1])
628            x,y=np.meshgrid(x,y)
629            x,y=x.flatten(),y.flatten()
630            alt=self._altOfData*np.ones((n,))
631    
632            # convert units
633            err=2e-5*np.array(err)
634            grav=1e-5*np.array(grav)
635            gravdata=np.column_stack((x,y,alt,grav,err))
636            return gravdata
637    
638    ##############################################################################
639    class SourceFeature(object):
640        """
641        A feature adds a density distribution to (parts of) a domain of a synthetic
642        data source, for example a layer of a specific rock type or a simulated
643        ore body.
644        """
645        def getDensity(self):
646            """
647            Returns the density for the area covered by mask. It can be constant
648            or a data object with spatial dependency.
649            """
650            raise NotImplementedError
651        def getMask(self, x):
652            """
653            Returns the mask of the area of interest for this feature. That is,
654            mask is non-zero where the density returned by getDensity() should be
655            applied, zero elsewhere.
656            """
657            raise NotImplementedError
658    
659    class SmoothAnomaly(SourceFeature):
660      def __init__(self, lx, ly, lz, x, y, depth, rho_inner, rho_outer):      def __init__(self, lx, ly, lz, x, y, depth, rho_inner, rho_outer):
661          self.x=x          self.x=x
662          self.y=y          self.y=y
# Line 473  class SmoothAnomaly(object): Line 667  class SmoothAnomaly(object):
667          self.rho_inner=rho_inner          self.rho_inner=rho_inner
668          self.rho_outer=rho_outer          self.rho_outer=rho_outer
669          self.rho=None          self.rho=None
670        
671        def getDensity(self):
672            return self.rho
673    
674      def getMask(self, x):      def getMask(self, x):
675          DIM=x.getDomain().getDim()          DIM=x.getDomain().getDim()
676          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]-(sup(x[DIM-1])-self.depth-self.lz/2)) * whereNonPositive(x[DIM-1]-(sup(x[DIM-1])-self.depth+self.lz/2)) \
# Line 487  class SmoothAnomaly(object): Line 685  class SmoothAnomaly(object):
685              if self.rho_inner<0: self.rho=-self.rho              if self.rho_inner<0: self.rho=-self.rho
686          return m              return m    
687    
688    ##############################################################################
689  class SyntheticDataSource(DataSource):  class SyntheticDataSource(DataSource):
690      def __init__(self, DIM, NE, l, h, features):      def __init__(self, DIM, NE, l, h, features):
691          super(SyntheticDataSource,self).__init__()          super(SyntheticDataSource,self).__init__()
# Line 539  class SyntheticDataSource(DataSource): Line 738  class SyntheticDataSource(DataSource):
738          rho_ref=0.          rho_ref=0.
739          for f in self._features:          for f in self._features:
740              m=f.getMask(self._x)              m=f.getMask(self._x)
741              rho_ref = rho_ref * (1-m) + f.rho * m              rho_ref = rho_ref * (1-m) + f.getDensity() * m
742          self._rho=rho_ref          self._rho=rho_ref
743    
744          return dom          return dom

Legend:
Removed from v.3957  
changed lines
  Added in v.3958

  ViewVC Help
Powered by ViewVC 1.1.26