.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/plot_hycom_gdp.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_examples_plot_hycom_gdp.py: Compare Hycom3d with a GDP drifter ================================== In this notebook, we show : * how to decode dataset so that it is easy to access generic coordinates and variables, * how to compute depths from layer thicknesses, * how to interpolate currents from U and V positions to T position on an arakawa C grid, * how to perform a 4D interpolation with a variable depth coordinate to random positions, * how to make an horizontal slice of a 4D variable with a variable depth coordinate. .. GENERATED FROM PYTHON SOURCE LINES 17-19 Initialisations --------------- .. GENERATED FROM PYTHON SOURCE LINES 19-34 .. code-block:: Python import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs import xoa from xoa.grid import dz2depth, shift from xoa.regrid import grid2loc, regrid1d import xoa.meta as xmeta import xoa.geo as xgeo from xoa.plot import plot_flow, plot_double_minimap xr.set_options(display_style="text") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 35-36 Register the :ref:`xoa ` accessors: .. GENERATED FROM PYTHON SOURCE LINES 36-39 .. code-block:: Python xoa.register_accessors() .. GENERATED FROM PYTHON SOURCE LINES 40-41 Set the internal Hycom naming specifications as the current ones .. GENERATED FROM PYTHON SOURCE LINES 41-45 .. code-block:: Python hycom_cfg_file = xoa.get_meta_config_file("hycom") print(hycom_cfg_file) xmeta.set_meta_specs(hycom_cfg_file) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/xoa/conda/develop/lib/python3.14/site-packages/xoa/meta/configs/hycom.cfg .. GENERATED FROM PYTHON SOURCE LINES 46-51 Note that, since this file is internal, you can do it simply with:: xmeta.set_meta_specs("hycom") Here is what these CF specifications contain .. GENERATED FROM PYTHON SOURCE LINES 51-55 .. code-block:: Python with open(hycom_cfg_file) as f: print(f.read()) .. rst-class:: sphx-glr-script-out .. code-block:: none [register] name=hycom [[attrs]] nsigma = "[0-9]?" [sglocator] name_format="{root}_{loc}" valid_locations=u,v [vertical] positive=down type=dz [data_vars] [[u]] loc=u [[[add_coords_loc]]] lon=True lat=True x=True y=True [[v]] loc=v [[[add_coords_loc]]] lon=True lat=True x=True y=True [[bathy]] name=bathymetry add_loc=True [[dz]] name=h [coords] [[x]] name=X [[y]] name=Y [[z]] name=lev .. GENERATED FROM PYTHON SOURCE LINES 56-58 Note that most of these specifications are not needed to decode the dataset, i.e to find known variables and coordinates, since the default specifications help for that. .. GENERATED FROM PYTHON SOURCE LINES 60-73 Read and decode the datasets ---------------------------- Hycom velocities ~~~~~~~~~~~~~~~~ U and V are stored in separate files at U and V locations on the Hycom Arakawa C-Grid. We choose to rename the dataset contents so that we can use generic names to access it. We call for that the :meth:`~xarray.Dataset.xoa.decode` accessor method. Note that the staggered grid location is automatically appended to horizontal coordinates and dimensions of the U and V velocity components since they are placed at U and V locations. This prevent any conflict with the bathymetry and layer thickness horizontal coordinates and dimensions, which point to the T location. .. GENERATED FROM PYTHON SOURCE LINES 75-76 U velocity component .. GENERATED FROM PYTHON SOURCE LINES 76-79 .. code-block:: Python u_hycom = xoa.open_data_sample("MODELS/HYCOM/hycom.gdp.u.nc").xoa.decode().u.squeeze(drop=True) print(u_hycom) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 172kB [42900 values with dtype=float32] Coordinates: * time (time) datetime64[ns] 400B 2021-03-01 ... 2021-03-25T12:00:00 * z (z) int32 24B 1 2 3 4 5 6 lon_u (y_u, x_u) float32 572B ... lat_u (y_u, x_u) float32 572B ... Dimensions without coordinates: y_u, x_u Attributes: long_name: Vitesse u units: m/s standard_name: sea_water_x_velocity .. GENERATED FROM PYTHON SOURCE LINES 80-81 V velocity component .. GENERATED FROM PYTHON SOURCE LINES 81-83 .. code-block:: Python v_hycom = xoa.open_data_sample("MODELS/HYCOM/hycom.gdp.v.nc").xoa.decode().v.squeeze(drop=True) .. GENERATED FROM PYTHON SOURCE LINES 84-85 Layer thicknesses dataset .. GENERATED FROM PYTHON SOURCE LINES 85-88 .. code-block:: Python hycom = xoa.open_data_sample("MODELS/HYCOM/hycom.gdp.h.nc").xoa.decode().squeeze(drop=True) print(hycom) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 174kB Dimensions: (y: 11, x: 13, time: 50, z: 6) Coordinates: lon (y, x) float32 572B ... lat (y, x) float32 572B ... * time (time) datetime64[ns] 400B 2021-03-01 ... 2021-03-25T12:00:00 * z (z) int32 24B 1 2 3 4 5 6 Dimensions without coordinates: y, x Data variables: bathy (y, x) float32 572B ... dz (time, z, y, x) float32 172kB ... Attributes: (12/26) history: Fri Apr 9 12:15:22 2021: ncrcat -d X,472,508,... Conventions: COARDS nhybrd: 20 nsigma: 19 kapref: 2 hybflg: 0 ... ... dp00ft: [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ... dp00xt: [21.874 21.414 21.414 13.197 7.85 3.85 4.... ds00t: [3.555 3.026 3.026 1.547 0.828 0.828 0.968 0.8... ds00ft: [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ... ds00xt: [21.874 21.414 21.414 7.197 3.85 3.85 4.... nco_openmp_thread_number: 1 .. GENERATED FROM PYTHON SOURCE LINES 89-92 We now compute the depths from the layer thicknesses at the T location thanks to the :func:`xoa.grid.dz2depth` function: we integrate from a null SSH and interpolate the depth from W locations to T locations. .. GENERATED FROM PYTHON SOURCE LINES 92-95 .. code-block:: Python hycom.coords["depth"] = dz2depth(hycom.dz, centered=True) .. GENERATED FROM PYTHON SOURCE LINES 96-100 Then we interpolate the velocity components to the T location with the :func:`xoa.grid.shift` function. The call to :meth:`reloc ` helps removing the staggered grid location prefixes. .. GENERATED FROM PYTHON SOURCE LINES 100-106 .. code-block:: Python ut_hycom = shift(u_hycom, {"x": "left", "y": "left"}).xoa.reloc(u=False) vt_hycom = shift(v_hycom, {"x": "left", "y": "left"}).xoa.reloc(v=False) hycom["u"] = ut_hycom.assign_coords(**hycom.coords) hycom["v"] = vt_hycom.assign_coords(**hycom.coords) .. GENERATED FROM PYTHON SOURCE LINES 107-108 So, finally we obtain: .. GENERATED FROM PYTHON SOURCE LINES 108-111 .. code-block:: Python print(hycom) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 689kB Dimensions: (y: 11, x: 13, time: 50, z: 6) Coordinates: lon (y, x) float32 572B ... lat (y, x) float32 572B ... depth (time, z, y, x) float32 172kB 1.511 1.511 1.51 ... 20.09 20.06 * time (time) datetime64[ns] 400B 2021-03-01 ... 2021-03-25T12:00:00 * z (z) int32 24B 1 2 3 4 5 6 Dimensions without coordinates: y, x Data variables: bathy (y, x) float32 572B ... dz (time, z, y, x) float32 172kB 3.021 3.021 3.02 ... 5.22 5.216 5.209 u (time, z, y, x) float32 172kB 0.03078 0.02984 ... 0.1292 0.1364 v (time, z, y, x) float32 172kB 0.0497 0.07527 ... 0.0498 0.05054 Attributes: (12/26) history: Fri Apr 9 12:15:22 2021: ncrcat -d X,472,508,... Conventions: COARDS nhybrd: 20 nsigma: 19 kapref: 2 hybflg: 0 ... ... dp00ft: [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ... dp00xt: [21.874 21.414 21.414 13.197 7.85 3.85 4.... ds00t: [3.555 3.026 3.026 1.547 0.828 0.828 0.968 0.8... ds00ft: [1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15 ... ds00xt: [21.874 21.414 21.414 7.197 3.85 3.85 4.... nco_openmp_thread_number: 1 .. GENERATED FROM PYTHON SOURCE LINES 112-116 GDP drifter ~~~~~~~~~~~ The drifter comes as a `csv` file and we read it as :class:`pandas.DataFrame` instance. .. GENERATED FROM PYTHON SOURCE LINES 116-122 .. code-block:: Python csv_name = xoa.get_data_sample("OBS/DRIFTERS/gdp-6203641.csv") drifter = pd.read_csv( csv_name, header=0, skiprows=[1], parse_dates=[0], index_col=0, usecols=[2, 3, 4, 5, 6, 7] ) .. GENERATED FROM PYTHON SOURCE LINES 123-124 Since the sampling is not that nice, we resample it to 3-hour intervals. .. GENERATED FROM PYTHON SOURCE LINES 124-127 .. code-block:: Python drifter = drifter.resample("3h").mean() .. GENERATED FROM PYTHON SOURCE LINES 128-129 We switch to default naming conventions for the drifter. .. GENERATED FROM PYTHON SOURCE LINES 129-133 .. code-block:: Python xmeta.set_meta_specs("default") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 134-135 We convert it to and :class:`xarray.Dataset`, fix the time coordinate and decode it. .. GENERATED FROM PYTHON SOURCE LINES 135-138 .. code-block:: Python drifter = drifter.to_xarray().assign_coords(time=drifter.index.values).xoa.decode() .. GENERATED FROM PYTHON SOURCE LINES 139-140 We drop missing values. .. GENERATED FROM PYTHON SOURCE LINES 140-143 .. code-block:: Python drifter = drifter.where(~drifter.lon.isnull() & ~drifter.lat.isnull(), drop=True) .. GENERATED FROM PYTHON SOURCE LINES 144-145 We add a constant depth of 15 m. .. GENERATED FROM PYTHON SOURCE LINES 145-148 .. code-block:: Python drifter.coords["depth"] = drifter.lon * 0 + 15 .. GENERATED FROM PYTHON SOURCE LINES 149-150 Here is what we obtain. .. GENERATED FROM PYTHON SOURCE LINES 150-153 .. code-block:: Python print(drifter) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 9kB Dimensions: (time: 159) Coordinates: * time (time) datetime64[us] 1kB 2021-03-03T18:00:00 ... 2021-03-24T09:... lat (time) float64 1kB 44.68 44.72 44.72 44.72 ... 44.72 44.71 44.71 lon (time) float64 1kB -2.544 -2.583 -2.592 ... -2.523 -2.514 -2.505 depth (time) float64 1kB 15.0 15.0 15.0 15.0 15.0 ... 15.0 15.0 15.0 15.0 Data variables: sst (time) float64 1kB 13.09 12.79 12.77 12.78 ... 12.56 12.53 12.61 slp (time) float64 1kB 1.027e+03 1.024e+03 ... 1.025e+03 1.026e+03 lon360 (time) float64 1kB 357.5 357.4 357.4 357.4 ... 357.5 357.5 357.5 .. GENERATED FROM PYTHON SOURCE LINES 154-158 Compute and interpolate velocities ---------------------------------- We compute the drifter velocity components .. GENERATED FROM PYTHON SOURCE LINES 158-167 .. code-block:: Python drifter["u"] = ( drifter.lon.differentiate("time", datetime_unit="s") * xgeo.EARTH_RADIUS * np.pi / 180 ) drifter["u"] *= np.cos(np.radians(drifter.lat.values)) drifter["v"] = ( drifter.lat.differentiate("time", datetime_unit="s") * xgeo.EARTH_RADIUS * np.pi / 180 ) .. GENERATED FROM PYTHON SOURCE LINES 168-169 We make a 4D linear interpolation the Hycom velocity to the drifter positions with :func:`xoa.regrid.grid2loc`. .. GENERATED FROM PYTHON SOURCE LINES 169-173 .. code-block:: Python uloc = grid2loc(hycom["u"], drifter) vloc = grid2loc(hycom["v"], drifter) .. GENERATED FROM PYTHON SOURCE LINES 174-177 Instead of just showing the velocities along the drifter positions, we can plot the mean model velocities over the same period as a background. So we interpolate them at 15 m with :func:`xoa.regrid.regrid1d` and compute a time average. .. GENERATED FROM PYTHON SOURCE LINES 177-182 .. code-block:: Python d15 = xr.DataArray([15.0], dims="depth", name="depth") uh15 = regrid1d(hycom["u"], d15).squeeze(drop=True).mean(dim="time") vh15 = regrid1d(hycom["v"], d15).squeeze(drop=True).mean(dim="time") .. GENERATED FROM PYTHON SOURCE LINES 183-188 Plot ---- The geographic extent is easily computed with the :func:`xoa.geo.get_extent` function. The background currents are plotted with the :func:`xoa.plot.plot_flow` function. .. GENERATED FROM PYTHON SOURCE LINES 188-222 .. code-block:: Python pmerc = ccrs.Mercator() pcarr = ccrs.PlateCarree() fig, ax = plt.subplots(figsize=(8, 7), subplot_kw={"facecolor": "teal", "projection": pmerc}) plt.subplots_adjust(right=0.85) ax.gridlines(draw_labels=True, dms=True) ax.set_extent(xgeo.get_extent(uh15)) kwqv = dict(scale_units="dots", scale=0.1 / 20, units="dots", transform=pcarr) plot_flow(uh15, vh15, axes=ax, transform=pcarr, color="w", alpha=(0.3, 1), linewidth=0.6) qv = ax.quiver( uloc.lon.values, uloc.lat.values, uloc.values, vloc.values, color="w", width=2, label="Model", **kwqv, ) ax.plot(drifter.lon.values, drifter.lat.values, '-', color="C1", transform=pcarr, lw=0.5) ax.quiver( drifter.lon.values, drifter.lat.values, drifter.u.values, drifter.v.values, color="C1", label="Drifter", width=2, **kwqv, ) plt.quiverkey(qv, 0.1, 1.06, 0.1, r"0.1 $m\,s^{-1}$", color="k", alpha=1, labelpos="E") plt.legend() plot_double_minimap(drifter) .. image-sg:: /examples/images/sphx_glr_plot_hycom_gdp_001.png :alt: plot hycom gdp :srcset: /examples/images/sphx_glr_plot_hycom_gdp_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (, ) .. GENERATED FROM PYTHON SOURCE LINES 223-227 The discrepancies between the lagrangian and mean eulerian currents highlight the variability due to the mesocale activity. The differences between the observed and modeled currents are partly due to the lack of data assimilation in the model. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 27.017 seconds) .. _sphx_glr_download_examples_plot_hycom_gdp.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/shom-fr/xoa/master?urlpath=lab/tree/notebooks/examples/plot_hycom_gdp.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_hycom_gdp.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_hycom_gdp.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_hycom_gdp.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_