.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/grid/ex_xarray_interpolator.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_xarray_interpolator.py: .. _example_xarray_interpolator: Interpolating an xarray DataArray ================================= :class:`pyinterp.backends.xarray.RegularGridInterpolator` is the one-call entry point to interpolate a 2-, 3- or 4-D ``xarray.DataArray``: it identifies the axes, builds the matching grid and dispatches to the bivariate, trivariate or quadrivariate interpolator. In this example you will learn how to build it, query it with a dictionary of coordinates (including dates), choose and tune the interpolation method, and what happens at the dateline and outside the grid. .. GENERATED FROM PYTHON SOURCE LINES 14-23 .. code-block:: Python import matplotlib.pyplot import numpy import xarray import pyinterp import pyinterp.backends.xarray import pyinterp.tests .. GENERATED FROM PYTHON SOURCE LINES 24-34 Building the Interpolator ------------------------- We load a global grid of mean sea surface (MSS). The variable is stored as ``(lat, lon)``, but the interpolator does not care about the order of the dimensions: the longitude and latitude are recognized from the ``units`` attribute of their coordinates (``degrees_east`` and ``degrees_north``, see :class:`pyinterp.cf.AxisLongitudeUnit` and :class:`pyinterp.cf.AxisLatitudeUnit`), and the data is transposed to ``(longitude, latitude)`` internally. .. GENERATED FROM PYTHON SOURCE LINES 34-39 .. code-block:: Python ds = pyinterp.tests.load_grid2d() interpolator = pyinterp.backends.xarray.RegularGridInterpolator(ds.mss) print(f"dims={ds.mss.dims}, ndim={interpolator.ndim}") print(interpolator.grid) .. rst-class:: sphx-glr-script-out .. code-block:: none dims=('lat', 'lon'), ndim=2 Grid2D(shape=(1440, 712), dtype=float64, nbytes=7 MB) .. GENERATED FROM PYTHON SOURCE LINES 40-48 Querying with a Dictionary of Coordinates ----------------------------------------- The interpolator is called with a dictionary that maps every dimension name to the coordinates of the points to evaluate. The values can be scalars, lists or arrays of any shape: they are broadcast against each other, and the result has their broadcast shape, here the shape of the 2-D target grid. The MSS is stored with longitudes in [0, 360), but we query it in [-180, 180). .. GENERATED FROM PYTHON SOURCE LINES 48-56 .. code-block:: Python mx, my = numpy.meshgrid( numpy.arange(-180, 180, 1.0) + 1 / 3, numpy.arange(-80, 80, 1.0) + 1 / 3, indexing="ij", ) coords = {"lon": mx, "lat": my} bilinear = interpolator(coords) .. GENERATED FROM PYTHON SOURCE LINES 57-69 Choosing the Method ------------------- The ``method`` argument selects the interpolator. ``bilinear`` (default), ``nearest`` and ``idw`` only use the cell surrounding the point. Windowed methods (``bicubic``, ``c_spline``, ``akima``, ``steffen``...) fit a function on a window of neighboring nodes; its size is set by ``half_window_size_x`` and ``half_window_size_y``, and ``boundary_mode`` decides whether the window may shrink near the edges of the grid (``"shrink"``) or whether the result is undefined there (``"undef"``, the default). The other keyword arguments, such as ``num_threads``, are forwarded to :func:`pyinterp.bivariate`. .. GENERATED FROM PYTHON SOURCE LINES 69-77 .. code-block:: Python bicubic = interpolator( coords, method="bicubic", half_window_size_x=3, half_window_size_y=3, boundary_mode="shrink", ) .. GENERATED FROM PYTHON SOURCE LINES 78-82 The two methods agree to within a millimeter; the largest differences follow the strongest gradients of the field. Note that a window containing a single masked (NaN) node gives NaN, so the windowed method loses a few points along the coasts. .. GENERATED FROM PYTHON SOURCE LINES 82-101 .. code-block:: Python print( f"NaN: bilinear={numpy.isnan(bilinear).sum()}, " f"bicubic={numpy.isnan(bicubic).sum()}" ) fig, (ax1, ax2) = matplotlib.pyplot.subplots( 2, 1, figsize=(10, 6), sharex=True, layout="constrained" ) pcm = ax1.pcolormesh(mx, my, bilinear, cmap="viridis", vmin=-0.1, vmax=0.1) fig.colorbar(pcm, ax=ax1, label="MSS (m)") ax1.set_title("Bilinear interpolation, longitudes in [-180, 180)") pcm = ax2.pcolormesh( mx, my, (bicubic - bilinear) * 1e3, cmap="RdBu_r", vmin=-0.5, vmax=0.5 ) fig.colorbar(pcm, ax=ax2, label="mm") ax2.set_title("Bicubic minus bilinear") ax2.set_xlabel("Longitude (degrees)") for ax in (ax1, ax2): ax.set_ylabel("Latitude (degrees)") .. image-sg:: /auto_examples/grid/images/sphx_glr_ex_xarray_interpolator_001.png :alt: Bilinear interpolation, longitudes in [-180, 180), Bicubic minus bilinear :srcset: /auto_examples/grid/images/sphx_glr_ex_xarray_interpolator_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none NaN: bilinear=19133, bicubic=21880 .. GENERATED FROM PYTHON SOURCE LINES 102-109 Longitude Periodicity --------------------- A longitude axis is always created as periodic (360 degrees). Longitudes are therefore never out of bounds: the same location written as -0.1, 359.9 or 719.9 degrees gives the same value, and points between the last column (359.75) and 360 are interpolated across the seam of the grid. .. GENERATED FROM PYTHON SOURCE LINES 109-111 .. code-block:: Python print(interpolator({"lon": [-0.1, 359.9, 719.9], "lat": -40.1})) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.017 0.017 0.017] .. GENERATED FROM PYTHON SOURCE LINES 112-117 Outside the Grid ---------------- The latitudes cover [-89, 88.75]. A point beyond this range returns NaN by default; set ``bounds_error=True`` to raise a :class:`ValueError` instead. .. GENERATED FROM PYTHON SOURCE LINES 117-124 .. code-block:: Python outside = {"lon": -30.0, "lat": 89.5} print(interpolator(outside)) try: interpolator(outside, bounds_error=True) except ValueError as exc: print(f"ValueError: {exc}") .. rst-class:: sphx-glr-script-out .. code-block:: none nan ValueError: 89.5 is out of bounds for axis y [-89, ..., 88.75] .. GENERATED FROM PYTHON SOURCE LINES 125-133 Adding a Time Dimension ----------------------- With a 3-D array the interpolator becomes trivariate. Here the total column water vapor is stored as ``(time, latitude, longitude)`` with daily values. A coordinate of type ``datetime64`` is detected as a temporal axis and moved to the third position, so the query dates are simply passed as ``datetime64`` values, at any resolution. .. GENERATED FROM PYTHON SOURCE LINES 133-148 .. code-block:: Python tcw = pyinterp.tests.load_grid3d().tcw interpolator = pyinterp.backends.xarray.RegularGridInterpolator(tcw) print(interpolator.grid) dates = numpy.arange( numpy.datetime64("2002-07-02T00", "h"), numpy.datetime64("2002-07-12T00", "h"), numpy.timedelta64(3, "h"), ) track = { "longitude": 10.3, "latitude": 45.1, "time": dates, } .. rst-class:: sphx-glr-script-out .. code-block:: none TemporalGrid3D(shape=(144, 73, 31), dtype=float64, nbytes=2 MB) .. GENERATED FROM PYTHON SOURCE LINES 149-154 The default ``bilinear`` method is linear along the time axis. The ``third_axis`` option (``"linear"`` or ``"nearest"``) selects the interpolation along time; for a 4-D array, ``fourth_axis`` plays the same role for the fourth dimension. These options, like ``half_window_size_x``, select the windowed version of ``bilinear``. .. GENERATED FROM PYTHON SOURCE LINES 154-164 .. code-block:: Python linear = interpolator(track) nearest = interpolator(track, method="bilinear", third_axis="nearest") fig, ax = matplotlib.pyplot.subplots(figsize=(10, 4), layout="constrained") ax.plot(dates, linear, label="bilinear (linear in time)") ax.plot(dates, nearest, label="bilinear, third_axis='nearest'") ax.set_ylabel("Total column water vapor (kg m$^{-2}$)") ax.set_title("Time series at 10.3E, 45.1N") ax.legend() .. image-sg:: /auto_examples/grid/images/sphx_glr_ex_xarray_interpolator_002.png :alt: Time series at 10.3E, 45.1N :srcset: /auto_examples/grid/images/sphx_glr_ex_xarray_interpolator_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 165-168 A 4-D array works the same way: with ``(time, level, latitude, longitude)`` dimensions, the axes are reordered to ``(longitude, latitude, time, level)``, and the dictionary must name all four dimensions. .. GENERATED FROM PYTHON SOURCE LINES 168-172 .. code-block:: Python pressure = pyinterp.tests.load_grid4d().pressure interpolator = pyinterp.backends.xarray.RegularGridInterpolator(pressure) print(interpolator.grid) .. rst-class:: sphx-glr-script-out .. code-block:: none TemporalGrid4D(shape=(12, 6, 2, 2), dtype=float32, nbytes=1 KB) .. GENERATED FROM PYTHON SOURCE LINES 173-183 Labelled Results with the ``pyinterp`` Accessor ----------------------------------------------- Once pyinterp and xarray are imported, data arrays and datasets have a ``pyinterp`` accessor. Its :meth:`~pyinterp.backends.xarray.DataArrayAccessor.interp` method follows :meth:`xarray.DataArray.interp`: the coordinates are passed as keyword arguments, and the result is a ``DataArray`` labelled with them, which keeps the name and the attributes of the variable. One-dimensional coordinates define a new grid, so there is no ``meshgrid`` to build. .. GENERATED FROM PYTHON SOURCE LINES 183-190 .. code-block:: Python mss = ds.mss.pyinterp.interp( lon=numpy.arange(-180, 180, 0.5), lat=numpy.arange(-80, 80, 0.5), method="bicubic", ) print(mss) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 2MB array([[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ..., [0.002, 0.002, 0.002, ..., 0.001, 0.002, 0.002], [0.002, 0.002, 0.002, ..., 0.002, 0.002, 0.002], [0.002, 0.002, 0.002, ..., 0.002, 0.002, 0.002]], shape=(320, 720)) Coordinates: * lat (lat) float64 3kB -80.0 -79.5 -79.0 -78.5 ... 78.0 78.5 79.0 79.5 * lon (lon) float64 6kB -180.0 -179.5 -179.0 -178.5 ... 178.5 179.0 179.5 Attributes: long_name: Mean sea surface height above reference ellipsoid units: m grid_mapping: crs source: SMO15_1s60_V2.dat cell_methods: lon: point lat: point .. GENERATED FROM PYTHON SOURCE LINES 191-195 Coordinates given as data arrays sharing a dimension are interpolated point by point, along a track for instance. The keyword arguments that do not name a dimension are the options of the interpolator, and dates can be given as strings. .. GENERATED FROM PYTHON SOURCE LINES 195-203 .. code-block:: Python track = tcw.pyinterp.interp( longitude=xarray.DataArray(numpy.linspace(-30, 30, 50), dims="obs"), latitude=xarray.DataArray(numpy.linspace(-20, 40, 50), dims="obs"), time="2002-07-05T12:00", bounds_error=True, ) print(track) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 400B array([24.47455481, 27.78406211, 30.50603547, 25.12550513, 19.3178315 , 25.70561549, 32.64532029, 35.09478524, 36.64133985, 37.316209 , 37.92638815, 37.06678047, 35.08439326, 33.29223999, 31.53000256, 31.32761604, 31.54333934, 34.86044902, 39.56814296, 43.9054787 , 49.14237962, 52.02903909, 52.62723436, 52.20315115, 50.58550299, 49.70437078, 49.69870273, 48.10531923, 45.46202381, 42.17753016, 38.70367639, 32.1830379 , 20.17924025, 14.31828244, 12.80989271, 11.3287591 , 14.6402907 , 18.40654166, 20.3557168 , 21.63105744, 23.3943745 , 23.74053944, 23.09037355, 22.39021747, 22.59515799, 22.1344152 , 20.58533696, 18.52913462, 17.40691458, 18.09264074]) Coordinates: latitude (obs) float64 400B -20.0 -18.78 -17.55 ... 37.55 38.78 40.0 longitude (obs) float64 400B -30.0 -28.78 -27.55 ... 27.55 28.78 30.0 time datetime64[ns] 8B 2002-07-05T12:00:00 Dimensions without coordinates: obs Attributes: units: kg m**-2 long_name: Total column water .. GENERATED FROM PYTHON SOURCE LINES 204-223 When to Use It -------------- The ``pyinterp`` accessor is the shortest path when the data is already an ``xarray.DataArray`` and a labelled result is wanted; it builds the grid at each call, so keep a :class:`~pyinterp.backends.xarray.RegularGridInterpolator` to interpolate the same variable many times. The interpolator is one class for any number of dimensions, coordinates passed by name, dates handled for you. It is a thin wrapper around :class:`~pyinterp.backends.xarray.Grid2D`, :class:`~pyinterp.backends.xarray.Grid3D` and :class:`~pyinterp.backends.xarray.Grid4D`; use these classes and their ``bivariate``, ``trivariate`` or ``quadrivariate`` methods, or build a :class:`pyinterp.Grid` and call :func:`pyinterp.bivariate` and friends, when you work with plain ``numpy`` arrays, need full control over the axes (a non-periodic longitude, for instance), or reuse a configuration object. These lower-level interfaces are shown in :ref:`example_2d_interpolation`, :ref:`example_3d_interpolation` and :ref:`example_4d_interpolation`. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.163 seconds) .. _sphx_glr_download_auto_examples_grid_ex_xarray_interpolator.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_xarray_interpolator.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ex_xarray_interpolator.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ex_xarray_interpolator.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ex_xarray_interpolator.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_