.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/nwp-radar_qpf_ms.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_nwp-radar_qpf_ms.py: Radar Advection with NWP Blend (Malaysia) =============================================================== This example shows the technique of blending Numerical Weather Forecast (NWP) with COM-SWIRLS output to generate operational deterministic QPF up to 3 and 6 hours ahead. NWP data was generated using the WRF Model while radar data was taken from the Malaysian radar observation network .. GENERATED FROM PYTHON SOURCE LINES 13-16 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 18-19 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 19-57 .. code-block:: default # Python package to allow system command line functions import os # Python package to manage warning message import warnings # Python package for timestamp import pandas as pd # Python package for xarrays to read and handle netcdf data import xarray as xr # Python package for numerical calculations import numpy as np # Python package for reading map shape file import cartopy.io.shapereader as shpreader # Python package for land/sea features import cartopy.feature as cfeature # Python package for projection import cartopy.crs as ccrs # Python package for creating plots from matplotlib import pyplot as plt # Python package for output import grid from matplotlib.gridspec import GridSpec # Python package for colorbars from matplotlib.colors import BoundaryNorm, ListedColormap # Python package for scalar data to RGBA mapping from matplotlib.cm import ScalarMappable # Python com-swirls package to standardize attributes from swirlspy.utils import standardize_attr, FrameType # Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection from swirlspy.qpf import rover, sla # Python com-swirls package to blend nwp and nowcast (RaINS) from swirlspy.blending import rains, nwp_bias_correction # directory constants from swirlspy.tests.samples import DATA_DIR from swirlspy.tests.outputs import OUTPUT_DIR warnings.filterwarnings("ignore") start_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 58-66 Initializing --------------------------------------------------- This section demonstrates the extraction of radar & nwp data from netcdf into python Step 1: Define your input data directory and output directory .. GENERATED FROM PYTHON SOURCE LINES 66-72 .. code-block:: default # Supply the directory of radar and nwp data data_dir = os.path.abspath( os.path.join(DATA_DIR, 'netcdf_ms') ) .. GENERATED FROM PYTHON SOURCE LINES 73-75 Step 2: Define a basetime .. GENERATED FROM PYTHON SOURCE LINES 75-79 .. code-block:: default # Supply basetime basetime = pd.Timestamp('201908090900') .. GENERATED FROM PYTHON SOURCE LINES 80-82 Step 3: Read data files from the radar data using xarray() .. GENERATED FROM PYTHON SOURCE LINES 82-105 .. code-block:: default # Radar data listed from the basetime[0] --> 3 hours before the basetime[17] (descending time) interval = 10 # Interval of radar data radar_datas = [] for i in range(0, 18): t = basetime - pd.Timedelta(minutes=i * interval) # Radar data nomenclature filename = os.path.join( data_dir, t.strftime("radar_d03_%Y-%m-%d_%H_%M_00.rapids.nc") ) reflec = xr.open_dataset(filename) radar_datas.append(reflec) # Concatenate list by time reflec_concat = xr.concat(radar_datas, dim='time') # Extracting the radar data: The radar dBZ variable is named 'Zradar', therefore, we extract 'Zradar' radar = reflec_concat['Zradar'] # Reversing such that time goes from earliest to latest; 3 hours before basetime[0] --> basetime[17] radar = radar.sortby('time', ascending=True) .. GENERATED FROM PYTHON SOURCE LINES 106-108 Step 4: Reading nwp netcdf data into xarray .. GENERATED FROM PYTHON SOURCE LINES 108-134 .. code-block:: default # NWP data listed from the basetime[0] --> 3 and 6 hours AFTER basetime[17] (nowcast) nwp_3hr_datas = [] nwp_6hr_datas = [] for i in range(0, 36): t = basetime + pd.Timedelta(minutes=i * interval) # Radar nwp nomenclature filename = os.path.join( data_dir, t.strftime("wrfout_d03_%Y-%m-%d_%H_%M_00.rapids.nc") ) reflec = xr.open_dataset(filename) nwp_6hr_datas.append(reflec) if i < 18: nwp_3hr_datas.append(reflec) # Concatenating the nwp reflectivity list of data by time reflec_3hr_concat = xr.concat(nwp_3hr_datas, dim='time') reflec_6hr_concat = xr.concat(nwp_6hr_datas, dim='time') # Extracting the nwp data: The nwp dBZ variable is called 'zwrf', extracting 'zwrf' nwp_3hr = reflec_3hr_concat['zwrf'] nwp_6hr = reflec_6hr_concat['zwrf'] initialising_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 135-139 Step 5: Bias correction of nwp data The objective of bias correction is to match the nwp percentile to the radar percentile This is also known as frequency matching. .. GENERATED FROM PYTHON SOURCE LINES 139-145 .. code-block:: default nwp_3hr_corrected = nwp_bias_correction(radar, nwp_3hr) nwp_6hr_corrected = nwp_bias_correction(radar, nwp_6hr) bias_correction_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 146-152 Nowcast (SWIRLS-Radar-Advection) --------------------------------------------------- The swirls radar advection was performed using the observed radar data Firstly, some attributes necessary for com-swirls input variable is added Secondly, rover function is invoked to compute the motion field Thirdly, semi-lagrangian advection is performed to advect the radar data using the rover motion field .. GENERATED FROM PYTHON SOURCE LINES 152-167 .. code-block:: default # Adding in some attributes that is step_size <10 mins in pandas.Timedelta>, zero_value <9999.> frame_type radar.attrs['step_size'] = pd.Timedelta(minutes=10) standardize_attr(radar, frame_type=FrameType.dBZ, zero_value=np.nan) # Rover motion field computation motion = rover(radar) rover_time = pd.Timestamp.now() # Semi-Lagrangian Advection swirls = sla(radar, motion, 35) # Radar time goes from earliest to latest sla_time = pd.Timestamp.now() .. rst-class:: sphx-glr-script-out .. code-block:: none RUNNING 'rover' FOR EXTRAPOLATION..... .. GENERATED FROM PYTHON SOURCE LINES 168-172 Blending (RaINS) --------------------------------------------------- Blending Numerical Weather Forecast (NWP) with COM-SWIRLS output to generate operational deterministic QPF up to 3 hours ahead. .. GENERATED FROM PYTHON SOURCE LINES 172-178 .. code-block:: default blended_3hr = rains(nwp_3hr_corrected, swirls[0:18]) blended_6hr = rains(nwp_6hr_corrected, swirls) rains_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 179-182 Plotting result --------------------------------------------------- Step 1: Defining the dBZ levels, colorbar parameters and projection .. GENERATED FROM PYTHON SOURCE LINES 182-233 .. code-block:: default # levels of colorbar (dBZ) levels = [-32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 53, 56, 58, 60, 62] # hko colormap for dBZ at each levels cmap = ListedColormap([ '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433', '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200', '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0' ]) # boundary norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) # scalar data to RGBA mapping scalar_map = ScalarMappable(cmap=cmap, norm=norm) scalar_map.set_array([]) # Defining plot parameters map_shape_file = os.path.abspath(os.path.join( DATA_DIR, 'shape/se_asia_s' )) # coastline and province se_asia = cfeature.ShapelyFeature( list(shpreader.Reader(map_shape_file).geometries()), ccrs.PlateCarree() ) # output area extents = [99, 120, 0.5, 7.25] # base_map plotting function def plot_base(ax: plt.Axes): ax.set_extent(extents, crs=ccrs.PlateCarree()) # fake the ocean color ax.imshow(np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]), origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -180, 180], zorder=-1) # coastline, state, color ax.add_feature(se_asia, facecolor=cfeature.COLORS['land'], edgecolor='none', zorder=0) # overlay coastline, state without color ax.add_feature(se_asia, facecolor='none', edgecolor='gray', linewidth=0.5) ax.set_title('') .. GENERATED FROM PYTHON SOURCE LINES 234-236 Step 2: Filtering values <= 5dbZ are not plotted .. GENERATED FROM PYTHON SOURCE LINES 236-241 .. code-block:: default nwp_3hr_corrected = nwp_3hr_corrected.where(nwp_3hr_corrected > 5, np.nan) blended_3hr = blended_3hr.where(blended_3hr > 5, np.nan) nwp_6hr_corrected = nwp_6hr_corrected.where(nwp_6hr_corrected > 5, np.nan) blended_6hr = blended_6hr.where(blended_6hr > 5, np.nan) .. GENERATED FROM PYTHON SOURCE LINES 242-244 Step 3: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 3 hours ahead .. GENERATED FROM PYTHON SOURCE LINES 244-316 .. code-block:: default fig: plt.Figure = plt.figure( figsize=(3 * 5 + 1, 3 * 2), frameon=False ) gs = GridSpec( 3, 3, figure=fig, wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845 ) for row in range(3): time_index = (row + 1) * 6 - 1 timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm') for col in range(3): ax: plt.Axes = fig.add_subplot( gs[row, col], projection=ccrs.PlateCarree() ) if col % 3 == 0: z = swirls[time_index].values lats = swirls[time_index].latitude lons = swirls[time_index].longitude title = 'SWIRLS Reflectivity' elif col % 3 == 1: z = nwp_3hr_corrected[time_index].values lats = nwp_3hr_corrected[time_index].lat lons = nwp_3hr_corrected[time_index].lon title = 'NWP Reflectivity' elif col % 3 == 2: z = blended_3hr[time_index].values lats = blended_3hr[time_index].latitude lons = blended_3hr[time_index].longitude title = 'Blended Reflectivity' # plot base map plot_base(ax) # plot reflectivity ax.contourf( lons, lats, z, 60, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm, levels=levels ) ax.set_title( f"{title}\n" + f"Initial @ {basetime.strftime('%H:%MZ')}", loc='left', fontsize=9 ) ax.set_title('') ax.set_title( f"Initial {basetime.strftime('%Y-%m-%d')} \n" + f"Forecast Valid @ {timelabel.strftime('%H:%MZ')} ", loc='right', fontsize=9 ) cbar_ax = fig.add_axes([0.85, 0.105, 0.01, 0.845]) cbar = fig.colorbar( scalar_map, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g' ) cbar.ax.set_ylabel('Reflectivity (dBZ)', rotation=90) fig.savefig( os.path.join( OUTPUT_DIR, "swirls_nwp_blend_ms_3hr.png" ), dpi=450, bbox_inches="tight", pad_inches=0.1 ) .. image-sg:: /auto_examples/images/sphx_glr_nwp-radar_qpf_ms_001.png :alt: SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z :srcset: /auto_examples/images/sphx_glr_nwp-radar_qpf_ms_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 317-319 Step 4: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 6 hours ahead .. GENERATED FROM PYTHON SOURCE LINES 319-393 .. code-block:: default fig: plt.Figure = plt.figure( figsize=(3 * 5 + 1, 6 * 2), frameon=False ) gs = GridSpec( 6, 3, figure=fig, wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845 ) for row in range(6): time_index = (row + 1) * 6 - 1 timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm') for col in range(3): ax: plt.Axes = fig.add_subplot( gs[row, col], projection=ccrs.PlateCarree() ) if col % 3 == 0: z = swirls[time_index].values lats = swirls[time_index].latitude lons = swirls[time_index].longitude title = 'SWIRLS Reflectivity' elif col % 3 == 1: z = nwp_6hr_corrected[time_index].values lats = nwp_6hr_corrected[time_index].lat lons = nwp_6hr_corrected[time_index].lon title = 'NWP Reflectivity' elif col % 3 == 2: z = blended_6hr[time_index].values lats = blended_6hr[time_index].latitude lons = blended_6hr[time_index].longitude title = 'Blended Reflectivity' # plot base map plot_base(ax) # plot reflectivity ax.contourf( lons, lats, z, 60, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm, levels=levels ) ax.set_title( f"{title}\n" + f"Initial @ {basetime.strftime('%H:%MZ')}", loc='left', fontsize=9 ) ax.set_title('') ax.set_title( f"Initial {basetime.strftime('%Y-%m-%d')} \n" + f"Forecast Valid @ {timelabel.strftime('%H:%MZ')} ", loc='right', fontsize=9 ) cbar_ax = fig.add_axes([0.85, 0.105, 0.01, 0.845]) cbar = fig.colorbar( scalar_map, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g' ) cbar.ax.set_ylabel('Reflectivity (dBZ)', rotation=90) fig.savefig( os.path.join( OUTPUT_DIR, "swirls_nwp_blend_ms_6hr.png" ), dpi=450, bbox_inches="tight", pad_inches=0.1 ) radar_image_time = pd.Timestamp.now() .. image-sg:: /auto_examples/images/sphx_glr_nwp-radar_qpf_ms_002.png :alt: SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 10:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 11:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 12:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 13:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 13:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 13:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 14:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 14:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 14:00Z , SWIRLS Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 15:00Z , NWP Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 15:00Z , Blended Reflectivity Initial @ 09:00Z, Initial 2019-08-09 Forecast Valid @ 15:00Z :srcset: /auto_examples/images/sphx_glr_nwp-radar_qpf_ms_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 394-397 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 397-412 .. code-block:: default print(f"Start time: {start_time}") print(f"Initialising time: {initialising_time}") print(f"NWP bias correction time: {bias_correction_time}") print(f"SLA time: {sla_time}") print(f"RaINS time: {rains_time}") print(f"Plotting radar image time: {radar_image_time}") print(f"Time to initialise: {initialising_time - start_time}") print( f"Time to run NWP bias correction: {bias_correction_time - initialising_time}") print(f"Time to run rover: {rover_time - bias_correction_time}") print(f"Time to perform SLA: {sla_time - rover_time}") print(f"Time to perform RaINS: {rains_time - sla_time}") print(f"Time to plot radar image: {radar_image_time - rains_time}") .. rst-class:: sphx-glr-script-out .. code-block:: none Start time: 2024-04-18 02:53:20.666814 Initialising time: 2024-04-18 02:53:24.618740 NWP bias correction time: 2024-04-18 02:54:08.382142 SLA time: 2024-04-18 02:59:55.747661 RaINS time: 2024-04-18 03:00:06.088907 Plotting radar image time: 2024-04-18 03:01:41.300048 Time to initialise: 0 days 00:00:03.951926 Time to run NWP bias correction: 0 days 00:00:43.763402 Time to run rover: 0 days 00:00:01.101595 Time to perform SLA: 0 days 00:05:46.263924 Time to perform RaINS: 0 days 00:00:10.341246 Time to plot radar image: 0 days 00:01:35.211141 .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 8 minutes 23.145 seconds) .. _sphx_glr_download_auto_examples_nwp-radar_qpf_ms.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: nwp-radar_qpf_ms.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: nwp-radar_qpf_ms.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_