{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nSTEPS (Hong Kong)\n========================================================\nThis example demonstrates how to make use of STEPS to perform 3-hour rainfall\nnowcasting with radar data in Hong Kong\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 time calculations\nimport pandas as pd\n# Python package for numerical calculations\nimport numpy as np\n# Python package for xarrays to read and handle netcdf data\nimport xarray as xr\n# Python package for text formatting\nimport textwrap\n# Python package for projection description\nfrom pyresample import get_area_def\n# Python package for projection\nimport cartopy.crs as ccrs\n# Python package for land/sea features\nimport cartopy.feature as cfeature\n# Python package for reading map shape file\nimport cartopy.io.shapereader as shpreader\n# Python package for creating plots\nfrom matplotlib import pyplot as plt\n# Python package for output grid format\nfrom matplotlib.gridspec import GridSpec\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\nfrom matplotlib.cm import ScalarMappable\n\n# swirlspy data parser function\nfrom swirlspy.rad.iris import read_iris_grid\n# swirlspy test data source locat utils function\nfrom swirlspy.qpe.utils import timestamps_ending, locate_file\n# swirlspy regrid function\nfrom swirlspy.core.resample import grid_resample\n# swirlspy standardize data function\nfrom swirlspy.utils import FrameType, standardize_attr, conversion\n# swirlspy pysteps integrated package\nfrom swirlspy.qpf import steps, dense_lucaskanade\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")" ] }, { "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": [ "radar_dir = os.path.abspath(\n os.path.join(DATA_DIR, '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 = 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 DATA_DIR,\n '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\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\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\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 OUTPUT_DIR,\n \"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------------------------------------------------\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\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 OUTPUT_DIR,\n \"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\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 }