.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_linear_blending.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_plot_linear_blending.py: Linear blending =============== This tutorial shows how to construct a simple linear blending between a STEPS ensemble nowcast and a Numerical Weather Prediction (NWP) rainfall forecast. The used datasets are from the Bureau of Meteorology, Australia. .. GENERATED FROM PYTHON SOURCE LINES 11-23 .. code-block:: Python import os from datetime import datetime from matplotlib import pyplot as plt import pysteps from pysteps import io, rcparams, nowcasts, blending from pysteps.utils import conversion from pysteps.visualization import plot_precip_field .. GENERATED FROM PYTHON SOURCE LINES 24-34 Read the radar images and the NWP forecast ------------------------------------------ First, we import a sequence of 3 images of 10-minute radar composites and the corresponding NWP rainfall forecast that was available at that time. You need the pysteps-data archive downloaded and the pystepsrc file configured with the data_source paths pointing to data folders. Additionally, the pysteps-nwp-importers plugin needs to be installed, see https://github.com/pySTEPS/pysteps-nwp-importers. .. GENERATED FROM PYTHON SOURCE LINES 34-42 .. code-block:: Python # Selected case date_radar = datetime.strptime("202010310400", "%Y%m%d%H%M") # The last NWP forecast was issued at 00:00 date_nwp = datetime.strptime("202010310000", "%Y%m%d%H%M") radar_data_source = rcparams.data_sources["bom"] nwp_data_source = rcparams.data_sources["bom_nwp"] .. GENERATED FROM PYTHON SOURCE LINES 43-45 Load the data from the archive ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 45-80 .. code-block:: Python root_path = radar_data_source["root_path"] path_fmt = "prcp-c10/66/%Y/%m/%d" fn_pattern = "66_%Y%m%d_%H%M00.prcp-c10" fn_ext = radar_data_source["fn_ext"] importer_name = radar_data_source["importer"] importer_kwargs = radar_data_source["importer_kwargs"] timestep = 10.0 # Find the radar files in the archive fns = io.find_by_date( date_radar, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2 ) # Read the radar composites importer = io.get_method(importer_name, "importer") radar_precip, _, radar_metadata = io.read_timeseries(fns, importer, **importer_kwargs) # Import the NWP data filename = os.path.join( nwp_data_source["root_path"], datetime.strftime(date_nwp, nwp_data_source["path_fmt"]), datetime.strftime(date_nwp, nwp_data_source["fn_pattern"]) + "." + nwp_data_source["fn_ext"], ) nwp_importer = io.get_method("bom_nwp", "importer") nwp_precip, _, nwp_metadata = nwp_importer(filename) # Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00) # End of the forecast is 18 time steps (+3 hours) in advance. precip_nwp = nwp_precip[24:43, :, :] .. rst-class:: sphx-glr-script-out .. code-block:: none Rainfall values are accumulated. Disaggregating by time step .. GENERATED FROM PYTHON SOURCE LINES 81-83 Pre-processing steps -------------------- .. GENERATED FROM PYTHON SOURCE LINES 83-125 .. code-block:: Python # Make sure the units are in mm/h converter = pysteps.utils.get_method("mm/h") radar_precip, radar_metadata = converter(radar_precip, radar_metadata) precip_nwp, nwp_metadata = converter(precip_nwp, nwp_metadata) # Threshold the data radar_precip[radar_precip < 0.1] = 0.0 precip_nwp[precip_nwp < 0.1] = 0.0 # Plot the radar rainfall field and the first time step of the NWP forecast. # For the initial time step (t=0), the NWP rainfall forecast is not that different # from the observed radar rainfall, but it misses some of the locations and # shapes of the observed rainfall fields. Therefore, the NWP rainfall forecast will # initially get a low weight in the blending process. date_str = datetime.strftime(date_radar, "%Y-%m-%d %H:%M") plt.figure(figsize=(10, 5)) plt.subplot(121) plot_precip_field( radar_precip[-1, :, :], geodata=radar_metadata, title=f"Radar observation at {date_str}", ) plt.subplot(122) plot_precip_field( precip_nwp[0, :, :], geodata=nwp_metadata, title=f"NWP forecast at {date_str}" ) plt.tight_layout() plt.show() # Only keep the NWP forecasts from 2020-10-31 04:05 onwards, because the first # forecast lead time starts at 04:05. precip_nwp = precip_nwp[1:] # Transform the radar data to dB - this transformation is useful for the motion # field estimation and the subsequent nowcasts. The NWP forecast is not # transformed, because the linear blending code sets everything back in mm/h # after the nowcast. transformer = pysteps.utils.get_method("dB") radar_precip, radar_metadata = transformer(radar_precip, radar_metadata, threshold=0.1) .. image-sg:: /auto_examples/images/sphx_glr_plot_linear_blending_001.png :alt: Radar observation at 2020-10-31 04:00, NWP forecast at 2020-10-31 04:00 :srcset: /auto_examples/images/sphx_glr_plot_linear_blending_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 126-128 Determine the velocity field for the radar rainfall nowcast ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 128-133 .. code-block:: Python oflow_method = pysteps.motion.get_method("lucaskanade") velocity_radar = oflow_method(radar_precip) .. GENERATED FROM PYTHON SOURCE LINES 134-136 The linear blending of nowcast and NWP rainfall forecast -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 136-152 .. code-block:: Python # Calculate the blended precipitation field precip_blended = blending.linear_blending.forecast( precip=radar_precip[-1, :, :], precip_metadata=radar_metadata, velocity=velocity_radar, timesteps=18, timestep=10, nowcast_method="extrapolation", # simple advection nowcast precip_nwp=precip_nwp, precip_nwp_metadata=nwp_metadata, start_blending=60, # in minutes (this is an arbritrary choice) end_blending=120, # in minutes (this is an arbritrary choice) ) .. GENERATED FROM PYTHON SOURCE LINES 153-163 The salient blending of nowcast and NWP rainfall forecast --------------------------------------------------------- This method follows the saliency-based blending procedure described in :cite:`Hwang2015`. The blending is based on intensities and forecast times. The blended product preserves pixel intensities with time if they are strong enough based on their ranked salience. Saliency is the property of an object to be outstanding with respect to its surroundings. The ranked salience is calculated by first determining the difference in the normalized intensity of the nowcasts and NWP. Next, the pixel intensities are ranked, in which equally comparable values receive the same ranking number. .. GENERATED FROM PYTHON SOURCE LINES 163-180 .. code-block:: Python # Calculate the salient blended precipitation field precip_salient_blended = blending.linear_blending.forecast( precip=radar_precip[-1, :, :], precip_metadata=radar_metadata, velocity=velocity_radar, timesteps=18, timestep=10, nowcast_method="extrapolation", # simple advection nowcast precip_nwp=precip_nwp, precip_nwp_metadata=nwp_metadata, start_blending=60, # in minutes (this is an arbritrary choice) end_blending=120, # in minutes (this is an arbritrary choice) saliency=True, ) .. GENERATED FROM PYTHON SOURCE LINES 181-183 Visualize the output -------------------- .. GENERATED FROM PYTHON SOURCE LINES 185-186 Calculate the radar rainfall nowcasts for visualization .. GENERATED FROM PYTHON SOURCE LINES 186-197 .. code-block:: Python nowcast_method_func = nowcasts.get_method("extrapolation") precip_nowcast = nowcast_method_func( precip=radar_precip[-1, :, :], velocity=velocity_radar, timesteps=18, ) # Make sure that precip_nowcast are in mm/h precip_nowcast, _ = conversion.to_rainrate(precip_nowcast, metadata=radar_metadata) .. GENERATED FROM PYTHON SOURCE LINES 198-207 The linear blending starts at 60 min, so during the first 60 minutes the blended forecast only consists of the extrapolation forecast (consisting of an extrapolation nowcast). Between 60 and 120 min, the NWP forecast gradually gets more weight, whereas the extrapolation forecasts gradually gets less weight. In addition, the saliency-based blending takes also the difference in pixel intensities into account, which are preserved over time if they are strong enough based on their ranked salience. Furthermore, pixels with relative low intensities get a lower weight and stay smaller in the saliency-based blending compared to linear blending. After 120 min, the blended forecast entirely consists of the NWP rainfall forecast. .. GENERATED FROM PYTHON SOURCE LINES 207-256 .. code-block:: Python fig = plt.figure(figsize=(8, 12)) leadtimes_min = [30, 60, 80, 100, 120] n_leadtimes = len(leadtimes_min) for n, leadtime in enumerate(leadtimes_min): # Extrapolation plt.subplot(n_leadtimes, 4, n * 4 + 1) plot_precip_field( precip_nowcast[int(leadtime / timestep) - 1, :, :], geodata=radar_metadata, title=f"Nowcast + {leadtime} min", axis="off", colorbar=False, ) # Nowcast with blending into NWP plt.subplot(n_leadtimes, 4, n * 4 + 2) plot_precip_field( precip_blended[int(leadtime / timestep) - 1, :, :], geodata=radar_metadata, title=f"Linear + {leadtime} min", axis="off", colorbar=False, ) # Nowcast with salient blending into NWP plt.subplot(n_leadtimes, 4, n * 4 + 3) plot_precip_field( precip_salient_blended[int(leadtime / timestep) - 1, :, :], geodata=radar_metadata, title=f"Salient + {leadtime} min", axis="off", colorbar=False, ) # Raw NWP forecast plt.subplot(n_leadtimes, 4, n * 4 + 4) plot_precip_field( precip_nwp[int(leadtime / timestep) - 1, :, :], geodata=nwp_metadata, title=f"NWP + {leadtime} min", axis="off", colorbar=False, ) plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_linear_blending_002.png :alt: Nowcast + 30 min, Linear + 30 min, Salient + 30 min, NWP + 30 min, Nowcast + 60 min, Linear + 60 min, Salient + 60 min, NWP + 60 min, Nowcast + 80 min, Linear + 80 min, Salient + 80 min, NWP + 80 min, Nowcast + 100 min, Linear + 100 min, Salient + 100 min, NWP + 100 min, Nowcast + 120 min, Linear + 120 min, Salient + 120 min, NWP + 120 min :srcset: /auto_examples/images/sphx_glr_plot_linear_blending_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 257-259 Note that the NaN values of the extrapolation forecast are replaced with NWP data in the blended forecast, even before the blending starts. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.258 seconds) .. _sphx_glr_download_auto_examples_plot_linear_blending.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_linear_blending.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_linear_blending.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_linear_blending.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_