Python scipy.interpolate.interpn() Examples

The following are 11 code examples of scipy.interpolate.interpn(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.interpolate , or try the search function .
Example #1
Source File: carbMesh.py    From badlands with GNU General Public License v3.0 6 votes vote down vote up
def _build_basement(self, tXY, bPts, regX, regY):
        """
        Using Pandas library to read the basement map file and define consolidated and
        soft sediment region.
        """
        self.tXY = tXY

        # Read basement file
        self.tinBase = numpy.ones(len(tXY))
        Bmap = pandas.read_csv(str(self.baseMap), sep=r'\s+', engine='c',
                    header=None, na_filter=False, dtype=numpy.float, low_memory=False)

        rectBase = numpy.reshape(Bmap.values,(len(regX), len(regY)),order='F')
        self.tinBase[bPts:] = interpolate.interpn( (regX, regY), rectBase,
                                                        tXY[bPts:,:], method='linear')

        return 
Example #2
Source File: elevationTIN.py    From badlands with GNU General Public License v3.0 6 votes vote down vote up
def getElevation(rX, rY, rZ, coords, interp="linear"):
    """
    This function interpolates elevation from the regular grid to the triamgular mesh
    using **SciPy** *interpn* funciton.

    Args:
        rX: numpy arrays containing the X coordinates from the regular grid.
        rY: numpy arrays containing the Y coordinates from the regular grid.
        rZ: numpy arrays containing the Z coordinates from the regular grid.
        coords: numpy float-type array containing X, Y coordinates for the TIN nodes.
        interp: interpolation method as in *SciPy interpn function* (default: 'linear')

    Returns:
        - elev - numpy array containing the updated elevations for the local domain.
    """

    # Set new elevation to 0
    elev = numpy.zeros(len(coords[:, 0]))

    # Get the TIN points elevation values using the regular grid dataset
    elev = interpn((rX, rY), rZ, (coords[:, :2]), method=interp)

    return elev 
Example #3
Source File: scaling.py    From CNNArt with Apache License 2.0 5 votes vote down vote up
def fScaleOnePatch(dPatch, randPatchSize, PatchSize):
    xaxis = np.linspace(0, PatchSize[0], randPatchSize[0])
    yaxis = np.linspace(0, PatchSize[1], randPatchSize[1])
    zaxis = np.linspace(0, PatchSize[2], randPatchSize[2])
    inter_train0 = np.mgrid[0:PatchSize[0], 0:PatchSize[1], 0:PatchSize[2]]
    inter_train1 = np.rollaxis(inter_train0, 0, 4)
    inter_train = np.reshape(inter_train1, [inter_train0.size // 3, 3])
    scaleddPatch = interpolate.interpn((xaxis, yaxis, zaxis), dPatch, inter_train, method='linear', bounds_error=False, fill_value=0)
    reshdPatch = np.reshape(scaleddPatch, [PatchSize[0], PatchSize[1], PatchSize[2]])
    return reshdPatch 
Example #4
Source File: ipol.py    From wradlib with MIT License 5 votes vote down vote up
def __call__(self, values, **kwargs):
        """Evaluate interpolator for values given at the source points.

        Parameters
        ----------
        values : :class:`numpy:numpy.ndarray` of floats
            Values at the source points which to interpolate, shape (num src pts, ...)

        Returns
        -------
        result : :class:`numpy:numpy.ndarray` of floats
            Target values with shape (num trg pts, ...)

        """

        # override bounds_error
        kwargs["bounds_error"] = kwargs.pop("bounds_error", False)
        kwargs["method"] = kwargs.pop("method", self.method)

        # need to flip ydim if grid origin is 'upper'
        if self.upper:
            values = np.flip(values, self.ydim)

        result = sinterp.interpn(self.ipol_grid, values, self.ipol_points, **kwargs)

        return result.reshape(self.points.shape[:-1]) 
Example #5
Source File: tabular.py    From Carnets with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def evaluate(self, *inputs):
        """
        Return the interpolated values at the input coordinates.

        Parameters
        ----------
        inputs : list of scalars or ndarrays
            Input coordinates. The number of inputs must be equal
            to the dimensions of the lookup table.
        """
        if isinstance(inputs, u.Quantity):
            inputs = inputs.value
        shape = inputs[0].shape
        inputs = [inp.flatten() for inp in inputs[: self.n_inputs]]
        inputs = np.array(inputs).T
        if not has_scipy:  # pragma: no cover
            raise ImportError("This model requires scipy >= v0.14")
        result = interpn(self.points, self.lookup_table, inputs,
                         method=self.method, bounds_error=self.bounds_error,
                         fill_value=self.fill_value)

        # return_units not respected when points has no units
        if (isinstance(self.lookup_table, u.Quantity) and
                not isinstance(self.points[0], u.Quantity)):
            result = result * self.lookup_table.unit

        if self.n_outputs == 1:
            result = result.reshape(shape)
        else:
            result = [r.reshape(shape) for r in result]
        return result 
Example #6
Source File: forceSim.py    From badlands with GNU General Public License v3.0 4 votes vote down vote up
def get_Rain(self, time, elev, inIDs):
        """
        Get rain value for a given period and perform interpolation from regular grid to unstructured TIN one.

        Parameters
        ----------
        time : float
            Requested time interval rain map to load.

        elev : float
            Unstructured grid (TIN) Z coordinates.

        inDs : integer
            List of unstructured vertices contained in each partition.

        Returns:
            - tinRain - numpy array containing the updated rainfall for the local domain.
        """

        events = numpy.where((self.T_rain[:, 1] - time) <= 0)[0]
        event = len(events)
        if not (time >= self.T_rain[event, 0]) and not (time < self.T_rain[event, 1]):
            raise ValueError("Problem finding the rain map to load!")

        if self.orographic[event]:
            if self.rzmax[event] <= 0:
                tinRain = self._build_OrographicRain_map(event, elev, inIDs)
                self.next_rain = min(time + self.ortime[event], self.T_rain[event, 1])
            else:
                tinRain = numpy.zeros(len(self.tXY[inIDs, 0]), dtype=float)
                rainslope = (self.rmax[event] - self.rmin[event]) / self.rzmax[event]
                tinRain = rainslope * elev[inIDs] + self.rmin[event]
                tinRain[tinRain < 0.0] = 0.0
                tinRain[tinRain > self.rmax[event]] = self.rmax[event]
                self.next_rain = min(time + self.ortime[event], self.T_rain[event, 1])
        elif self.Map_rain[event] == None:
            tinRain = numpy.zeros(len(self.tXY[inIDs, 0]), dtype=float)
            tinRain = self.rainVal[event]
            self.next_rain = self.T_rain[event, 1]
        else:
            rainMap = pandas.read_csv(
                str(self.Map_rain[event]),
                sep=r"\s+",
                engine="c",
                header=None,
                na_filter=False,
                dtype=numpy.float,
                low_memory=False,
            )

            rectRain = numpy.reshape(
                rainMap.values, (len(self.regX), len(self.regY)), order="F"
            )
            tinRain = interpolate.interpn(
                (self.regX, self.regY), rectRain, self.tXY[inIDs, :], method="linear"
            )
            self.next_rain = self.T_rain[event, 1]

        return tinRain 
Example #7
Source File: forceSim.py    From badlands with GNU General Public License v3.0 4 votes vote down vote up
def _build_OrographicRain_map(self, event, elev, inIDs):
        """
        Build rain map using Smith & Barstad (2004) model for a given period and perform interpolation from regular grid to unstructured TIN one.

        Args:
            event: float rain event number.
            elev : float unstructured grid (TIN) Z coordinates.
            inDs : integernumpy array of unstructured vertices contained in each partition.

        Returns:
            - tinRain - numpy array containing the updated rainfall for the local domain.
        """

        # Interpolate elevation on regular grid
        distances, indices = self.tree.query(self.xyi, k=8)
        if len(elev[indices].shape) == 3:
            elev_vals = elev[indices][:, :, 0]
        else:
            elev_vals = elev[indices]

        distances[distances < 0.0001] = 0.0001
        with numpy.errstate(divide="ignore"):
            oelev = numpy.average(elev_vals, weights=(1.0 / distances), axis=1)
        onIDs = numpy.where(distances[:, 0] <= 0.0001)[0]
        if len(onIDs) > 0:
            oelev[onIDs] = elev[indices[onIDs, 0]]
        oelev -= self.sealevel
        oelev = oelev.clip(0)
        regZ = numpy.reshape(oelev, (len(self.regX), len(self.regY)), order="F")
        # Use Smith & Barstad model
        rectRain = ormodel.compute(
            regZ,
            self.dx,
            self.windx[event],
            self.windy[event],
            self.rmin[event],
            self.rmax[event],
            self.rbgd[event],
            self.nm[event],
            self.cw[event],
            self.hw[event],
            self.tauc[event],
            self.tauf[event],
        )

        # Apply smoothing here
        smthRain = gaussian_filter(rectRain, sigma=3)

        # Interpolate
        tinRain = interpolate.interpn(
            (self.regX, self.regY), smthRain, self.tXY[inIDs, :], method="linear"
        )

        return tinRain 
Example #8
Source File: forceSim.py    From badlands with GNU General Public License v3.0 4 votes vote down vote up
def load_Tecto_map(self, time, inIDs):
        """
        Load vertical displacement map for a given period and perform interpolation from regular grid to unstructured TIN one.

        Args:
            time : float requested time interval rain map to load.
            inDs : integer numpy array of unstructured vertices contained in each partition.

        Returns:
            - tinDisp - numpy array containing the updated displacement rate for the local domain.

        """

        events = numpy.where((self.T_disp[:, 1] - time) <= 0)[0]
        event = len(events)

        if not (time >= self.T_disp[event, 0]) and not (time < self.T_disp[event, 1]):
            raise ValueError("Problem finding the displacements map to load!")

        self.next_disp = self.T_disp[event, 1]

        if self.injected_disps is not None or self.Map_disp[event] != None:
            if self.injected_disps is not None:
                dispMap = self.injected_disps
            else:
                dispMap = pandas.read_csv(
                    str(self.Map_disp[event]),
                    sep=r"\s+",
                    engine="c",
                    header=None,
                    na_filter=False,
                    dtype=numpy.float,
                    low_memory=False,
                ).values

            rectDisp = numpy.reshape(
                dispMap, (len(self.regX), len(self.regY)), order="F"
            )
            tinDisp = interpolate.interpn(
                (self.regX, self.regY), rectDisp, self.tXY[inIDs, :], method="linear"
            )
            dt = self.T_disp[event, 1] - self.T_disp[event, 0]
            if dt <= 0:
                raise ValueError(
                    "Problem computing the displacements rate for event %d." % event
                )
            tinDisp = tinDisp / dt
        else:
            tinDisp = numpy.zeros(len(self.tXY[inIDs, 0]), dtype=float)

        return tinDisp 
Example #9
Source File: scaling.py    From CNNArt with Apache License 2.0 4 votes vote down vote up
def fscaling3D(X_train, X_test, scpatchSize, iscalefactor):
    afterSize = np.ceil(np.multiply(scpatchSize, iscalefactor)).astype(int)

    # Prepare for the using of scipy.interpolation: create the coordinates of grid
    if iscalefactor == 1:
        return X_train, X_test, scpatchSize
    else:
        xaxis = np.linspace(0, afterSize[0], scpatchSize[0])
        yaxis = np.linspace(0, afterSize[1], scpatchSize[1])
        zaxis = np.linspace(0, afterSize[2], scpatchSize[2])

    dAllx_train = None
    dAllx_test = None

    for ifold in range(len(X_train)):
        lenTrain = X_train[ifold].shape[0]
        lenTest = X_test[ifold].shape[0]

        start = time.clock()

        # no batch
        inter_train0 = np.mgrid[0:lenTrain, 0:afterSize[0], 0:afterSize[1], 0:afterSize[2]]
        inter_train1 = np.rollaxis(inter_train0, 0, 5)
        inter_train = np.reshape(inter_train1, [inter_train0.size // 4, 4])  # 4 for the dimension of coordinates

        zaxis_train = np.arange(lenTrain)

        upedTrain = interpolate.interpn((zaxis_train, xaxis, yaxis, zaxis),
                                        X_train[ifold],
                                        inter_train, method='linear', bounds_error=False, fill_value=0)
        dFoldx_train = np.reshape(upedTrain, [1, lenTrain, afterSize[0], afterSize[1], afterSize[2]])


        inter_test0 = np.mgrid[0:lenTest, 0:afterSize[0], 0:afterSize[1], 0:afterSize[2]]
        inter_test1 = np.rollaxis(inter_test0, 0, 5)
        inter_test = np.reshape(inter_test1, [inter_test0.size // 4, 4])  # 4 for the dimension of coordinates

        zaxis_test = np.arange(lenTest)

        upedTest = interpolate.interpn((zaxis_test, xaxis, yaxis, zaxis),
                                       X_test[ifold],
                                       inter_test, method='linear', bounds_error=False, fill_value=0)
        dFoldx_test = np.reshape(upedTest, [1, lenTest, afterSize[0], afterSize[1], afterSize[2]])

        stop = time.clock()
        print(stop-start)

        if dAllx_train is None:
            dAllx_train = dFoldx_train
        else:
            dAllx_train = np.concatenate((dAllx_train, dFoldx_train), axis=0)

        if dAllx_test is None:
            dAllx_test = dFoldx_test
        else:
            dAllx_test = np.concatenate((dAllx_test, dFoldx_test), axis=0)

    return dAllx_train, dAllx_test, afterSize 
Example #10
Source File: DEM.py    From sarpy with MIT License 4 votes vote down vote up
def elevate(self, coord, method='linear',lonlat=False):
        '''
        Given a geographic coordinate, return an elevation from the DEM
        Coord may be a sinlge coord (c  = [lat,lon]) or multiples (c = [[lat,lon],[lat,lon],...]
        '''
        if not 'dem' in dir(self):  # Check to see if dem is an attribute
            self.logger.warning('There are no DEMs to interpolate from.')  # If not, no dems have been read in
            return np.zeros(coord.shape[0]) + -1234.5
        if np.max(self.origindd[:,1]) == 179.0 and np.min(self.origindd[:,1]) == -180.0:  # If DEM range crosses the dateline
            coord[(coord < 0)] += 360.0  # Convert to -180 to 180 -> 0 to 360
        # interpolated_elevation = interpn((1-D LON array, 1-D LAT array), 2-D elev array, coord array)
        coord = np.array(coord)  # Ensure coord is an array for the following line to work
        if len(coord.shape) == 1:  # If only one point is specified
            coord = coord.reshape((1,2))  # Then make sure its a 2-D array of 1 by 2 [[x,y]]
        if not lonlat:  # If coords are given at LAT/LON (y,x), then switch to LON/LAT (x,y)
            coord = self.geo_swap(coord)  # Convert [[lat,lon],...] to [[lon, lat], ...]
        # The following is to ensure interpn evaluates all the good, valid coordinates instead
        # of throwing the baby out with the bath water.
        elev = np.zeros(coord.shape[0]) + -1234.5  # Create a list of dummy elevations the same length as input list
        elev2 = elev
        # Get all valid coordinates
        in_bounds = np.where((coord[:,1] > np.min(self.lats_1D)) & (coord[:,1] < np.max(self.lats_1D)) &
                             (coord[:,0] > np.min(self.lons_1D)) & (coord[:,0] < np.max(self.lons_1D)))[0]
        
        self.logger.debug('Coord LAT range: {}'.format((np.min(coord[:,1]),np.max(coord[:,1]))))
        self.logger.debug('Coord LON range: {}'.format((np.min(coord[:,0]),np.max(coord[:,0]))))
        self.logger.debug('DEM LAT range: {}'.format((np.min(self.lats_1D),np.max(self.lats_1D))))
        self.logger.debug('DEM LON range: {}'.format((np.min(self.lons_1D),np.max(self.lons_1D))))
        self.logger.debug('Coord shape: {}'.format(coord.shape))
        self.logger.debug('Elev size: {}'.format(elev.size))
        self.logger.debug('In_bounds size: {}'.format(in_bounds.size))
        if in_bounds.size < elev.size:
            self.logger.warning('Some points may be outside of DEM boundary. Check coordinate list.')
        if in_bounds.size > 0:  # If there are any valid points, then do try the interpolation
            try:
                self.logger.info('Interpolating elevation points.')
                elev[in_bounds] = interpn((self.lons_1D, self.lats_1D), self.dem, coord[in_bounds,:], method=method)
                #f = spint.interp2d(self.lon_1D, self.lats_1D, self.dem)
                #elev2[in_bounds] = f(coord[in_bounds,:])
            except Exception as err:
                self.logger.critical('Interpolation error: {}'.format(err))
        good_heights = np.where(elev > -1234.5)[0]  # Do stats on valid points only.
        if good_heights.size > 0:  # If there are good points then print stats
            emin = np.round(np.min(elev[good_heights]),2)
            emean = np.round(np.mean(elev[good_heights]),2)
            emax = np.round(np.max(elev[good_heights]),2)
            self.logger.info('Elevation stats (min/mean/max): {}/{}/{}'.format(emin,emean,emax))
        else:
            self.logger.info('No valid points found.')
        return elev  # END OF ELEVATE 
Example #11
Source File: atomic_diffraction_generator_utils.py    From diffsims with GNU General Public License v3.0 4 votes vote down vote up
def grid2sphere(arr, x, dx, C):
    """
    Projects 3d array onto a sphere

    Parameters
    ----------
    arr : `numpy.ndarray` [`float`], (nx, ny, nz)
        Input function to be projected
    x : `list` [`numpy.ndarray` [`float`]], of shapes [(nx,), (ny,), (nz,)]
        Vectors defining mesh of <arr>
    dx : `list` [`numpy.ndarray` [`float`]], of shapes [(3,), (3,), (3,)]
        Basis in which to orient sphere. Centre of sphere will be at `C*dx[2]`
        and mesh of output array will be defined by the first two vectors
    C : `float`
        Radius of sphere

    Returns
    -------
    out : `numpy.ndarray` [`float`], (nx, ny)
        If y is the point on the line between `i*dx[0]+j*dx[1]` and `C*dx[2]`
        which also lies on the sphere of radius `C` from `C*dx[2]` then:
            `out[i,j] = arr(y)`
        Interpolation on arr is linear.
    """
    if C in (None, 0) or x[2].size == 1:
        if arr.ndim == 2:
            return arr
        elif arr.shape[2] == 1:
            return arr[:, :, 0]

    y = to_mesh((x[0], x[1], array([0])), dx).reshape(-1, 3)

    if C is not None:  # project on line to centre
        w = 1 / (1 + (y ** 2).sum(-1) / C ** 2)
        y *= w[:, None]
        if dx is None:
            y[:, 2] = C * (1 - w)
        else:
            y += C * (1 - w)[:, None] * dx[2]

    out = interpn(x, arr, y, method="linear", bounds_error=False, fill_value=0)

    return out.reshape(x[0].size, x[1].size)