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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_examples_plot_shoot_netcdf.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_plot_shoot_netcdf.py:


Open and analyse shoot output netcdf
====================================

.. GENERATED FROM PYTHON SOURCE LINES 6-19

.. code-block:: Python


    import cartopy.crs as ccrs
    import cmocean as cm
    import matplotlib.pyplot as plt
    import numpy as np
    import xarray as xr

    from shoot.num import points_in_polygon
    from shoot.plot import create_map, pcarr
    from shoot.samples import get_sample_file

    pmerc = ccrs.Mercator()








.. GENERATED FROM PYTHON SOURCE LINES 20-22

Usefull function
----------------

.. GENERATED FROM PYTHON SOURCE LINES 22-72

.. code-block:: Python



    def create_map_ax(
        lons,
        lats,
        ax,
        margin=0.0,
        square=False,
        coastlines=True,
        emodnet=False,
        title=None,
        **kwargs,
    ):
        """Create a simple decorated cartopy map"""
        # lons, lats = lons.values, lats.values
        xmin, xmax = np.min(lons), np.max(lons)
        ymin, ymax = np.min(lats), np.max(lats)
        dx, dy = xmax - xmin, ymax - ymin
        x0, y0 = 0.5 * (xmin + xmax), 0.5 * (ymin + ymax)
        if square:
            aspect = dx / dy * np.cos(np.radians(y0))
            if aspect > 1:
                dy *= aspect
            else:
                dx /= aspect
        xmargin = margin * dx
        ymargin = margin * dy
        xmin = x0 - 0.5 * dx - xmargin
        xmax = x0 + 0.5 * dx + xmargin
        ymin = y0 - 0.5 * dy - ymargin
        ymax = y0 + 0.5 * dy + ymargin

        ax.set_extent([xmin, xmax, ymin, ymax])
        ax.gridlines(
            draw_labels=["bottom", "left"],
            linewidth=1,
            color="k",
            alpha=0.5,
            linestyle="--",
            rotate_labels=False,
        )
        if coastlines:
            ax.coastlines()
        if emodnet:
            ax.add_wms("https://ows.emodnet-bathymetry.eu/wms", "emodnet:mean_atlas_land")
        if title:
            ax.set_title(title)
        return fig, ax









.. GENERATED FROM PYTHON SOURCE LINES 73-75

Initialisation
--------------

.. GENERATED FROM PYTHON SOURCE LINES 75-87

.. code-block:: Python


    #  Import data

    root_path = "OBS/SATELLITE/jan2024_ionian_sea_duacs.nc"
    path = get_sample_file(root_path)
    ds = xr.open_dataset(path)
    root_track = "EDDIES/track_ionian_sea_duacs_jan2024.nc"
    track_path = get_sample_file(root_track)
    tracks = xr.open_dataset(track_path)
    tracks







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <pre>&lt;xarray.Dataset&gt; Size: 1MB
    Dimensions:               (obs: 820, nb_sample: 50, eddies: 71)
    Dimensions without coordinates: obs, nb_sample, eddies
    Data variables: (12/30)
        time                  (obs) datetime64[ns] 7kB ...
        i_cen                 (obs) int64 7kB ...
        j_cen                 (obs) int64 7kB ...
        x_cen                 (obs) float64 7kB ...
        y_cen                 (obs) float64 7kB ...
        track_id              (obs) int64 7kB ...
        ...                    ...
        life_time             (eddies) float64 568B ...
        x_start               (eddies) float64 568B ...
        y_start               (eddies) float64 568B ...
        x_end                 (eddies) float64 568B ...
        y_end                 (eddies) float64 568B ...
        eddy_type             (eddies) &lt;U11 3kB ...
    Attributes:
        window_center:  50 km
        window_fit:     120 km
        min_radius:     20 km
        project:        SHOOT
        institution:    SHOM
        contact:        jean.baptiste.roustan@shom.fr</pre>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 88-90

Plots
-----

.. GENERATED FROM PYTHON SOURCE LINES 90-150

.. code-block:: Python


    # Plot tracking

    fig, ax = create_map(ds.longitude, ds.latitude, figsize=(8, 5))
    colors = {"cyclone": "b", "anticyclone": "r"}

    dss = ds.isel(time=-1)

    cb = dss.adt.plot(
        x="longitude",
        y="latitude",
        cmap="cmo.dense",
        ax=ax,
        add_colorbar=False,
        transform=pcarr,
    )

    plt.colorbar(cb, label="ADT [m]")

    nj = 1
    plt.quiver(
        dss.longitude[::nj].values,
        dss.latitude[::nj].values,
        dss.ugos[::nj, ::nj].values,
        dss.vgos[::nj, ::nj].values,
        transform=pcarr,
    )


    day = np.where(tracks.time.values == dss.time.values)[0]
    track_day = tracks.isel(obs=day)

    for i in range(len(track_day.obs)):
        tmp = track_day.isel(obs=i)
        plt.scatter(
            tmp.x_cen,
            tmp.y_cen,
            c=colors[str(tmp.eddy_type.isel(eddies=tmp.track_id.values).values)],
            s=10,
            transform=pcarr,
        )

        plt.plot(
            tmp.x_vmax_contour,
            tmp.y_vmax_contour,
            "--",
            c=colors[str(tmp.eddy_type.isel(eddies=tmp.track_id.values).values)],
            transform=pcarr,
        )

        plt.plot(
            tmp.x_eff_contour,
            tmp.y_eff_contour,
            c=colors[str(tmp.eddy_type.isel(eddies=tmp.track_id.values).values)],
            transform=pcarr,
        )

    plt.title(np.datetime_as_string(dss.time, unit="D"))
    plt.tight_layout()




.. image-sg:: /examples/images/sphx_glr_plot_shoot_netcdf_001.png
   :alt: 2024-01-31
   :srcset: /examples/images/sphx_glr_plot_shoot_netcdf_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 151-152

#### Heatmaps

.. GENERATED FROM PYTHON SOURCE LINES 152-230

.. code-block:: Python


    ## to be modified by the user
    lons = np.arange(-6, 36, 0.1)
    lats = np.arange(30, 44.5, 0.1)

    count_s_anti = np.zeros((len(lats), len(lons)))
    count_s_cyc = np.zeros((len(lats), len(lons)))

    ind_anti = np.where(tracks.eddy_type[tracks.track_id.values].values == "anticyclone")[0]
    ind_cyc = np.where(tracks.eddy_type[tracks.track_id.values].values == "cyclone")[0]

    shoot_anti = tracks.isel(obs=ind_anti)
    shoot_cyc = tracks.isel(obs=ind_cyc)

    ## Should work with previously made modifications

    Xlons, Xlats = np.meshgrid(lons, lats)


    for i in range(len(shoot_anti.obs)):
        line = np.array(
            [
                shoot_anti.isel(obs=i).x_vmax_contour,
                shoot_anti.isel(obs=i).y_vmax_contour,
            ]
        ).T
        for j in range(len(lons)):  # On parcoure le tableau via les longitude
            points = np.array([Xlons[:, j], Xlats[:, j]]).T
            in_poly = points_in_polygon(points, line)
            # print(in_poly)
            count_s_anti[:, j] += in_poly

    for i in range(len(shoot_cyc.obs)):
        line = np.array(
            [
                shoot_cyc.isel(obs=i).x_vmax_contour,
                shoot_cyc.isel(obs=i).y_vmax_contour,
            ]
        ).T
        for j in range(len(lons)):  # On parcoure le tableau via les longitude
            points = np.array([Xlons[:, j], Xlats[:, j]]).T
            in_poly = points_in_polygon(points, line)
            # print(in_poly)
            count_s_cyc[:, j] += in_poly


    count_s_anti[count_s_anti == 0] = np.nan
    count_s_cyc[count_s_cyc == 0] = np.nan


    fig, axs = plt.subplots(1, 2, subplot_kw=dict(projection=pmerc), figsize=(10, 5))
    # create_map_ax(np.arange(-6, 36), np.arange(30, 44.5), axs[0])
    create_map_ax(np.arange(15, 30), np.arange(30, 40), axs[0])
    plt.sca(axs[0])
    plt.title("Anticyclone occurence 2023-2024 (MED1.8)")
    plt.pcolormesh(
        lons,
        lats,
        count_s_anti / 510,
        transform=pcarr,
        vmin=0,
        vmax=0.7,
        cmap=cm.cm.thermal,
    )

    # create_map_ax(np.arange(-6, 36), np.arange(30, 44.5), axs[1])
    create_map_ax(np.arange(15, 30), np.arange(30, 40), axs[1])
    plt.sca(axs[1])
    plt.title("Cyclone occurence 2023-2024 (MED1.8)")
    plt.pcolormesh(
        lons,
        lats,
        count_s_cyc / 510,
        transform=pcarr,
        vmin=0,
        vmax=0.7,
        cmap=cm.cm.thermal,
    )



.. image-sg:: /examples/images/sphx_glr_plot_shoot_netcdf_002.png
   :alt: Anticyclone occurence 2023-2024 (MED1.8), Cyclone occurence 2023-2024 (MED1.8)
   :srcset: /examples/images/sphx_glr_plot_shoot_netcdf_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <cartopy.mpl.geocollection.GeoQuadMesh object at 0x70269b43de90>




.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 8.434 seconds)


.. _sphx_glr_download_examples_plot_shoot_netcdf.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_shoot_netcdf.zip <plot_shoot_netcdf.zip>`


.. only:: html

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

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