.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/grid/ex_unstructured.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_grid_ex_unstructured.py: Unstructured Grid Interpolation =============================== This example illustrates the various interpolation methods available for unstructured grids, such as data from satellite tracks or other non-uniform sources. The core of unstructured grid interpolation in ``pyinterp`` is the :py:class:`pyinterp.RTree` class, which uses an R-tree data structure for efficient spatial queries. First, we'll create a synthetic dataset to represent our unstructured grid. .. GENERATED FROM PYTHON SOURCE LINES 16-23 .. code-block:: Python import cartopy.crs import cartopy.feature import matplotlib.pyplot import numpy import pyinterp .. GENERATED FROM PYTHON SOURCE LINES 24-26 The R-tree can be configured to work with different geodetic systems. The default is WGS-84. Here, we'll use the default. .. GENERATED FROM PYTHON SOURCE LINES 26-29 .. code-block:: Python mesh = pyinterp.RTree() .. GENERATED FROM PYTHON SOURCE LINES 30-32 We will create a synthetic field of data to simulate measurements from an unstructured grid. .. GENERATED FROM PYTHON SOURCE LINES 32-40 .. code-block:: Python def field(lon, lat): """A synthetic field of data.""" return ( numpy.sin(numpy.radians(lon) * 3) * numpy.cos(numpy.radians(lat) * 2) + 0.5 * numpy.sin(numpy.radians(lon) * 5) * numpy.sin(numpy.radians(lat) * 4)) .. GENERATED FROM PYTHON SOURCE LINES 41-45 Now, we generate random longitude and latitude points and populate our R-tree with these points and their corresponding data values. The :py:meth:`~pyinterp.RTree.packing` method is the most efficient way to build the tree from a complete dataset at once. .. GENERATED FROM PYTHON SOURCE LINES 45-55 .. code-block:: Python N_POINTS = 2000 X0, X1 = 80, 170 Y0, Y1 = -45, 30 generator = numpy.random.Generator(numpy.random.PCG64(0)) lons = generator.uniform(low=X0, high=X1, size=(N_POINTS, )) lats = generator.uniform(low=Y0, high=Y1, size=(N_POINTS, )) data = field(lons, lats) mesh.packing(numpy.vstack((lons, lats)).T, data) .. GENERATED FROM PYTHON SOURCE LINES 56-58 For our interpolation, we need a grid where we want to estimate the values. We'll create a regular grid for this purpose. .. GENERATED FROM PYTHON SOURCE LINES 58-63 .. code-block:: Python STEP = 0.5 mx, my = numpy.meshgrid(numpy.arange(X0, X1 + STEP, STEP), numpy.arange(Y0, Y1 + STEP, STEP), indexing='ij') .. GENERATED FROM PYTHON SOURCE LINES 64-127 Interpolation Methods --------------------- ``pyinterp`` offers several methods for interpolating data from an unstructured grid. We will now apply and visualize each of them. Inverse Distance Weighting (IDW), Radial Basis Function (RBF), and Kriging are all interpolation methods used to estimate a value for a target location based on the values of surrounding sample points. However, each method approaches this estimation differently. IDW uses a weighted average of the surrounding sample points, where the weight assigned to each point is inversely proportional to its distance from the target location. The further away a sample point is from the target location, the less influence it has on the estimated value. This method is relatively simple to implement and computationally efficient, but it can produce over-smoothed results in areas with a lot of sample points and under-smoothed results in areas with few sample points. The Window Function method, much like IDW, uses a weighted average of the surrounding sample points. However, instead of the weight being solely determined by the inverse distance, it is determined by a kernel function (the "window"). This function gives the most weight to points at the center of the window and progressively less weight to points further away. This can sometimes provide a smoother result than IDW. RBF, on the other hand, models the spatial relationship between sample points and the target location by using a mathematical function (radial basis function) that is based on the distance between the points. The radial basis function is usually Gaussian, multiquadric, or inverse multiquadric. The estimated value at the target location is obtained by summing up the weighted contributions of all sample points. This method is more flexible than IDW as it can produce both smooth and non-smooth surfaces depending on the chosen radial basis function. However, it can be computationally expensive and may produce artifacts in areas with sparse or clustered sample points. Kriging is a geostatistical interpolation method that uses a variogram to model the spatial correlation between sample points. The variogram describes how the similarity between sample points decreases as the distance between them increases. The estimated value at the target location is obtained by a weighted average of the surrounding sample points, where the weights are determined by the variogram model. This method is considered to be the most accurate interpolation method as it provides an optimal and unbiased estimate of the value at the target location. However, it is also the most computationally expensive and requires a good understanding of the underlying spatial structure of the data. In summary, the choice of interpolation method is a trade-off between simplicity, computational cost, and the quality of the results. * **IDW** is a good starting point as it is simple and fast. * The **Window Function** can provide a smoother alternative to IDW. * **RBF** offers more flexibility but at a higher computational cost and with a risk of artifacts. * **Kriging** is the most accurate method but also the most complex and computationally intensive, requiring a good understanding of the data's spatial structure. In this notebook, we will compare the results of these four interpolation methods on a synthetic dataset. Inverse Distance Weighting (IDW) ******************************** .. GENERATED FROM PYTHON SOURCE LINES 127-136 .. code-block:: Python idw, _ = mesh.inverse_distance_weighting( numpy.vstack((mx.ravel(), my.ravel())).T, within=False, k=11, radius=None, num_threads=0, ) idw = idw.reshape(mx.shape) .. GENERATED FROM PYTHON SOURCE LINES 137-139 Radial Basis Function (RBF) *************************** .. GENERATED FROM PYTHON SOURCE LINES 139-148 .. code-block:: Python rbf, _ = mesh.radial_basis_function(numpy.vstack((mx.ravel(), my.ravel())).T, within=False, k=11, rbf='multiquadric', epsilon=None, smooth=0, num_threads=0) rbf = rbf.reshape(mx.shape) .. GENERATED FROM PYTHON SOURCE LINES 149-151 Window Function *************** .. GENERATED FROM PYTHON SOURCE LINES 151-159 .. code-block:: Python wf, _ = mesh.window_function(numpy.vstack((mx.ravel(), my.ravel())).T, within=False, k=11, wf='parzen', radius=None, num_threads=0) wf = wf.reshape(mx.shape) .. GENERATED FROM PYTHON SOURCE LINES 160-162 Kriging ******* .. GENERATED FROM PYTHON SOURCE LINES 162-171 .. code-block:: Python kriging, _ = mesh.kriging(numpy.vstack((mx.ravel(), my.ravel())).T, within=False, k=11, radius=None, covariance='gaussian', drift_function='linear', num_threads=0) kriging = kriging.reshape(mx.shape) .. GENERATED FROM PYTHON SOURCE LINES 172-177 Visualization of Results ------------------------ Finally, let's visualize the original scattered data and the results of the different interpolation methods on a map. .. GENERATED FROM PYTHON SOURCE LINES 177-226 .. code-block:: Python fig = matplotlib.pyplot.figure(figsize=(20, 10)) fig.patch.set_alpha(0.0) gs = fig.add_gridspec(2, 4) ax1 = fig.add_subplot(gs[0, 0], projection=cartopy.crs.PlateCarree()) ax2 = fig.add_subplot(gs[0, 1], projection=cartopy.crs.PlateCarree()) ax3 = fig.add_subplot(gs[0, 2], projection=cartopy.crs.PlateCarree()) ax4 = fig.add_subplot(gs[0, 3], projection=cartopy.crs.PlateCarree()) ax5 = fig.add_subplot(gs[1, 0], projection=cartopy.crs.PlateCarree()) ax6 = fig.add_subplot(gs[1, 1], projection=cartopy.crs.PlateCarree()) ax7 = fig.add_subplot(gs[1, 2], projection=cartopy.crs.PlateCarree()) ax8 = fig.add_subplot(gs[1, 3], projection=cartopy.crs.PlateCarree()) # Common plotting function def plot_grid(ax, grid, title, cmap, vmin, vmax): """Helper function to plot interpolated grids.""" pcm = ax.pcolormesh(mx, my, grid, cmap=cmap, transform=cartopy.crs.PlateCarree(), vmin=vmin, vmax=vmax) ax.coastlines() ax.add_feature(cartopy.feature.BORDERS, linestyle='-') ax.set_title(title) return pcm # Plotting each interpolation result true = field(mx, my) plot_grid(ax1, idw, 'Inverse Distance Weighting', 'viridis', -1, 1) plot_grid(ax2, rbf, 'Radial Basis Function', 'viridis', -1, 1) plot_grid(ax3, wf, 'Window Function', 'viridis', -1, 1) plot_grid(ax4, kriging, 'Kriging', 'viridis', -1, 1) # Plotting errors error_idw = numpy.where(true != 0, (idw - true) / true, 0) error_rbf = numpy.where(true != 0, (rbf - true) / true, 0) error_wf = numpy.where(true != 0, (wf - true) / true, 0) error_kriging = numpy.where(true != 0, (kriging - true) / true, 0) plot_grid(ax5, error_idw, 'IDW Relative Error', 'RdBu_r', -0.5, 0.5) plot_grid(ax6, error_rbf, 'RBF Relative Error', 'RdBu_r', -0.5, 0.5) plot_grid(ax7, error_wf, 'Window Function Relative Error', 'RdBu_r', -0.5, 0.5) plot_grid(ax8, error_kriging, 'Kriging Relative Error', 'RdBu_r', -0.5, 0.5) fig.tight_layout() matplotlib.pyplot.show() .. image-sg:: /auto_examples/grid/images/sphx_glr_ex_unstructured_001.png :alt: Inverse Distance Weighting, Radial Basis Function, Window Function, Kriging, IDW Relative Error, RBF Relative Error, Window Function Relative Error, Kriging Relative Error :srcset: /auto_examples/grid/images/sphx_glr_ex_unstructured_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.525 seconds) .. _sphx_glr_download_auto_examples_grid_ex_unstructured.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/CNES/pangeo-pyinterp/master?urlpath=lab/tree/notebooks/auto_examples/grid/ex_unstructured.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ex_unstructured.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ex_unstructured.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ex_unstructured.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_