glidertools.mapping.interp_obj

glidertools.mapping.interp_obj(x, y, z, xi, yi, partial_sill=0.1, nugget=0.01, lenscale_x=20, lenscale_y=20, detrend=True, max_points_per_quad=55, min_points_per_quad=8, return_error=False, n_cpus=None, verbose=True, parallel_chunk_size=512)

Performs objective interpolation (or Kriging) of a 2D field.

The objective interpolation breaks the problem into smaller fields by iteratively breaking the problem into quadrants. Each quadrant is then interpolated (also using intformation from its neighbours). The interpolation is inverse distance weighted using a gaussian kernel (or radial basis function). The kernel has a width of 12 hours if the x-dimension is time, otherwise scaled by the x-variable unit. The kernel is in meters assuming that depth is the y-coord. This can be changed with keyword arguements. An error estimate can also be calculated if requested.

The following link provides good background on the Kriging procedure: http://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/how-kriging-works.htm

Parameters:

x (np.array | pd.series) – horizontal coordinates of the input data (same length as y, z) can be types float or datetime64
y (np.array | pd.series) – vertical coordinates of the input data (same length as x, z)
z (np.array | pd.series) – values to be interoplated (same length as x, y)
xi (np.array) – horizontal coordinates of the interpolation grid (must be 1D) can be types float or datetime64
yi (np.array | pd.series) – vertical coordinates of the interpolation grid (must be 1D)
nugget (float [0.01]) – the error estimate due to sampling inaccuracy also known as the nugget in Kriging literature. This should be taken from the semivariogram
partial_sill (float [0.1]) – represents the spatial covariance of the variable being interpolated. Should be estimated from the semivariogram. See Kriging literature for more information
lenscale_x (float [20]) – horizontal length scale horizontal coordinate variable If dtype(x) is np.datetime64 (any format) then lenscale units is in hours. Otherwise if type(x).
lenscale_y (float [20]) – horizontal length scale horizontal coordinate variable.
max_points_per_quad (int [55]) – the data is divided into quadrants using a quadtree approach - iteratively dividing data into smaller quadrants using x and y coordinates. The algorithm stops splitting the data into quadrants when there are no quadrants exceeding the limit set with max_points_per_quad is. This is done to reduce the computational cost of the function.
min_points_per_quad (int [8]) – sets the minimum number of points allowed in a neighbouring quadrant when creating the interpolation function for a particular quadrant. If the number of points is less than specified, the algorithm looks for neighbours of the neighbours to include more points in the interpolation.
n_cpus (int [n - 1]) – use parallel computing. The quadrant calculations are spread across CPUs. Must be positive and > 0
parallel_chunk_size (int [512]) – the number of leaves that will be processed in parallel in one go. This is a memory saving feature. If your dataset is very large, parallel processing will use up a lot of memmory. Increasing the chunk size increases the memory requirements.
verbose (bool [True]) – will print out information about the interpolation

Returns:

Contains the following arrays: - z: interpolated values - variance: error estimate of the interpolation - weights: the quadtree weighting used to calculate the estimates - nugget: the nugget used in the interpolation - partial_sill: value used for the interpolation

Return type:

xr.Dataset

Note

The data may have semi-discrete artifacts. This is also present in the MATLAB output.

Example

>>> xi = np.arange(time.values.min(), time.values.max(), 30,
                   dtype='datetime64[m]')
>>> yi = np.linspace(depth.min(), depth.max(), 1.)
>>> interpolated = gt.mapping.interp_obj(
        time, depth, var, xi, yi,
        nugget=.0035, partial_sill=0.02,
        lenscale_x=80, lenscale_y=80,
        detrend=True)