{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# STEPS (Hong Kong)\nThis example demonstrates how to make use of STEPS to perform 3-hour rainfall\nnowcasting with radar data in Hong Kong\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\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": [ "import os\nimport textwrap\n\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\n\nfrom pyresample import utils\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\nimport cartopy.io.shapereader as shpreader\nfrom matplotlib.gridspec import GridSpec\nfrom matplotlib.colors import BoundaryNorm, LinearSegmentedColormap, ListedColormap\nfrom matplotlib.cm import ScalarMappable\n\nfrom swirlspy.rad.iris import read_iris_grid\nfrom swirlspy.qpe.utils import locate_file, timestamps_ending\nfrom swirlspy.core.resample import grid_resample\n\nfrom swirlspy.utils import FrameType\nfrom swirlspy.utils import standardize_attr, FrameType\nfrom swirlspy.utils import conversion\nfrom swirlspy.qpf import steps\nfrom swirlspy.qpf import dense_lucaskanade" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the working directory and nowcast parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# working_dir = os.path.join(os.getcwd(), 'swirlspy/examples')\nworking_dir = os.getcwd()\nradar_dir = os.path.abspath(\n os.path.join(working_dir, '../tests/samples/iris/ppi')\n)\n\n# Set nowcast parameters\nn_timesteps = int(3 * 60 / 6) # 3 hours, each timestamp is 6 minutes\nn_ens_members = 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the User Grid\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "area_id = \"hk1980_250km\"\ndescription = (\"A 250 m resolution rectangular grid \"\n \"centred at HKO and extending to 250 km \"\n \"in each direction in HK1980 easting/northing coordinates\")\nproj_id = 'hk1980'\nprojection = ('+proj=tmerc +lat_0=22.31213333333334 '\n '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '\n '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'\n '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '\n '+no_defs')\nx_size = 500\ny_size = 500\narea_extent = (587000, 569000, 1087000, 1069000)\narea_def_tgt = utils.get_area_def(\n area_id, description, proj_id, projection, x_size, y_size, area_extent\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the plotting function:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining plot parameters\nmap_shape_file = os.path.abspath(os.path.join(\n working_dir,\n '../tests/samples/shape/hk'\n))\n\n# coastline and province\nmap_with_province = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n)\n\n\ndef plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection):\n ax.set_extent(extents, crs=crs)\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, province and state, color\n ax.add_feature(map_with_province,\n facecolor=cfeature.COLORS['land'], edgecolor='none', zorder=0)\n # overlay coastline, province and state without color\n ax.add_feature(map_with_province, facecolor='none',\n edgecolor='gray', linewidth=0.5)\n\n ax.set_title('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading radar data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Log the start time\nstart_time = pd.Timestamp.now()\n\n# Define the basetime\nbasetime = pd.Timestamp('201902190800')\n\n# Generate the timestamps and locate the files\nlocated_files = []\nradar_ts = timestamps_ending(\n basetime,\n duration=pd.Timedelta(minutes=12),\n exclude_end=False\n)\nfor timestamp in radar_ts:\n located_files.append(locate_file(radar_dir, timestamp))\n\n# Read in the data\nreflectivity_list = [] # stores reflec from read_iris_grid()\nfor filename in located_files:\n reflec = read_iris_grid(filename)\n reflectivity_list.append(reflec)\n\n# Reproject the radar data to the user-defined grid\narea_def_src = reflectivity_list[0].attrs['area_def']\nreproj_reflectivity_list = []\nfor reflec in reflectivity_list:\n reproj_reflec = grid_resample(\n reflec, area_def_src, area_def_tgt,\n coord_label=['x', 'y']\n )\n reproj_reflectivity_list.append(reproj_reflec)\n\n# Fill in all fields of the xarray of reflectivity data\nframes = xr.concat(reproj_reflectivity_list,\n dim='time').sortby(['y'], ascending=False)\nstandardize_attr(frames, frame_type=FrameType.dBZ)\n\n\n# Convert from reflectivity to rainfall rain\nframes = conversion.to_rainfall_rate(frames, True, a=58.53, b=1.56)\n\n# Set the fill value\nframes.attrs['zero_value'] = -15.0\n\n# apply threshold to -10dBR i.e. 0.1mm/h\nthreshold = -10.0\nframes.values[frames.values < threshold] = frames.attrs['zero_value']\n\n# Set missing values with the fill value\nframes.values[~np.isfinite(frames.values)] = frames.attrs['zero_value']\n\n# Log the time for record\ninitialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running Lucas Kanade Optical flow and S-PROG\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Estimate the motion field with Lucas Kanade\nmotion = dense_lucaskanade(frames)\n\n# Log the time for record\nmotion_time = pd.Timestamp.now()\n\n# Nowcast using STEP\nforcast_frames = steps(\n frames,\n motion,\n n_timesteps,\n n_ens_members=n_ens_members,\n n_cascade_levels=8,\n R_thr=threshold,\n kmperpixel=2.0,\n decomp_method=\"fft\",\n bandpass_filter_method=\"gaussian\",\n noise_method=\"nonparametric\",\n probmatching_method=\"mean\",\n vel_pert_method=\"bps\",\n mask_method=\"incremental\",\n seed=24\n)\n\nsteps_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining the colour scale\nlevels = [\n -32768,\n 10, 15, 20, 24, 28, 32,\n 34, 38, 41, 44, 47, 50,\n 53, 56, 58, 60, 62\n]\ncmap = ListedColormap([\n '#FFFFFF00', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n])\n\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\nmappable = ScalarMappable(cmap=cmap, norm=norm)\nmappable.set_array([])\n\n# Defining the crs\ncrs = area_def_tgt.to_cartopy_crs()\n\n# Defining area\nx = frames.coords['x'].values\ny = frames.coords['y'].values\nextents = [\n x[0], y[0],\n x[-1], y[-1]\n]\n\n# Generate a time steps for every hour\ntime_steps = [\n (basetime + pd.Timedelta(minutes=6*i)) \n for i in range(n_timesteps + 1) if i % 10 == 0\n]\n\nref_frames = conversion.to_reflectivity(forcast_frames, True)\nref_frames.data[ref_frames.data < 0.1] = np.nan\n\nqx = motion.coords['x'].values[::25]\nqy = motion.coords['y'].values[::25]\nqu = motion.values[0, ::25, ::25]\nqv = motion.values[1, ::25, ::25]\n\nn_rows = len(time_steps)\nfig: plt.Figure = plt.figure(\n figsize=(n_ens_members * 3.5 + 1, n_rows * 3.5), frameon=False)\ngs = GridSpec(n_rows, n_ens_members, figure=fig,\n wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845)\n\nfor row in range(n_rows):\n for col in range(n_ens_members):\n ax: plt.Axes = fig.add_subplot(gs[row, col], projection=crs)\n\n ensemble = ref_frames.coords['ensembles'].values[col]\n t = time_steps[row]\n\n # plot base map\n plot_base(ax, extents, crs)\n\n # plot reflectivity\n member = ref_frames.sel(ensembles=ensemble)\n frame = member.sel(time=t)\n im = ax.imshow(frame.values, cmap=cmap, norm=norm, interpolation='nearest',\n extent=extents)\n\n # plot motion vector\n ax.quiver(\n qx, qy, qu, qv, pivot='mid'\n )\n\n ax.text(\n extents[0],\n extents[1],\n textwrap.dedent(\n \"\"\"\n Reflectivity\n Based @ {baseTime}\n \"\"\"\n ).format(\n baseTime=basetime.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='left',\n linespacing=1\n )\n ax.text(\n extents[2] - (extents[2] - extents[0]) * 0.03,\n extents[1],\n textwrap.dedent(\n \"\"\"\n {validDate}\n Valid @ {validTime}\n \"\"\"\n ).format(\n validDate=basetime.strftime('%Y-%m-%d'),\n validTime=t.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='right',\n linespacing=1\n )\n\ncbar_ax = fig.add_axes([0.875, 0.075, 0.03, 0.85])\ncbar = fig.colorbar(\n mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g')\ncbar.ax.set_ylabel(ref_frames.attrs['values_name'], rotation=90)\n\n\nfig.savefig(\n os.path.join(\n working_dir,\n \"../tests/outputs/steps-reflectivity.png\"\n ),\n bbox_inches='tight'\n)\n\nradar_image_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accumulating hourly rainfall for 3-hour forecast\n\nHourly accumulated rainfall is calculated every 30 minutes, the first\nendtime is the basetime i.e. T+30min.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Convert from rainfall rate to rainfall\nrf_frames = conversion.to_rainfall_depth(ref_frames, a=58.53, b=1.56)\n\n# Compute hourly accumulated rainfall every 30 minutes.\nacc_rf_frames = []\nfor ens in rf_frames.coords['ensembles']:\n af = conversion.acc_rainfall_depth(\n rf_frames.sel(ensembles=ens).drop('ensembles'), basetime +\n pd.Timedelta(minutes=60), basetime + pd.Timedelta(hours=3)\n )\n af_ensembles = af.assign_coords(ensembles=ens)\n acc_rf_frames.append(af_ensembles.expand_dims('ensembles'))\nacc_rf_frames = xr.concat(acc_rf_frames, dim='ensembles')\n\n# Replace zero value with NaN\nacc_rf_frames.data[acc_rf_frames.data <=\n acc_rf_frames.attrs['zero_value']] = np.nan\n\nacc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining colour scale\nlevels = [\n 0, 0.5, 2, 5, 10, 20,\n 30, 40, 50, 70, 100, 150,\n 200, 300, 400, 500, 600, 700\n]\ncmap = ListedColormap([\n '#FFFFFF00', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646',\n '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063',\n '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd'\n])\n\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\nmappable = ScalarMappable(cmap=cmap, norm=norm)\nmappable.set_array([])\n\ntime_steps = acc_rf_frames.coords['time'].values\nn_rows = len(time_steps)\n\nfig: plt.Figure = plt.figure(\n figsize=(n_ens_members * 4 + 1, n_rows * 4), frameon=False)\ngs = GridSpec(n_rows, n_ens_members, figure=fig,\n wspace=0.03, hspace=-0.25, top=0.95, bottom=0.05, left=0.17, right=0.845)\n\nfor row in range(n_rows):\n for col in range(n_ens_members):\n ax = fig.add_subplot(gs[row, col], projection=crs)\n\n ensemble = acc_rf_frames.coords['ensembles'].values[col]\n t = time_steps[row]\n\n # plot base map\n plot_base(ax, extents, crs)\n\n # plot accumulated rainfall depth\n t = pd.Timestamp(t)\n member = acc_rf_frames.sel(ensembles=ensemble)\n frame = member.sel(time=t)\n im = ax.imshow(frame.values, cmap=cmap, norm=norm, interpolation='nearest',\n extent=extents)\n\n ax.text(\n extents[0],\n extents[1],\n textwrap.dedent(\n \"\"\"\n Hourly Rainfall\n Based @ {baseTime}\n \"\"\"\n ).format(\n baseTime=basetime.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='left',\n linespacing=1\n )\n ax.text(\n extents[2] - (extents[2] - extents[0]) * 0.03,\n extents[1],\n textwrap.dedent(\n \"\"\"\n {validDate}\n Valid @ {validTime}\n \"\"\"\n ).format(\n validDate=basetime.strftime('%Y-%m-%d'),\n validTime=t.strftime('%H:%MH')\n ).strip(),\n fontsize=10,\n va='bottom',\n ha='right',\n linespacing=1\n )\n\ncbar_ax = fig.add_axes([0.875, 0.095, 0.03, 0.8])\ncbar = fig.colorbar(\n mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g')\ncbar.ax.set_ylabel(acc_rf_frames.attrs['values_name'], rotation=90)\n\nfig.savefig(\n os.path.join(\n working_dir,\n \"../tests/outputs/steps-rainfall.png\"\n ),\n bbox_inches='tight'\n)\n\nptime = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking run time of each component\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\"Motion field time: {motion_time}\")\nprint(f\"STEPS time: {steps_time}\")\nprint(f\"Plotting radar image time: {radar_image_time}\")\nprint(f\"Accumulating rainfall time: {acc_time}\")\nprint(f\"Plotting rainfall maps: {ptime}\")\n\nprint(f\"Time to initialise: {initialising_time - start_time}\")\nprint(f\"Time to run motion field: {motion_time - initialising_time}\")\nprint(f\"Time to perform STEPS: {steps_time - motion_time}\")\nprint(f\"Time to plot radar image: {radar_image_time - steps_time}\")\nprint(f\"Time to accumulate rainfall: {acc_time - radar_image_time}\")\nprint(f\"Time to plot rainfall maps: {ptime - acc_time}\")\n\nprint(f\"Total: {ptime - start_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 }