{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# QPF (Laos)\nThis example demonstrates how to perform\noperational deterministic QPF up to three hours using \nreflectivity data.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Definitions\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# 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 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 package to allow system command line functions\nimport os\n\n# working directory\nworking_dir = os.getcwd()\nos.chdir(working_dir)\n\n# output directory\noutput_dir = os.path.abspath(os.path.join(working_dir, '../tests/outputs'))\n\n# data files\ndata_paths = [\n os.path.abspath(os.path.join(working_dir, '../tests/samples/laos_h8/laos_20190731070000.nc')),\n os.path.abspath(os.path.join(working_dir, '../tests/samples/laos_h8/laos_20190731071000.nc')),\n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Initialising\n\nThis section demonstrates preprocessing\nreflectivity data.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Read data from saved files\nreflec_datas = []\nfor file_path in data_paths:\n d = xr.open_dataarray(file_path)\n reflec_datas.append(d)\n\n# Concatenate list on time axis\ndata_frames = xr.concat(reflec_datas, dim='time')\n\n# Make sure data is ordered by time\ndata_frames = data_frames.sortby('time', ascending=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Nowcast (SWIRLS-Radar-Advection)\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 <numpy.nan> and frame_type <FrameType.dBZ>\nstandardize_attr(data_frames)\n\n# Rover motion field computation\nmotion = rover(data_frames)\n\nrover_time = pd.Timestamp.now()\n\n# Semi-Lagrangian Advection\nswirls = sla(data_frames, motion, 17) # Radar time goes from earliest to latest"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting result\nStep 1: Import plotting library and necessary library\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# 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\nplt.switch_backend('agg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Step 2: Defining the dBZ levels, colorbar parameters and projection\n\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 '#FFFFFF00', '#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 working_dir,\n '../tests/samples/shape/se_asia'\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.5, 108., 13.75, 22.75]\n\n# base_map plotting function\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 3: Filtering\nvalues <= 15dbZ are not plotted\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"swirls = swirls.where(swirls > 15, np.nan)\n\n# Defining motion quivers\nqx = motion.coords['x'].values[::5]\nqy = motion.coords['y'].values[::5]\nqu = motion.values[0, ::5, ::5]\nqv = motion.values[1, ::5, ::5]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Step 4: 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 * 3.5 + 1, 3),\n frameon=False\n)\ngs = GridSpec(\n 1, 3, figure=fig,\n wspace=0, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nbasetime = pd.Timestamp(swirls.time.values[0])\ninterval = swirls.attrs['step_size']\n\nfor col in range(3):\n time_index = (col + 1) * 6 - 1\n timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm')\n\n ax: plt.Axes = fig.add_subplot(\n gs[0, col],\n projection=ccrs.PlateCarree()\n )\n z = swirls[time_index].values\n y = swirls[time_index].y\n x = swirls[time_index].x\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\n # plot motion field\n ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\n\n ax.set_title(\n f\"Nowcast\\n\" +\n f\"Initial @ {basetime.strftime('%H:%MZ')}\",\n loc='left', fontsize=8.75\n )\n ax.set_title('')\n ax.set_title(\n f\"Initial {basetime.strftime('%Y-%m-%d')} \\n\" +\n f\"Valid @ {timelabel.strftime('%H:%MZ')} \",\n loc='right', fontsize=8.75\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 \"qpf_laos.png\"\n ),\n dpi=450,\n bbox_inches=\"tight\",\n pad_inches=0.25\n)\n\nradar_image_time = pd.Timestamp.now()"
]
}
],
"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
}