/[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 4196 by caltinay, Wed Feb 13 00:58:37 2013 UTC revision 4197 by caltinay, Wed Feb 13 05:53:19 2013 UTC
# Line 32  from esys.escript import ReducedFunction Line 32  from esys.escript import ReducedFunction
32  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
33  from esys.escript.linearPDEs import LinearSinglePDE  from esys.escript.linearPDEs import LinearSinglePDE
34  from esys.escript.util import *  from esys.escript.util import *
35  from esys.ripley import Brick, Rectangle, ripleycpp  from esys.ripley import ripleycpp
36    
37  try:  try:
38      from scipy.io.netcdf import netcdf_file      from scipy.io.netcdf import netcdf_file
# Line 45  def LatLonToUTM(lon, lat, wkt_string=Non Line 45  def LatLonToUTM(lon, lat, wkt_string=Non
45      Converts one or more longitude,latitude pairs to the corresponding x,y      Converts one or more longitude,latitude pairs to the corresponding x,y
46      coordinates in the Universal Transverse Mercator projection.      coordinates in the Universal Transverse Mercator projection.
47    
48      :note: If the ``pyproj`` module is not installed a warning is printed and      :note: The ``pyproj`` module is required for this method. If it is not
49             the input values are scaled by a constant and returned.             found an exception is raised.
50      :note: If `wkt_string` is not given or invalid or the ``gdal`` module is      :note: If `wkt_string` is not given or invalid or the ``gdal`` module is
51             not available to convert the string, then the input values are             not available to convert the string, then the input values are
52             assumed to be given in the Clarke 1866 projection.             assumed to be using the Clarke 1866 ellipsoid.
53    
54      :param lon: longitude value(s)      :param lon: longitude value(s)
55      :type lon: ``float``, ``list``, ``tuple``, or ``numpy.array``      :type lon: ``float``, ``list``, ``tuple``, or ``numpy.array``
56      :param lat: latitude value(s)      :param lat: latitude value(s)
57      :type lat: ``float``, ``list``, ``tuple``, or ``numpy.array``      :type lat: ``float``, ``list``, ``tuple``, or ``numpy.array``
58        :param wkt_string: Well-known text (WKT) string describing the coordinate
59                           system used. The ``gdal`` module is used to convert
60                           the string to the corresponding Proj4 string.
61        :type wkt_string: ``str``
62      :rtype: ``numpy.array``      :rtype: ``numpy.array``
63      """      """
64    
     # not really optimal: if pyproj is not installed we return the input  
     # values scaled by a constant.  
65      try:      try:
66          import pyproj          import pyproj
67      except:      except:
68          print("Warning, pyproj not available. Domain extents will be wrong")          raise ImportError("The pyproj module is required for coordinate conversion. Please install and try again.")
         return np.array(lon)*1000., np.array(lat)*1000.  
69    
70      # determine UTM zone from the input data      # determine UTM zone from the input data
71      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 75  def LatLonToUTM(lon, lat, wkt_string=Non Line 76  def LatLonToUTM(lon, lat, wkt_string=Non
76          p_src = pyproj.Proj(srs.ExportToProj4())          p_src = pyproj.Proj(srs.ExportToProj4())
77      except:      except:
78          p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')          p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
79      # we assume southern hemisphere here  
80      p_dest = pyproj.Proj('+proj=utm +zone=%d +south +units=m'%zone)      # check for hemisphere
81        if np.median(np.array(lat)) < 0.:
82            south='+south '
83        else:
84            south=''
85        p_dest = pyproj.Proj('+proj=utm +zone=%d %s+units=m'%(zone,south))
86      x,y=pyproj.transform(p_src, p_dest, lon, lat)      x,y=pyproj.transform(p_src, p_dest, lon, lat)
87      return x,y      return x,y
88    
# Line 326  class NetCdfData(DataSource): Line 332  class NetCdfData(DataSource):
332      """      """
333      Data Source for gridded netCDF data that use CF/COARDS conventions.      Data Source for gridded netCDF data that use CF/COARDS conventions.
334      """      """
335      def __init__(self, datatype, filename, altitude=0.):      def __init__(self, data_type, filename, altitude=0., data_variable=None,
336                           error=None, scale_factor=None, null_value=None):
337          """          """
338          :param filename: file name for survey data in netCDF format          :param filename: file name for survey data in netCDF format
339          :type filename: ``str``          :type filename: ``str``
340          :param datatype: type of data, must be `GRAVITY` or `MAGNETIC`          :param data_type: type of data, must be `GRAVITY` or `MAGNETIC`
341          :type datatype: ``int``          :type data_type: ``int``
342          :param altitude: altitude of measurements in meters          :param altitude: altitude of measurements in meters
343          :type altitude: ``float``          :type altitude: ``float``
344            :param data_variable: name of the netCDF variable that holds the data.
345                                  If not provided an attempt is made to determine
346                                  the variable and an exception thrown on failure.
347            :type data_variable: ``str``
348            :param error: either the name of the netCDF variable that holds the
349                          uncertainties of the measurements or a constant value
350                          to use for the uncertainties. If a constant value is
351                          supplied, it is scaled by the same factor as the
352                          measurements. If not provided the error is assumed to
353                          be 2 units for all measurements (i.e. 0.2 mGal and 2 nT
354                          for gravity and magnetic, respectively)
355            :type error: ``str`` or ``float``
356            :param scale_factor: the measurements and error values are scaled by
357                                 this factor. By default, gravity data is assumed
358                                 to be given in 1e-6 m/s^2 (0.1 mGal), while
359                                 magnetic data is assumed to be in 1e-9 T (1 nT).
360            :type scale_factor: ``float``
361            :param null_value: value that is used in the file to mark undefined
362                               areas. This information is usually included in the
363                               file.
364            :type null_value: ``float``
365          """          """
366          super(NetCdfData,self).__init__()          super(NetCdfData,self).__init__()
367          self.__filename=filename          self.__filename=filename
368          if not datatype in [self.GRAVITY,self.MAGNETIC]:          if not data_type in [self.GRAVITY,self.MAGNETIC]:
369              raise ValueError("Invalid value for datatype parameter")              raise ValueError("Invalid value for datatype parameter")
370          self.__datatype=datatype          self.__data_type = data_type
371          self.__altitude=altitude          self.__altitude = altitude
372          self.__readMetadata()          self.__data_name = data_variable
373            self.__scale_factor = scale_factor
374            self.__null_value = null_value
375            self.__readMetadata(error)
376    
377      def __readMetadata(self):      def __readMetadata(self, error):
378          self.logger.debug("Checking Data Source: %s"%self.__filename)          self.logger.debug("Checking Data Source: %s"%self.__filename)
379          f=netcdf_file(self.__filename, 'r')          f=netcdf_file(self.__filename, 'r')
380    
381            ### longitude- / X-dimension and variable
382          NX=0          NX=0
383          for n in ['lon','longitude','x']:          for n in ['lon','longitude','x']:
384              if n in f.dimensions:              if n in f.dimensions:
385                  NX=f.dimensions[n]                  NX=f.dimensions[n]
386                    lon_name=n
387                  break                  break
388          if NX==0:          if NX==0:
389              raise RuntimeError("Could not determine extents of data")              raise RuntimeError("Could not determine extents of data")
390    
391            # CF/COARDS states that coordinate variables have the same name as
392            # the dimensions
393            if not lon_name in f.variables:
394                raise RuntimeError("Could not determine longitude variable")
395            longitude=f.variables.pop(lon_name)
396    
397            ### latitude- / Y-dimension and variable
398          NY=0          NY=0
399          for n in ['lat','latitude','y']:          for n in ['lat','latitude','y']:
400              if n in f.dimensions:              if n in f.dimensions:
401                  NY=f.dimensions[n]                  NY=f.dimensions[n]
402                    lat_name=n
403                  break                  break
404          if NY==0:          if NY==0:
405              raise RuntimeError("Could not determine extents of data")              raise RuntimeError("Could not determine extents of data")
406            if not lat_name in f.variables:
         # 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:  
407              raise RuntimeError("Could not determine latitude variable")              raise RuntimeError("Could not determine latitude variable")
408            latitude=f.variables.pop(lat_name)
409    
410          # try to figure out data variable name          ### data variable
411          data_name=None          if self.__data_name is not None:
412          if len(f.variables)==1:              try:
413              data_name=f.variables.keys()[0]                  dims = f.variables[self.__data_name].dimensions
414                    if not ((lat_name in dims) and (lon_name in dims)):
415                        raise ValueError("Invalid data variable name supplied")
416                except KeyError:
417                    raise ValueError("Invalid data variable name supplied")
418          else:          else:
419              for n in f.variables.keys():              for n in f.variables.keys():
420                  dims=f.variables[n].dimensions                  dims=f.variables[n].dimensions
421                  if (lat_name in dims) and (lon_name in dims):                  if (lat_name in dims) and (lon_name in dims):
422                      data_name=n                      self.__data_name=n
423                      break                      break
424          if data_name is None:          if self.__data_name is None:
425              raise RuntimeError("Could not determine data variable")              raise RuntimeError("Could not determine data variable")
426    
427          # try to determine value for unused data          datavar = f.variables[self.__data_name]
428          # NaN is always filtered out in ripley  
429          if hasattr(f.variables[data_name], 'missing_value'):          ### error value/variable
430              maskval = float(f.variables[data_name].missing_value)          self.__error_name = None
431          elif hasattr(f.variables[data_name], '_FillValue'):          if isinstance(error,str):
432              maskval = float(f.variables[data_name]._FillValue)              try:
433                    dims = f.variables[error].dimensions
434                    if not ((lat_name in dims) and (lon_name in dims)):
435                        raise ValueError("Invalid error variable name supplied")
436                except KeyError:
437                    raise ValueError("Invalid error variable name supplied")
438                self.__error_name = error
439            elif isinstance(error,float) or isinstance(error,int):
440                self.__error_value = float(error)
441            elif error is None:
442                self.__error_value = 2.
443          else:          else:
444              self.logger.debug("missing_value attribute not found, using default.")              raise TypeError("Invalid type of error parameter")
             maskval = 99999  
445    
446          # try to determine units of data - this is disabled for now          ### mask/null value
447            # note that NaN is always filtered out in ripley
448            if self.__null_value is None:
449                if hasattr(datavar, 'missing_value'):
450                    self.__null_value = float(datavar.missing_value)
451                elif hasattr(datavar, '_FillValue'):
452                    self.__null_value = float(datavar._FillValue)
453                else:
454                    self.logger.debug("Could not determine null value, using default.")
455                    self.__null_value = 99999
456    
457            # try to determine units of data - this is disabled until we obtain a
458            # file with valid information
459          #if hasattr(f.variables[data_name], 'units'):          #if hasattr(f.variables[data_name], 'units'):
460          #   units=f.variables[data_name].units          #   units=f.variables[data_name].units
461          if self.__datatype == self.GRAVITY:  
462              self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")          ### scale factor
463              self.__scalefactor = 1e-6          if self.__scale_factor is None:
464          else:              if self.__data_type == self.GRAVITY:
465              self.logger.info("Assuming magnetic data units are 'nT'.")                  self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
466              self.__scalefactor = 1e-9                  self.__scale_factor = 1e-6
467                else:
468                    self.logger.info("Assuming magnetic data units are 'nT'.")
469                    self.__scale_factor = 1e-9
470    
471          # see if there is a wkt string to convert coordinates          # see if there is a wkt string to convert coordinates
472          try:          try:
473              wkt_string=f.variables[data_name].esri_pe_string              wkt_string=datavar.esri_pe_string
474          except:          except:
475              wkt_string=None              wkt_string=None
476    
477          # we don't trust actual_range & geospatial_lon_min/max since subset          # actual_range & geospatial_lon_min/max do not always contain correct
478          # data does not seem to have these fields updated.          # values so we have to obtain the min/max in a less efficient way:
         # Getting min/max from the arrays is obviously not very efficient but..  
         #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]  
479          lon_range=longitude.data.min(),longitude.data.max()          lon_range=longitude.data.min(),longitude.data.max()
480          lat_range=latitude.data.min(),latitude.data.max()          lat_range=latitude.data.min(),latitude.data.max()
481          if lon_range[1]<180:          lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
             lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)  
482          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]          lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
483          f.close()          f.close()
484    
# Line 436  class NetCdfData(DataSource): Line 486  class NetCdfData(DataSource):
486          self.__origin=[lon_range[0],lat_range[0]]          self.__origin=[lon_range[0],lat_range[0]]
487          # we are rounding to avoid interpolation issues          # we are rounding to avoid interpolation issues
488          self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(2)]          self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(2)]
         self.__data_name=data_name  
         self.__maskval=maskval  
489    
490      def getDataExtents(self):      def getDataExtents(self):
491          """          """
# Line 446  class NetCdfData(DataSource): Line 494  class NetCdfData(DataSource):
494          return (list(self.__origin), list(self.__nPts), list(self.__delta))          return (list(self.__origin), list(self.__nPts), list(self.__delta))
495    
496      def getDataType(self):      def getDataType(self):
497          return self.__datatype          return self.__data_type
498    
499      def getSurveyData(self, domain, origin, NE, spacing):      def getSurveyData(self, domain, origin, NE, spacing):
500            FS=ReducedFunction(domain)
501          nValues=self.__nPts          nValues=self.__nPts
502          # determine base location of this dataset within the domain          # determine base location of this dataset within the domain
503          first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]          first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
# Line 456  class NetCdfData(DataSource): Line 505  class NetCdfData(DataSource):
505              first.append(int((self.__altitude-origin[2])/spacing[2]))              first.append(int((self.__altitude-origin[2])/spacing[2]))
506              nValues=nValues+[1]              nValues=nValues+[1]
507    
508          data=ripleycpp._readNcGrid(self.__filename, self.__data_name,          data = ripleycpp._readNcGrid(self.__filename, self.__data_name, FS,
509                    ReducedFunction(domain), first, nValues, (), self.__maskval)                                       first, nValues, (), self.__null_value)
510          sigma=whereNonZero(data-self.__maskval)          if self.__error_name is not None:
511          data=data*self.__scalefactor              sigma = ripleycpp._readNcGrid(self.__filename, self.__error_name,
512          sigma=sigma * 2. * self.__scalefactor                                            FS, first, nValues, (), 0.)
513            else:
514                sigma = self.__error_value * whereNonZero(data-self.__null_value)
515    
516            data = data * self.__scale_factor
517            sigma = sigma * self.__scale_factor
518          return data, sigma          return data, sigma
519    
520    

Legend:
Removed from v.4196  
changed lines
  Added in v.4197

  ViewVC Help
Powered by ViewVC 1.1.26