.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/advection_correction.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_auto_examples_advection_correction.py: Advection correction ==================== This tutorial shows how to use the optical flow routines of pysteps to implement the advection correction procedure described in Anagnostou and Krajewski (1999). Advection correction is a temporal interpolation procedure that is often used when estimating rainfall accumulations to correct for the shift of rainfall patterns between consecutive radar rainfall maps. This shift becomes particularly significant for long radar scanning cycles and in presence of fast moving precipitation features. .. note:: The code for the advection correction using pysteps was originally written by `Daniel Wolfensberger `_. .. GENERATED FROM PYTHON SOURCE LINES 18-28 .. code-block:: Python from datetime import datetime import matplotlib.pyplot as plt import numpy as np from pysteps import io, motion, rcparams from pysteps.utils import conversion, dimension from pysteps.visualization import plot_precip_field from scipy.ndimage import map_coordinates .. GENERATED FROM PYTHON SOURCE LINES 29-39 Read the radar input images --------------------------- First, we import a sequence of 36 images of 5-minute radar composites that we will use to produce a 3-hour rainfall accumulation map. We will keep only one frame every 10 minutes, to simulate a longer scanning cycle and thus better highlight the need for advection correction. You need the pysteps-data archive downloaded and the pystepsrc file configured with the data_source paths pointing to data folders. .. GENERATED FROM PYTHON SOURCE LINES 39-44 .. code-block:: Python # Selected case date = datetime.strptime("201607112100", "%Y%m%d%H%M") data_source = rcparams.data_sources["mch"] .. GENERATED FROM PYTHON SOURCE LINES 45-47 Load the data from the archive ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 47-75 .. code-block:: Python root_path = data_source["root_path"] path_fmt = data_source["path_fmt"] fn_pattern = data_source["fn_pattern"] fn_ext = data_source["fn_ext"] importer_name = data_source["importer"] importer_kwargs = data_source["importer_kwargs"] timestep = data_source["timestep"] # Find the input files from the archive fns = io.archive.find_by_date( date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=35 ) # Read the radar composites importer = io.get_method(importer_name, "importer") R, __, metadata = io.read_timeseries(fns, importer, **importer_kwargs) # Convert to mm/h R, metadata = conversion.to_rainrate(R, metadata) # Upscale to 2 km (simply to reduce the memory demand) R, metadata = dimension.aggregate_fields_space(R, metadata, 2000) # Keep only one frame every 10 minutes (i.e., every 2 timesteps) # (to highlight the need for advection correction) R = R[::2] .. GENERATED FROM PYTHON SOURCE LINES 76-85 Advection correction -------------------- Now we need to implement the advection correction for a pair of successive radar images. The procedure is based on the algorithm described in Anagnostou and Krajewski (Appendix A, 1999). To evaluate the advection occurred between two successive radar images, we are going to use the Lucas-Kanade optical flow routine available in pysteps. .. GENERATED FROM PYTHON SOURCE LINES 85-116 .. code-block:: Python def advection_correction(R, T=5, t=1): """ R = np.array([qpe_previous, qpe_current]) T = time between two observations (5 min) t = interpolation timestep (1 min) """ # Evaluate advection oflow_method = motion.get_method("LK") fd_kwargs = {"buffer_mask": 10} # avoid edge effects V = oflow_method(np.log(R), fd_kwargs=fd_kwargs) # Perform temporal interpolation Rd = np.zeros((R[0].shape)) x, y = np.meshgrid( np.arange(R[0].shape[1], dtype=float), np.arange(R[0].shape[0], dtype=float) ) for i in range(t, T + t, t): pos1 = (y - i / T * V[1], x - i / T * V[0]) R1 = map_coordinates(R[0], pos1, order=1) pos2 = (y + (T - i) / T * V[1], x + (T - i) / T * V[0]) R2 = map_coordinates(R[1], pos2, order=1) Rd += (T - i) * R1 + i * R2 return t / T**2 * Rd .. GENERATED FROM PYTHON SOURCE LINES 117-119 Finally, we apply the advection correction to the whole sequence of radar images and produce the rainfall accumulation map. .. GENERATED FROM PYTHON SOURCE LINES 119-125 .. code-block:: Python R_ac = R[0].copy() for i in range(R.shape[0] - 1): R_ac += advection_correction(R[i : (i + 2)], T=10, t=1) R_ac /= R.shape[0] .. GENERATED FROM PYTHON SOURCE LINES 126-135 Results ------- We compare the two accumulation maps. The first map on the left is computed without advection correction and we can therefore see that the shift between successive images 10 minutes apart produces irregular accumulations. Conversely, the rainfall accumulation of the right is produced using advection correction to account for this spatial shift. The final result is a smoother rainfall accumulation map. .. GENERATED FROM PYTHON SOURCE LINES 135-144 .. code-block:: Python plt.figure(figsize=(9, 4)) plt.subplot(121) plot_precip_field(R.mean(axis=0), title="3-h rainfall accumulation") plt.subplot(122) plot_precip_field(R_ac, title="Same with advection correction") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_advection_correction_001.png :alt: 3-h rainfall accumulation, Same with advection correction :srcset: /auto_examples/images/sphx_glr_advection_correction_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 145-152 Reference ~~~~~~~~~ Anagnostou, E. N., and W. F. Krajewski. 1999. "Real-Time Radar Rainfall Estimation. Part I: Algorithm Formulation." Journal of Atmospheric and Oceanic Technology 16: 189–97. https://doi.org/10.1175/1520-0426(1999)016<0189:RTRREP>2.0.CO;2 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.253 seconds) .. _sphx_glr_download_auto_examples_advection_correction.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: advection_correction.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: advection_correction.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: advection_correction.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_