.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/example_wind_drift_coefficient_from_trajectory.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_example_wind_drift_coefficient_from_trajectory.py: Retieving wind drift factor from trajectory =========================================== .. GENERATED FROM PYTHON SOURCE LINES 6-14 .. code-block:: Python from datetime import datetime, timedelta import numpy as np import matplotlib.pyplot as plt import cmocean from opendrift.models.oceandrift import OceanDrift from opendrift.models.physics_methods import wind_drift_factor_from_trajectory, distance_between_trajectories, skillscore_liu_weissberg .. GENERATED FROM PYTHON SOURCE LINES 15-20 A very simple drift model is: current + wind_drift_factor*wind where wind_drift_factor for surface drift is typically 0.033 if Stokes drift included, and 0.02 in addition to Stokes drift. This example illustrates how a best-fit wind_drift_factor can be calculated from an observed trajectory, using two different methods. .. GENERATED FROM PYTHON SOURCE LINES 22-23 First we simulate a synthetic drifter trajectory .. GENERATED FROM PYTHON SOURCE LINES 23-28 .. code-block:: Python ot = OceanDrift(loglevel=50) ot.add_readers_from_list([ot.test_data_folder() + '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc', ot.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) .. GENERATED FROM PYTHON SOURCE LINES 29-30 Adding some horizontal diffusivity as "noise" .. GENERATED FROM PYTHON SOURCE LINES 30-32 .. code-block:: Python ot.set_config('drift:horizontal_diffusivity', 10) .. GENERATED FROM PYTHON SOURCE LINES 33-34 Using a wind_drift_factor of 0.33 i.e. drift is current + 3.3% of wind speed .. GENERATED FROM PYTHON SOURCE LINES 34-37 .. code-block:: Python ot.seed_elements(lon=4, lat=60, number=1, time=ot.env.readers[list(ot.env.readers)[0]].start_time, wind_drift_factor=0.033) .. GENERATED FROM PYTHON SOURCE LINES 38-40 .. code-block:: Python ot.run(duration=timedelta(hours=12), time_step=600) .. GENERATED FROM PYTHON SOURCE LINES 41-42 Secondly, calculating the wind_drift_factor which reproduces the "observed" trajectory with minimal difference .. GENERATED FROM PYTHON SOURCE LINES 42-60 .. code-block:: Python drifter_lons = ot.history['lon'][0] drifter_lats = ot.history['lat'][0] drifter_times = ot.get_time_array()[0] drifter={'lon': drifter_lons, 'lat': drifter_lats, 'time': drifter_times, 'linewidth': 2, 'color': 'b', 'label': 'Synthetic drifter'} o = OceanDrift(loglevel=50) o.add_readers_from_list([o.test_data_folder() + '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc', o.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False) t = o.env.get_variables_along_trajectory(variables=['x_sea_water_velocity', 'y_sea_water_velocity', 'x_wind', 'y_wind'], lons=drifter_lons, lats=drifter_lats, times=drifter_times) wind_drift_factor, azimuth = wind_drift_factor_from_trajectory(t) o.seed_elements(lon=4, lat=60, number=1, time=ot.env.readers[list(ot.env.readers)[0]].start_time, wind_drift_factor=0.033) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) .. GENERATED FROM PYTHON SOURCE LINES 61-62 New simulation, this time without diffusivity/noise .. GENERATED FROM PYTHON SOURCE LINES 62-66 .. code-block:: Python o.run(duration=timedelta(hours=12), time_step=600) o.plot(fast=True, legend=True, drifter=drifter) .. image-sg:: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_001.png :alt: OpenDrift - OceanDrift 2015-11-16 00:00 to 2015-11-16 12:00 UTC (73 steps) :srcset: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (,
) .. GENERATED FROM PYTHON SOURCE LINES 67-71 The mean retrieved wind_drift_factor is 0.036, slichtly higher than the value 0.033 used for the simulation. The difference is due to the "noise" added by horizontal diffusion. Note that the retieved wind_drift_factor is linked to the wind and current used for the retrieval, other forcing datasets will yeld different value of the wind_drift_factor. .. GENERATED FROM PYTHON SOURCE LINES 71-82 .. code-block:: Python print(wind_drift_factor.mean()) plt.hist(wind_drift_factor, label='Retrieved wind_drift_factor') plt.axvline(x=0.033, label='Actual wind_drift_factor of 0.033', color='k') plt.axvline(x=wind_drift_factor.mean(), label='Mean retieved wind_drift_factor of %.3f' % wind_drift_factor.mean(), color='r') plt.ylabel('Number') plt.xlabel('Wind drift factor [fraction]') plt.legend(loc='lower left') plt.show() .. image-sg:: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_002.png :alt: example wind drift coefficient from trajectory :srcset: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.03639823255601244 .. GENERATED FROM PYTHON SOURCE LINES 83-85 A polar 2D histogram showing also the azimuth offset direction of the retrieved wind drift factor. See Sutherland et al. (2020), https://doi.org/10.1175/JTECH-D-20-0013.1 .. GENERATED FROM PYTHON SOURCE LINES 85-100 .. code-block:: Python wmax = wind_drift_factor.max() wbins = np.arange(0, wmax+.005, .005) abins = np.linspace(-180, 180, 30) hist, _, _ = np.histogram2d(azimuth, wind_drift_factor, bins=(abins, wbins)) A, W = np.meshgrid(abins, wbins) fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) ax.set_theta_zero_location('N', offset=0) ax.set_theta_direction(-1) pc = ax.pcolormesh(np.radians(A), W, hist.T, cmap=cmocean.cm.dense) plt.arrow(np.pi, wmax, 0, -wmax, width=.015, facecolor='k', zorder=100, head_width=.8, lw=2, head_length=.005, length_includes_head=True) plt.text(np.pi, wmax*.5, ' Wind direction', fontsize=18) plt.grid() plt.show() .. image-sg:: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_003.png :alt: example wind drift coefficient from trajectory :srcset: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 101-107 Alternative method, using skillscore ==================================== Here we calculate trajectories for a range of wind_drift_factors, and calculate the Liu-Weissberg skillscore for each. The optimized wind_drift_factor corresponds to the maximum skillscore. .. GENERATED FROM PYTHON SOURCE LINES 107-113 .. code-block:: Python o = OceanDrift(loglevel=50) o.add_readers_from_list([o.test_data_folder() + '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc', o.test_data_folder() + '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc'], lazy=False) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) /opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems proj = self._crs.to_proj4(version=version) .. GENERATED FROM PYTHON SOURCE LINES 114-115 Calulating trajectories for 100 different wind_drift_factors between 0 and 0.05 .. GENERATED FROM PYTHON SOURCE LINES 115-121 .. code-block:: Python wdf = np.linspace(0.0, 0.05, 100) o.seed_elements(lon=4, lat=60, time=ot.env.readers[list(ot.env.readers)[0]].start_time, wind_drift_factor=wdf, number=len(wdf)) o.run(duration=timedelta(hours=12), time_step=600) o.plot(linecolor='wind_drift_factor', drifter=drifter) .. image-sg:: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_004.png :alt: OpenDrift - OceanDrift 2015-11-16 00:00 to 2015-11-16 12:00 UTC (73 steps) :srcset: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (,
) .. GENERATED FROM PYTHON SOURCE LINES 122-123 Plotting and finding the wind_drift_factor which maximises the skillscore .. GENERATED FROM PYTHON SOURCE LINES 123-131 .. code-block:: Python skillscore = o.skillscore_trajectory(drifter_lons, drifter_lats, drifter_times, tolerance_threshold=1) ind = np.argmax(skillscore) plt.plot(wdf, skillscore) plt.xlabel('Wind drift factor [fraction]') plt.ylabel('Liu-Weissberg skillscore') plt.title('Maximum skillscore %.3f for wdf=%.3f' % (skillscore[ind], wdf[ind])) plt.show() .. image-sg:: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_005.png :alt: Maximum skillscore 0.962 for wdf=0.034 :srcset: /gallery/images/sphx_glr_example_wind_drift_coefficient_from_trajectory_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-138 We see that using skillscore from the full trajectories gives a better estimate than calculating the average of the wind_drift_factor from one position to the next (i.e. polar histogram above). This is even more clear if increasing the diffusivity (i.e. noise) above from 10 m2/s to 200 m2/s: The histogram method then gives 0.071, which is much to high (true is 0.033), and the histogram is noisy. The skillscore method is still robust, and gives a `wind_drift_factor` of 0.036, only slightly too high. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 34.735 seconds) .. _sphx_glr_download_gallery_example_wind_drift_coefficient_from_trajectory.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_wind_drift_coefficient_from_trajectory.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_wind_drift_coefficient_from_trajectory.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_