Example - 2D Grid Interpolation

The xmsinterp library contains classes for performing 2d (x,y) linear, natural neighbor, and idw interpolation: InterpLinear and InterpIdw. Natural neighbor interpolation is performed using the linear interpolation class. 3d (x,y,z) interpolation can be performed using the idw class. Points (x,y) must be given to create the interpolation class.

We will learn how to use the interpolation classes by using a set of 2D points (x,y) with a measurement at each point. These points come from a csv file. Then we will generate a grid that bounds these points and use different interpolation options to visualize the data.

In [1]:
import param
import panel as pn
import holoviews as hv
import numpy as np
import xmsinterp as interp
import pandas as pd
import xarray as xr
import hvplot.pandas
import hvplot.xarray
from holoviews.streams import Params
In [2]:
# read in the data to do the interpolation
# drop the id columns from the data frame so we can pass the remaining data as points to the interp class
df = pd.read_csv("plumedat.csv").drop(['id'], axis=1)
df.hvplot.scatter(x='x', y='y', c='c', height=600, cmap='jet', size=120, colorbar=True)
Out[2]:

Generate grid for visualization

Get the min, max (x, y) of the data and generate a grid of points.

In [3]:
# get the (x,y) range of the data and expand by 10%
def buffer_range(range_min, range_max, buf_by=0.1):
    buf = (range_max - range_min) * buf_by
    return (range_min - buf, range_max + buf)

x_range = buffer_range(df.x.min(), df.x.max())
y_range = buffer_range(df.y.min(), df.y.max())

# adjust these values for a higher/lower resolution grid
number_of_x_points=600
number_of_y_points=400
x = np.linspace(x_range[0], x_range[1], number_of_x_points)
y = np.linspace(y_range[0], y_range[1], number_of_y_points)
grid = np.meshgrid(x, y)
pts = np.stack([i.flatten() for i in grid], axis=1)

Interpolate to the grid using linear interpolation and display the results with matplotlib.

In [4]:
linear = interp.interpolate.InterpLinear(df.values)
z = linear.interp_to_pts(pts)
# I would prefer to be able to set up an interpolator like this:
# linear = interp.Interpolator(kind='linear', values=df.values)
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[4]:

Interpolate using natural neighbor. The first example shows a "constant" nodal function and the second example shows a "quadratic" nodals function. A nodal function can better model trends that exist in the measured values.

In [5]:
linear.set_use_natural_neighbor()
z = linear.interp_to_pts(pts)
# z = linear.interp_to_pts(pts, use='natural_neighbor')
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[5]:

Now do natural neighbor with a quadratic nodal function.

In [6]:
linear.set_use_natural_neighbor(nodal_func_type="quadratic", nd_func_pt_search_opt="nearest_pts")
z = linear.interp_to_pts(pts)
# z = linear.interp_to_pts(pts, use='natural_neighbor', nodal_func_type='quadratic', nd_func_pt_search_opt='nearest_pts')
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[6]:

Now we will use Inverse Distance Weighted (IDW) to interpolate to the grid. Notice that unlike Linear and Natural Neighbor interpolation, IDW interpolation can extrapolate beyond the bounds of the input data points.

In [7]:
idw = interp.interpolate.InterpIdw(df.values)
z = idw.interp_to_pts(pts)
# idw = interp.Interpolator(kind='idw', values=df.values)
# z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[7]:

Now interpolate using a nodal function to try to better capture trends in the data set.

In [8]:
idw.set_nodal_function(nodal_func_type="gradient_plane")
z = idw.interp_to_pts(pts)
# z = idw.interp_to_pts(pts, nodal_func_type='gradient_plane')
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[8]:

Notice, in the previous interpolation that we get negative values in some locations. If we assume that there should be no negatives values (like when measuring concentration) we can use the truncation options to adjust the interpolation.

In [9]:
idw.set_truncation_max_min(smax=150, smin=0)
z = idw.interp_to_pts(pts)
# z = idw.interp_to_pts(pts, smax=150, smin=0)
zz = np.reshape(z, (len(y), len(x)))
xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
xa.hvplot(x='x', y='y', height=500, cmap='viridis')
Out[9]:

Interactive interpolation using param, holoviews, and panel

Below we define a param class for the interpolation options. Then we use panel to show widgets and a dynamic plot. The user can change the widget values and then push the "Compute new plot" button to update the plot using the currently specified interpolation options.

In [10]:
class InterpOptions(param.Parameterized):
    interpolation_option = param.ObjectSelector(default='idw', objects=['linear', 'natural_neighbor', 'idw'],
                                                precedence=1)
    nodal_function = param.ObjectSelector(default='constant', objects=['constant', 'gradient_plane', 'quadratic'],
                                          precedence=1.1)
    truncation = param.Boolean(default=True, precedence=2)
    truncation_maximum = param.Number(default=df.loc[df['c'].idxmax()]['c'], precedence=2.1)
    truncation_minimum = param.Number(default=df.loc[df['c'].idxmin()]['c'], precedence=2.2)
    
    @param.depends('interpolation_option', watch=True)
    def _update_nodal_function(self):
        if self.interpolation_option == 'linear':
            self.params('nodal_function').precedence = -1
        else:
            self.params('nodal_function').precedence = 1.1

    @param.depends('truncation', watch=True)
    def _update_truncation_opts(self):
        if self.truncation:
            self.params('truncation_maximum').precedence = 2.1
            self.params('truncation_minimum').precedence = 2.2
        else:
            self.params('truncation_maximum').precedence = -1
            self.params('truncation_minimum').precedence = -1
    
    def generate_zz(self):
        # Get Interp Object
        interp_object = None
        if self.interpolation_option != 'idw':
            #interp_object = linear
            interp_object = interp.interpolate.InterpLinear(pts=df.values)
            if self.interpolation_option == 'natural_neighbor':
                interp_object.interp_to_pt((0,0)) # this will force triangle creation
                interp_object.set_use_natural_neighbor(nodal_func_type=self.nodal_function,
                                                       nd_func_pt_search_opt="nearest_pts")
        else:
            interp_object = interp.interpolate.InterpIdw(pts=df.values)
            interp_object.set_nodal_function(nodal_func_type=self.nodal_function)

        if self.truncation:
            interp_object.set_truncation_max_min(self.truncation_maximum, self.truncation_minimum)
        z = interp_object.interp_to_pts(pts)
        zz = np.reshape(z, (len(y), len(x)))
        xa = xr.DataArray(coords={'x': x, 'y': y}, dims=['x', 'y'], data=zz.T)
        return xa
        
    def get_image(self, *args, **kwargs):
        return self.generate_zz().hvplot(x='x', y='y', height=500, cmap='viridis')
In [11]:
opts = InterpOptions(name="")
dmap = hv.util.Dynamic(opts.get_image(), operation=opts.get_image, streams=[Params(opts)])
pn.Row(opts, dmap)
Out[11]: