/[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 4361 by gross, Wed Mar 27 07:58:34 2013 UTC revision 4362 by caltinay, Tue Apr 16 04:37:13 2013 UTC
# Line 22  __license__="""Licensed under the Open S Line 22  __license__="""Licensed under the Open S
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', 'ErMapperData', \  __all__ = ['DataSource', 'ErMapperData', 'NumpyData',\
26          'SyntheticDataBase', 'SyntheticFeatureData', 'SyntheticData',          'SyntheticDataBase', 'SyntheticFeatureData', 'SyntheticData',
27          'SmoothAnomaly']          'SmoothAnomaly']
28    
29  import logging  import logging
30  import numpy as np  import numpy as np
31    import tempfile
32  from esys.escript import ReducedFunction  from esys.escript import ReducedFunction
33  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
34  from esys.escript.linearPDEs import LinearSinglePDE  from esys.escript.linearPDEs import LinearSinglePDE
# Line 975  class SyntheticData(SyntheticDataBase): Line 976  class SyntheticData(SyntheticDataBase):
976              self._reference_data= k              self._reference_data= k
977          return self._reference_data          return self._reference_data
978    
979    ##############################################################################
980    class NumpyData(DataSource):
981        """
982        """
983    
984        def __init__(self, data_type, data, error=1., length=1.*U.km, null_value=-1.):
985            """
986            A data source that uses survey data from a ``numpy`` object or list
987            instead of a file.
988            The dimensionality is inferred from the shape of ``data`` (1- and
989            2-dimensional data is supported). The data origin is assumed to be
990            at the coordinate origin.
991    
992            :param data_type: data type indicator
993            :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
994            :param data: the survey data array. Note that for a cartesian coordinate
995                         system the shape of the data is considered to be
996                         (nz,ny,nx).
997            :type data: ``numpy.array`` or ``list``
998            :param error: constant value to use for the data uncertainties or a
999                          numpy object with uncertainties for every data point.
1000            :type error: ``float`` or ``list`` or ``ndarray``
1001            :param length: side length(s) of the data slice/volume. This can be
1002                           a scalar to indicate constant length in all dimensions
1003                           or an array/list of values in each coordinate dimension.
1004            :type length: ``float`` or ``list`` or ``ndarray``
1005            :param null_value: value that is used in the undefined regions of the
1006                               survey data object.
1007            :type null_value: ``float``
1008            """
1009            super(NumpyData, self).__init__()
1010            if not data_type in [self.GRAVITY, self.MAGNETIC]:
1011                raise ValueError("Invalid value for data_type parameter")
1012            self.__data_type = data_type
1013            self.__data = np.asarray(data, dtype=np.float32)
1014            DIM = len(self.__data.shape)
1015            if DIM not in (1,2):
1016                raise ValueError("NumpyData requires 1- or 2-dimensional data")
1017            self.__error = np.asarray(error, dtype=np.float32)
1018            if len(self.__error.shape) > 0 and \
1019                    self.__error.shape != self.__data.shape:
1020                raise ValueError("error argument must be scalar or match the shape of the data")
1021            if isinstance(length,float) or isinstance(length,int):
1022                length = [float(length)] * DIM
1023            else:
1024                length=np.asarray(length, dtype=float)
1025                if len(length.shape) != 1 or length.shape[0] != DIM:
1026                    raise ValueError("length must be scalar or an array with %d values"%DIM)
1027                length=length.tolist()
1028    
1029            self.__null_value = null_value
1030            self.__nPts = list(reversed(self.__data.shape))
1031            self.__delta = [length[i]/self.__nPts[i] for i in range(DIM)]
1032            self.__origin = [0.] * DIM # this could be an optional argument
1033    
1034        def getDataExtents(self):
1035            return (self.__origin, self.__nPts, self.__delta)
1036    
1037        def getDataType(self):
1038            return self.__data_type
1039    
1040        def getSurveyData(self, domain, origin, NE, spacing):
1041            FS = ReducedFunction(domain)
1042            nValues = self.__nPts
1043            dataDIM = len(nValues)
1044            # determine base location of this dataset within the domain
1045            first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(dataDIM)]
1046            # determine the resolution difference between domain and data.
1047            # If domain has twice the resolution we can double up the data etc.
1048            multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(dataDIM)]
1049            if domain.getDim() > dataDIM:
1050                first.append(int((-origin[-1])/spacing[-1]))
1051                multiplier=multiplier+[1]
1052                nValues=nValues+[1]
1053    
1054            _handle, numpyfile = tempfile.mkstemp()
1055            os.close(_handle)
1056            self.__data.tofile(numpyfile)
1057    
1058            data = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1059                                             multiplier, (), self.__null_value)
1060            if len(self.__error.shape) > 0:
1061                self.__error.tofile(numpyfile)
1062                sigma = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1063                                                 multiplier, (), 0.)
1064            else:
1065                sigma = self.__error.item() * whereNonZero(data-self.__null_value)
1066    
1067            os.unlink(numpyfile)
1068    
1069            return data, sigma
1070    
1071        def getUtmZone(self):
1072            """
1073            returns a dummy UTM zone since this class does not use real coordinate
1074            values.
1075            """
1076            return 0
1077    

Legend:
Removed from v.4361  
changed lines
  Added in v.4362

  ViewVC Help
Powered by ViewVC 1.1.26