Python scipy.interpolate.interpn() Examples

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 

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 

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]) 

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 

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 

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 

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 

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 

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 

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 

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)