.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/nwp-radar_qpf_vn.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_vn.py: Radar Advection with NWP Blend (Vietnam) =============================================================== 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. .. GENERATED FROM PYTHON SOURCE LINES 10-13 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-16 Import all required modules and methods: .. GENERATED FROM PYTHON SOURCE LINES 16-54 .. 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 55-63 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 63-69 .. code-block:: default # Supply the directory of radar and nwp data data_dir = os.path.abspath( os.path.join(DATA_DIR, 'vnmha') ) .. GENERATED FROM PYTHON SOURCE LINES 70-72 Step 2: Define a basetime .. GENERATED FROM PYTHON SOURCE LINES 72-76 .. code-block:: default # Supply basetime basetime = pd.Timestamp('202001240300') .. GENERATED FROM PYTHON SOURCE LINES 77-79 Step 3: Read data files from the radar data using xarray() .. GENERATED FROM PYTHON SOURCE LINES 79-103 .. 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(19): t = basetime - pd.Timedelta(minutes=i * interval) # Radar data nomenclature filename = os.path.join( data_dir, 'uf_vietnam', t.strftime("composite_uf_vn_%Y%m%d%H%M.nc") ) reflec = xr.open_dataset(filename).assign_coords(time=[t]) 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 '__xarray_dataarray_variable__', therefore, we extract '__xarray_dataarray_variable__' radar = reflec_concat['__xarray_dataarray_variable__'] # 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 104-106 Step 4: Reading nwp netcdf data into xarray .. GENERATED FROM PYTHON SOURCE LINES 106-131 .. code-block:: default # NWP data listed from the basetime[0] --> 3 hours AFTER basetime[18] (nowcast) nwp_datas = [] for i in range(1, 19): t = basetime + pd.Timedelta(minutes=i * interval) # Radar nwp nomenclature filename = os.path.join( data_dir, 'wrf_dbz_vietnam', t.strftime("wrfout_d01_%Y-%m-%d_%H_%M_00.nc") ) reflec = xr.open_dataset(filename).assign_coords(time=[t]) nwp_datas.append(reflec) # Concatenating the nwp reflectivity list of data by time reflec_concat = xr.concat(nwp_datas, dim='time') # Extracting the nwp data: The nwp dBZ variable is called 'zwrf', extracting 'zwrf' nwp = reflec_concat['zwrf'].rename({'lon': 'x', 'lat': 'y'}) nwp.sortby('y', ascending=False).sortby('time', ascending=True) # Remove "invalid" data nwp.values[nwp.values < 0] = np.nan initialising_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 132-136 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 136-143 .. code-block:: default nwp_corrected = nwp_bias_correction( radar, nwp, proj4_str='+proj=longlat +datum=WGS84 +no_defs' ) bias_correction_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 144-150 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 150-165 .. 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, 18) # 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 166-170 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 170-174 .. code-block:: default blended = rains(nwp_corrected, swirls[1:]) # skip swirls base frame rains_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 175-178 Plotting result --------------------------------------------------- Step 1: Defining the dBZ levels, colorbar parameters and projection .. GENERATED FROM PYTHON SOURCE LINES 178-232 .. 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 = [ blended.x.values[0], blended.x.values[-1], blended.y.values[0], blended.y.values[-1] ] # 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 233-235 Step 2: Filtering values <= 5dbZ are not plotted .. GENERATED FROM PYTHON SOURCE LINES 235-239 .. code-block:: default nwp_corrected = nwp_corrected.where(nwp_corrected > 5, np.nan) blended = blended.where(blended > 5, np.nan) swirls = swirls[1:] # skip base time frame .. GENERATED FROM PYTHON SOURCE LINES 240-242 Step 3: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 3 hours ahead .. GENERATED FROM PYTHON SOURCE LINES 242-344 .. code-block:: default plot_height = 3 plot_width = 2 n_rows = 2 n_column = 3 fig: plt.Figure = plt.figure( figsize=(plot_width * n_column + 1, plot_height * n_rows), frameon=False ) gs = GridSpec( n_rows, n_column, figure=fig, wspace=0.05, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845 ) for row in range(n_rows): time_index = (row + 1) * 6 - 1 timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm') for col in range(n_column): ax: plt.Axes = fig.add_subplot( gs[row, col], projection=ccrs.PlateCarree() ) if col % 3 == 0: z = swirls[time_index].values y = swirls[time_index].y x = swirls[time_index].x if row == 0: ax.text( 0, 1, 'SWIRLS Reflectivity', fontsize=8, ha='left', va='bottom', transform=ax.transAxes ) elif col % 3 == 1: z = nwp_corrected[time_index].values y = nwp_corrected[time_index].y x = nwp_corrected[time_index].x if row == 0: ax.text( 0, 1, 'NWP Reflectivity', fontsize=8, ha='left', va='bottom', transform=ax.transAxes ) elif col % 3 == 2: z = blended[time_index].values y = blended[time_index].y x = blended[time_index].x if row == 0: ax.text( 0, 1, 'Blended Reflectivity', fontsize=8, ha='left', va='bottom', transform=ax.transAxes ) if col == 0: ax.text( -0.25, 0, f"Initial @ {basetime.strftime('%H:%MZ')} \n" + f"Forecast Valid @ {timelabel.strftime('%H:%MZ')} ", fontsize=8, ha='left', va='bottom', transform=ax.transAxes, rotation=90 ) # plot base map plot_base(ax) # plot reflectivity ax.contourf( x, y, z, 60, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm, levels=levels ) 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_3hr.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_vn_001.png :alt: nwp radar qpf vn :srcset: /auto_examples/images/sphx_glr_nwp-radar_qpf_vn_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 345-348 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 348-363 .. 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-22 04:01:34.372580 Initialising time: 2024-04-22 04:01:35.322663 NWP bias correction time: 2024-04-22 04:01:36.920984 SLA time: 2024-04-22 04:01:54.221705 RaINS time: 2024-04-22 04:01:54.385413 Plotting radar image time: 2024-04-22 04:02:02.891491 Time to initialise: 0 days 00:00:00.950083 Time to run NWP bias correction: 0 days 00:00:01.598321 Time to run rover: 0 days 00:00:00.110280 Time to perform SLA: 0 days 00:00:17.190441 Time to perform RaINS: 0 days 00:00:00.163708 Time to plot radar image: 0 days 00:00:08.506078 .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 29.178 seconds) .. _sphx_glr_download_auto_examples_nwp-radar_qpf_vn.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_vn.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: nwp-radar_qpf_vn.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_