{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \nRadar Advection with NWP Blend (Vietnam) \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.\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, 'vnmha')\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('202001240300')" ] }, { "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(19):\n t = basetime - pd.Timedelta(minutes=i * interval)\n # Radar data nomenclature\n filename = os.path.join(\n data_dir,\n 'uf_vietnam',\n t.strftime(\"composite_uf_vn_%Y%m%d%H%M.nc\")\n )\n reflec = xr.open_dataset(filename).assign_coords(time=[t])\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 '__xarray_dataarray_variable__', therefore, we extract '__xarray_dataarray_variable__'\nradar = reflec_concat['__xarray_dataarray_variable__']\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 hours AFTER basetime[18] (nowcast)\nnwp_datas = []\nfor i in range(1, 19):\n t = basetime + pd.Timedelta(minutes=i * interval)\n # Radar nwp nomenclature\n filename = os.path.join(\n data_dir,\n 'wrf_dbz_vietnam',\n t.strftime(\"wrfout_d01_%Y-%m-%d_%H_%M_00.nc\")\n )\n reflec = xr.open_dataset(filename).assign_coords(time=[t])\n nwp_datas.append(reflec)\n\n# Concatenating the nwp reflectivity list of data by time\nreflec_concat = xr.concat(nwp_datas, dim='time')\n\n# Extracting the nwp data: The nwp dBZ variable is called 'zwrf', extracting 'zwrf'\nnwp = reflec_concat['zwrf'].rename({'lon': 'x', 'lat': 'y'})\nnwp.sortby('y', ascending=False).sortby('time', ascending=True)\n# Remove \"invalid\" data\nnwp.values[nwp.values < 0] = np.nan\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_corrected = nwp_bias_correction(\n radar, nwp, proj4_str='+proj=longlat +datum=WGS84 +no_defs'\n)\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, 18) # 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 = rains(nwp_corrected, swirls[1:]) # skip swirls base frame\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 = [\n blended.x.values[0], blended.x.values[-1],\n blended.y.values[0], blended.y.values[-1]\n]\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_corrected = nwp_corrected.where(nwp_corrected > 5, np.nan)\nblended = blended.where(blended > 5, np.nan)\nswirls = swirls[1:] # skip base time frame" ] }, { "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": [ "plot_height = 3\nplot_width = 2\nn_rows = 2\nn_column = 3\nfig: plt.Figure = plt.figure(\n figsize=(plot_width * n_column + 1, plot_height * n_rows),\n frameon=False\n)\ngs = GridSpec(\n n_rows, n_column, figure=fig,\n wspace=0.05, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nfor row in range(n_rows):\n time_index = (row + 1) * 6 - 1\n timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm')\n\n for col in range(n_column):\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 y = swirls[time_index].y\n x = swirls[time_index].x\n if row == 0:\n ax.text(\n 0, 1,\n 'SWIRLS Reflectivity',\n fontsize=8,\n ha='left',\n va='bottom',\n transform=ax.transAxes\n )\n elif col % 3 == 1:\n z = nwp_corrected[time_index].values\n y = nwp_corrected[time_index].y\n x = nwp_corrected[time_index].x\n if row == 0:\n ax.text(\n 0, 1,\n 'NWP Reflectivity',\n fontsize=8,\n ha='left',\n va='bottom',\n transform=ax.transAxes\n )\n elif col % 3 == 2:\n z = blended[time_index].values\n y = blended[time_index].y\n x = blended[time_index].x\n if row == 0:\n ax.text(\n 0, 1,\n 'Blended Reflectivity',\n fontsize=8,\n ha='left',\n va='bottom',\n transform=ax.transAxes\n )\n\n if col == 0:\n ax.text(\n -0.25, 0,\n f\"Initial @ {basetime.strftime('%H:%MZ')} \\n\" +\n f\"Forecast Valid @ {timelabel.strftime('%H:%MZ')} \",\n fontsize=8,\n ha='left',\n va='bottom',\n transform=ax.transAxes,\n rotation=90\n )\n\n # plot base map\n plot_base(ax)\n\n # plot reflectivity\n ax.contourf(\n x, y, z, 60,\n transform=ccrs.PlateCarree(),\n cmap=cmap, norm=norm, levels=levels\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_3hr.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 }