{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \nRadar Advection with NWP Blend (Malaysia) \n===============================================================\nThis example shows the technique of blending Numerical Weather \nForecast (NWP) with COM-SWIRLS output to generate operational \ndeterministic QPF up to 3 and 6 hours ahead. NWP data was generated\nusing the WRF Model while radar data was taken from the \nMalaysian radar observation network\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Definitions\n--------------------------------------------------------\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import all required modules and methods:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Python package to allow system command line functions\nimport os\n# Python package to manage warning message\nimport warnings\n# Python package for timestamp\nimport pandas as pd\n# Python package for xarrays to read and handle netcdf data\nimport xarray as xr\n# Python package for numerical calculations\nimport numpy as np\n# Python package for reading map shape file\nimport cartopy.io.shapereader as shpreader\n# Python package for land/sea features\nimport cartopy.feature as cfeature\n# Python package for projection\nimport cartopy.crs as ccrs\n# Python package for creating plots\nfrom matplotlib import pyplot as plt\n# Python package for output import grid\nfrom matplotlib.gridspec import GridSpec\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\n# Python package for scalar data to RGBA mapping\nfrom matplotlib.cm import ScalarMappable\n# Python com-swirls package to standardize attributes\nfrom swirlspy.utils import standardize_attr, FrameType\n# Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection\nfrom swirlspy.qpf import rover, sla\n# Python com-swirls package to blend nwp and nowcast (RaINS)\nfrom swirlspy.blending import rains, nwp_bias_correction\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")\nstart_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initializing\n---------------------------------------------------\n\nThis section demonstrates the extraction of radar & nwp data\nfrom netcdf into python\n\nStep 1: Define your input data directory and output directory\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply the directory of radar and nwp data\ndata_dir = os.path.abspath(\n os.path.join(DATA_DIR, 'netcdf_ms')\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Define a basetime\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply basetime\nbasetime = pd.Timestamp('201908090900')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Read data files from the radar data using xarray()\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Radar data listed from the basetime[0] --> 3 hours before the basetime[17] (descending time)\ninterval = 10 # Interval of radar data\nradar_datas = []\nfor i in range(0, 18):\n t = basetime - pd.Timedelta(minutes=i * interval)\n # Radar data nomenclature\n filename = os.path.join(\n data_dir,\n t.strftime(\"radar_d03_%Y-%m-%d_%H_%M_00.rapids.nc\")\n )\n reflec = xr.open_dataset(filename)\n radar_datas.append(reflec)\n\n# Concatenate list by time\nreflec_concat = xr.concat(radar_datas, dim='time')\n\n# Extracting the radar data: The radar dBZ variable is named 'Zradar', therefore, we extract 'Zradar'\nradar = reflec_concat['Zradar']\n\n# Reversing such that time goes from earliest to latest; 3 hours before basetime[0] --> basetime[17]\nradar = radar.sortby('time', ascending=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: Reading nwp netcdf data into xarray\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# NWP data listed from the basetime[0] --> 3 and 6 hours AFTER basetime[17] (nowcast)\nnwp_3hr_datas = []\nnwp_6hr_datas = []\nfor i in range(0, 36):\n t = basetime + pd.Timedelta(minutes=i * interval)\n # Radar nwp nomenclature\n filename = os.path.join(\n data_dir,\n t.strftime(\"wrfout_d03_%Y-%m-%d_%H_%M_00.rapids.nc\")\n )\n reflec = xr.open_dataset(filename)\n nwp_6hr_datas.append(reflec)\n if i < 18:\n nwp_3hr_datas.append(reflec)\n\n# Concatenating the nwp reflectivity list of data by time\nreflec_3hr_concat = xr.concat(nwp_3hr_datas, dim='time')\nreflec_6hr_concat = xr.concat(nwp_6hr_datas, dim='time')\n\n# Extracting the nwp data: The nwp dBZ variable is called 'zwrf', extracting 'zwrf'\nnwp_3hr = reflec_3hr_concat['zwrf']\nnwp_6hr = reflec_6hr_concat['zwrf']\n\ninitialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 5: Bias correction of nwp data\nThe objective of bias correction is to match the nwp percentile to the radar percentile\nThis is also known as frequency matching.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nwp_3hr_corrected = nwp_bias_correction(radar, nwp_3hr)\nnwp_6hr_corrected = nwp_bias_correction(radar, nwp_6hr)\n\nbias_correction_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nowcast (SWIRLS-Radar-Advection)\n---------------------------------------------------\nThe swirls radar advection was performed using the observed radar data\nFirstly, some attributes necessary for com-swirls input variable is added\nSecondly, rover function is invoked to compute the motion field\nThirdly, semi-lagrangian advection is performed to advect the radar data using the rover motion field\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Adding in some attributes that is step_size <10 mins in pandas.Timedelta>, zero_value <9999.> frame_type \nradar.attrs['step_size'] = pd.Timedelta(minutes=10)\nstandardize_attr(radar, frame_type=FrameType.dBZ, zero_value=np.nan)\n\n# Rover motion field computation\nmotion = rover(radar)\n\nrover_time = pd.Timestamp.now()\n\n# Semi-Lagrangian Advection\nswirls = sla(radar, motion, 35) # Radar time goes from earliest to latest\n\nsla_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Blending (RaINS)\n---------------------------------------------------\nBlending Numerical Weather Forecast (NWP) with COM-SWIRLS output to generate operational\ndeterministic QPF up to 3 hours ahead.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "blended_3hr = rains(nwp_3hr_corrected, swirls[0:18])\nblended_6hr = rains(nwp_6hr_corrected, swirls)\n\nrains_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting result\n---------------------------------------------------\nStep 1: Defining the dBZ levels, colorbar parameters and projection\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# levels of colorbar (dBZ)\nlevels = [-32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44,\n 47, 50, 53, 56, 58, 60, 62]\n# hko colormap for dBZ at each levels\ncmap = ListedColormap([\n '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n])\n# boundary\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n# scalar data to RGBA mapping\nscalar_map = ScalarMappable(cmap=cmap, norm=norm)\nscalar_map.set_array([])\n# Defining plot parameters\nmap_shape_file = os.path.abspath(os.path.join(\n DATA_DIR,\n 'shape/se_asia_s'\n))\n# coastline and province\nse_asia = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n)\n# output area\nextents = [99, 120, 0.5, 7.25]\n\n# base_map plotting function\n\n\ndef plot_base(ax: plt.Axes):\n ax.set_extent(extents, crs=ccrs.PlateCarree())\n\n # fake the ocean color\n ax.imshow(np.tile(np.array([[[178, 208, 254]]],\n dtype=np.uint8), [2, 2, 1]),\n origin='upper',\n transform=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180],\n zorder=-1)\n # coastline, state, color\n ax.add_feature(se_asia,\n facecolor=cfeature.COLORS['land'], edgecolor='none', zorder=0)\n # overlay coastline, state without color\n ax.add_feature(se_asia, facecolor='none',\n edgecolor='gray', linewidth=0.5)\n\n ax.set_title('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Filtering\nvalues <= 5dbZ are not plotted\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nwp_3hr_corrected = nwp_3hr_corrected.where(nwp_3hr_corrected > 5, np.nan)\nblended_3hr = blended_3hr.where(blended_3hr > 5, np.nan)\nnwp_6hr_corrected = nwp_6hr_corrected.where(nwp_6hr_corrected > 5, np.nan)\nblended_6hr = blended_6hr.where(blended_6hr > 5, np.nan)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 3 hours ahead\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig: plt.Figure = plt.figure(\n figsize=(3 * 5 + 1, 3 * 2),\n frameon=False\n)\ngs = GridSpec(\n 3, 3, figure=fig,\n wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nfor row in range(3):\n time_index = (row + 1) * 6 - 1\n timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm')\n\n for col in range(3):\n ax: plt.Axes = fig.add_subplot(\n gs[row, col],\n projection=ccrs.PlateCarree()\n )\n if col % 3 == 0:\n z = swirls[time_index].values\n lats = swirls[time_index].latitude\n lons = swirls[time_index].longitude\n title = 'SWIRLS Reflectivity'\n elif col % 3 == 1:\n z = nwp_3hr_corrected[time_index].values\n lats = nwp_3hr_corrected[time_index].lat\n lons = nwp_3hr_corrected[time_index].lon\n title = 'NWP Reflectivity'\n elif col % 3 == 2:\n z = blended_3hr[time_index].values\n lats = blended_3hr[time_index].latitude\n lons = blended_3hr[time_index].longitude\n title = 'Blended Reflectivity'\n\n # plot base map\n plot_base(ax)\n\n # plot reflectivity\n ax.contourf(\n lons, lats, z, 60,\n transform=ccrs.PlateCarree(),\n cmap=cmap, norm=norm, levels=levels\n )\n\n ax.set_title(\n f\"{title}\\n\" +\n f\"Initial @ {basetime.strftime('%H:%MZ')}\",\n loc='left', fontsize=9\n )\n ax.set_title('')\n ax.set_title(\n f\"Initial {basetime.strftime('%Y-%m-%d')} \\n\" +\n f\"Forecast Valid @ {timelabel.strftime('%H:%MZ')} \",\n loc='right', fontsize=9\n )\n\ncbar_ax = fig.add_axes([0.85, 0.105, 0.01, 0.845])\ncbar = fig.colorbar(\n scalar_map, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g'\n)\ncbar.ax.set_ylabel('Reflectivity (dBZ)', rotation=90)\n\nfig.savefig(\n os.path.join(\n OUTPUT_DIR,\n \"swirls_nwp_blend_ms_3hr.png\"\n ),\n dpi=450,\n bbox_inches=\"tight\",\n pad_inches=0.1\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 6 hours ahead\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig: plt.Figure = plt.figure(\n figsize=(3 * 5 + 1, 6 * 2),\n frameon=False\n)\ngs = GridSpec(\n 6, 3, figure=fig,\n wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nfor row in range(6):\n time_index = (row + 1) * 6 - 1\n timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm')\n\n for col in range(3):\n ax: plt.Axes = fig.add_subplot(\n gs[row, col],\n projection=ccrs.PlateCarree()\n )\n if col % 3 == 0:\n z = swirls[time_index].values\n lats = swirls[time_index].latitude\n lons = swirls[time_index].longitude\n title = 'SWIRLS Reflectivity'\n elif col % 3 == 1:\n z = nwp_6hr_corrected[time_index].values\n lats = nwp_6hr_corrected[time_index].lat\n lons = nwp_6hr_corrected[time_index].lon\n title = 'NWP Reflectivity'\n elif col % 3 == 2:\n z = blended_6hr[time_index].values\n lats = blended_6hr[time_index].latitude\n lons = blended_6hr[time_index].longitude\n title = 'Blended Reflectivity'\n\n # plot base map\n plot_base(ax)\n\n # plot reflectivity\n ax.contourf(\n lons, lats, z, 60,\n transform=ccrs.PlateCarree(),\n cmap=cmap, norm=norm, levels=levels\n )\n\n ax.set_title(\n f\"{title}\\n\" +\n f\"Initial @ {basetime.strftime('%H:%MZ')}\",\n loc='left', fontsize=9\n )\n ax.set_title('')\n ax.set_title(\n f\"Initial {basetime.strftime('%Y-%m-%d')} \\n\" +\n f\"Forecast Valid @ {timelabel.strftime('%H:%MZ')} \",\n loc='right', fontsize=9\n )\n\ncbar_ax = fig.add_axes([0.85, 0.105, 0.01, 0.845])\ncbar = fig.colorbar(\n scalar_map, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g'\n)\ncbar.ax.set_ylabel('Reflectivity (dBZ)', rotation=90)\n\nfig.savefig(\n os.path.join(\n OUTPUT_DIR,\n \"swirls_nwp_blend_ms_6hr.png\"\n ),\n dpi=450,\n bbox_inches=\"tight\",\n pad_inches=0.1\n)\n\nradar_image_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Checking run time of each component\n--------------------------------------------------------------------\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(f\"Start time: {start_time}\")\nprint(f\"Initialising time: {initialising_time}\")\nprint(f\"NWP bias correction time: {bias_correction_time}\")\nprint(f\"SLA time: {sla_time}\")\nprint(f\"RaINS time: {rains_time}\")\nprint(f\"Plotting radar image time: {radar_image_time}\")\n\nprint(f\"Time to initialise: {initialising_time - start_time}\")\nprint(\n f\"Time to run NWP bias correction: {bias_correction_time - initialising_time}\")\nprint(f\"Time to run rover: {rover_time - bias_correction_time}\")\nprint(f\"Time to perform SLA: {sla_time - rover_time}\")\nprint(f\"Time to perform RaINS: {rains_time - sla_time}\")\nprint(f\"Time to plot radar image: {radar_image_time - rains_time}\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.15" } }, "nbformat": 4, "nbformat_minor": 0 }