{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nSPROG (Hong Kong)\n========================================================\nThis example demonstrates how to use SPROG to forecast rainfall up to\nthree hours, using rain guage and radar data from Hong Kong.\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 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 sprog, 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 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" ] }, { "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 = (\n \"A 250 m resolution rectangular grid centred at HKO and extending\"\n \"to 250 km in each direction in HK1980 easting/northing coordinates\"\n)\nproj_id = 'hk1980'\nprojection = (\n '+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 +no_defs'\n)\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 base map:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Load the shape of Hong Kong\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\n# define the plot function\ndef plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection):\n ax.set_extent(extents, crs=crs)\n # fake the ocean color\n ax.imshow(\n np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]),\n origin='upper', transform=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180], zorder=-1\n )\n # coastline, province and state, color\n ax.add_feature(\n map_with_province, facecolor=cfeature.COLORS['land'],\n edgecolor='none', zorder=0\n )\n # overlay coastline, province and state without color\n ax.add_feature(\n map_with_province, facecolor='none', edgecolor='gray', linewidth=0.5\n )\n ax.set_title('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Log the start time for reference:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "start_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Loading Radar Data\n---------------------------------------------------\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Specify the basetime\nbasetime = pd.Timestamp('201902190800')\n\n# Generate a list of timestamps of the radar data files\nlocated_files = []\nradar_ts = timestamps_ending(\n basetime,\n duration=pd.Timedelta(minutes=60),\n exclude_end=False\n)\nfor timestamp in radar_ts:\n located_files.append(locate_file(radar_dir, timestamp))\n\n# Read in the radar 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# Standardize reflectivity xarrays\nraw_frames = xr.concat(reproj_reflectivity_list,\n dim='time').sortby(['y'], ascending=False)\nstandardize_attr(raw_frames, frame_type=FrameType.dBZ)\n\n# Transform from reflecitiy to rainfall rate\nframes = conversion.to_rainfall_rate(raw_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\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\nmotion = dense_lucaskanade(frames)\n\nmotion_time = pd.Timestamp.now()\n\n# Generate forecast rainrate field\nforcast_frames = sprog(\n frames,\n motion,\n n_timesteps,\n n_cascade_levels=8,\n R_thr=threshold,\n decomp_method=\"fft\",\n bandpass_filter_method=\"gaussian\",\n probmatching_method=\"mean\",\n)\n\nsprog_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating radar reflectivity maps\n----------------------------------\n\nDefine the color scale and format of the plot.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining colour scale and format.\nlevels = [\n -32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 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\nx_d = x[1] - x[0]\ny_d = y[1] - y[0]\nextents = [x[0], y[0], x[-1], y[-1]]\n\n# Generating 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\nref_frames = xr.concat([raw_frames[:-1, ...], ref_frames], dim='time')\nref_frames.attrs['values_name'] = 'Reflectivity 2km CAPPI'\nstandardize_attr(ref_frames)\n\nqx = motion.coords['x'].values[::5]\nqy = motion.coords['y'].values[::5]\nqu = motion.values[0, ::5, ::5]\nqv = motion.values[1, ::5, ::5]\n\nfig: plt.Figure = plt.figure(figsize=(8, 8), frameon=False)\ngs = GridSpec(\n 2, 2, figure=fig, wspace=0.03, hspace=-0.25, top=0.95,\n bottom=0.05, left=0.17, right=0.845\n)\n\nfor i, t in enumerate(time_steps):\n row = i // 2\n col = i % 2\n ax = fig.add_subplot(gs[row, col], projection=crs)\n\n # plot base map\n plot_base(ax, extents, crs)\n\n # plot reflectivity\n frame = ref_frames.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(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\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.125, 0.03, 0.75])\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\nfig.savefig(\n os.path.join(\n OUTPUT_DIR,\n \"sprog-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": [ "# Optional, convert to rainfall depth\nrf_frames = conversion.to_rainfall_depth(ref_frames, a=58.53, b=1.56)\n\n# Compute hourly accumulated rainfall every 60 minutes.\nacc_rf_frames = conversion.acc_rainfall_depth(\n rf_frames,\n basetime,\n basetime + pd.Timedelta(hours=3),\n pd.Timedelta(minutes=60)\n)\n\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 and format.\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\nfig: plt.Figure = plt.figure(figsize=(8, 8), frameon=False)\ngs = GridSpec(\n 2, 2, figure=fig, wspace=0.03, hspace=-0.25, top=0.95,\n bottom=0.05, left=0.17, right=0.845\n)\n\nfor i, t in enumerate(acc_rf_frames.coords['time'].values):\n row = i // 2\n col = i % 2\n ax = fig.add_subplot(gs[row, col], projection=crs)\n\n # plot base map\n plot_base(ax, extents, crs)\n\n # plot accumulated rainfall depth\n t = pd.Timestamp(t)\n frame = acc_rf_frames.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.125, 0.03, 0.75])\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 \"sprog-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\"S-PROG time: {sprog_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 S-PROG: {sprog_time - motion_time}\")\nprint(f\"Time to plot radar image: {radar_image_time - sprog_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 }