.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/pangeo_unstructured_grid.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_pangeo_unstructured_grid.py>`
        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_pangeo_unstructured_grid.py:


************************************
Interpolation of LLC4320 ocean model
************************************

Interpolation of LLC4320 ocean model

The interpolation of this object is based on a :py:class:`R*Tree
<pyinterp.RTree>` structure. To begin with, we start by building this
object. By default, this object considers the WGS-84 geodetic coordinate system.
But you can define another one using the class :py:class:`Spheroid
<pyinterp.geodetic.Spheroid>`.

.. GENERATED FROM PYTHON SOURCE LINES 14-20

.. code-block:: Python

    import cartopy.crs
    import cartopy.mpl.ticker
    import intake
    import matplotlib.pyplot
    import numpy


.. GENERATED FROM PYTHON SOURCE LINES 21-25

.. code-block:: Python

    import pyinterp

    mesh = pyinterp.RTree()


.. GENERATED FROM PYTHON SOURCE LINES 26-33

Then, we will insert points into the tree. The class allows you to add points
using two algorithms. The first one, called :py:meth:`packing
<pyinterp.RTree.packing>`, will enable you to enter the values in the tree at
once. This mechanism is the recommended solution to create an optimized
in-memory structure, both in terms of construction time and queries. When this
is not possible, you can insert new information into the tree as you go along
using the :py:meth:`insert <pyinterp.RTree.insert>` method.

.. GENERATED FROM PYTHON SOURCE LINES 33-37

.. code-block:: Python

    cat_url = 'https://raw.githubusercontent.com/pangeo-data/pangeo-datastore' \
        '/master/intake-catalogs/ocean/llc4320.yaml'
    cat = intake.open_catalog(cat_url)


.. GENERATED FROM PYTHON SOURCE LINES 38-39

Grid subsampling (original volume is too huge for this example)

.. GENERATED FROM PYTHON SOURCE LINES 39-41

.. code-block:: Python

    indices = slice(0, None, 8)


.. GENERATED FROM PYTHON SOURCE LINES 42-43

Reads longitudes and latitudes of the grid

.. GENERATED FROM PYTHON SOURCE LINES 43-47

.. code-block:: Python

    array = cat.LLC4320_grid.to_dask()
    lons = array['XC'].isel(i=indices, j=indices)
    lats = array['YC'].isel(i=indices, j=indices)


.. GENERATED FROM PYTHON SOURCE LINES 48-49

Reads SSH values for the first time step of the time series

.. GENERATED FROM PYTHON SOURCE LINES 49-52

.. code-block:: Python

    ssh = cat.LLC4320_SSH.to_dask()
    ssh = ssh['Eta'].isel(time=0, i=indices, j=indices)


.. GENERATED FROM PYTHON SOURCE LINES 53-54

Populates the search tree

.. GENERATED FROM PYTHON SOURCE LINES 54-58

.. code-block:: Python

    mesh.packing(
        numpy.vstack((lons.values.ravel(), lats.values.ravel())).T,
        ssh.values.ravel())


.. GENERATED FROM PYTHON SOURCE LINES 59-78

When the tree is created, you can interpolate data with two algorithms:

* :py:meth:`Inverse Distance Weighting
  <pyinterp.RTree.inverse_distance_weighting>` or IDW
* :py:meth:`Radial Basis Function
  <pyinterp.RTree.radial_basis_function>` or RBF

Yon can also search the :py:meth:`nearest neighbors <pyinterp.RTree.query>` on
the tree.

.. note::

    When comparing an RBF to IDW, IDW will never predict values higher than
    the maximum measured value or lower than the minimum measured value.
    However, RBFs can predict values higher than the maximum values and lower
    than the minimum measured values.

In this example, we will under-sample the source grid at 1/32 degree over an
area of the globe.

.. GENERATED FROM PYTHON SOURCE LINES 78-85

.. code-block:: Python

    x0, x1 = 80, 170
    y0, y1 = -45, 30
    res = 1 / 32.0
    mx, my = numpy.meshgrid(numpy.arange(x0, x1, res),
                            numpy.arange(y0, y1, res),
                            indexing='ij')


.. GENERATED FROM PYTHON SOURCE LINES 86-87

IDW interpolation

.. GENERATED FROM PYTHON SOURCE LINES 87-95

.. code-block:: Python

    idw_eta, neighbors = mesh.inverse_distance_weighting(
        numpy.vstack((mx.ravel(), my.ravel())).T,
        within=True,  # Extrapolation is forbidden
        radius=55000,  # In a radius of 5.5 Km
        k=8,  # We are looking for at most 8 neighbours
        num_threads=0)
    idw_eta = idw_eta.reshape(mx.shape)


.. GENERATED FROM PYTHON SOURCE LINES 96-97

RBF interpolation

.. GENERATED FROM PYTHON SOURCE LINES 97-104

.. code-block:: Python

    rbf_eta, neighbors = mesh.radial_basis_function(
        numpy.vstack((mx.ravel(), my.ravel())).T,
        within=True,  # Extrapolation is forbidden
        k=11,  # We are looking for at most 11 neighbours
        num_threads=0)
    rbf_eta = rbf_eta.reshape(mx.shape)


.. GENERATED FROM PYTHON SOURCE LINES 105-106

Let's visualize our interpolated data

.. GENERATED FROM PYTHON SOURCE LINES 106-139

.. code-block:: Python

    fig = matplotlib.pyplot.figure(figsize=(18, 9))
    lon_formatter = cartopy.mpl.ticker.LongitudeFormatter(
        zero_direction_label=True)
    lat_formatter = cartopy.mpl.ticker.LatitudeFormatter()
    ax = fig.add_subplot(121, projection=cartopy.crs.PlateCarree())
    ax.pcolormesh(mx,
                  my,
                  idw_eta,
                  cmap='terrain',
                  shading='auto',
                  transform=cartopy.crs.PlateCarree())
    ax.coastlines()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_xticks(numpy.arange(x0, x1, 10.0))
    ax.set_yticks(numpy.arange(y0, y1, 10))
    ax.set_title('Eta (IDW)')

    ax = fig.add_subplot(122, projection=cartopy.crs.PlateCarree())
    ax.pcolormesh(mx,
                  my,
                  rbf_eta,
                  cmap='terrain',
                  shading='auto',
                  transform=cartopy.crs.PlateCarree())
    ax.coastlines()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_xticks(numpy.arange(x0, x1, 10.0))
    ax.set_yticks(numpy.arange(y0, y1, 10))
    ax.set_title('Eta (RBF)')
    fig.show()


.. GENERATED FROM PYTHON SOURCE LINES 140-144

The image below illustrates the result of the IDW interpolation:

.. figure:: ../pictures/mit_gcm.png
    :align: center


.. _sphx_glr_download_auto_examples_pangeo_unstructured_grid.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/pangeo_unstructured_grid.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: pangeo_unstructured_grid.ipynb <pangeo_unstructured_grid.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: pangeo_unstructured_grid.py <pangeo_unstructured_grid.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_