{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# SPROG (Hong Kong)\nThis example demonstrates how to use SPROG to forecast rainfall up to\nthree hours, using rain guage and radar data from Hong Kong.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\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": [ "import os\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\nimport textwrap\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\nfrom matplotlib.colors import 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, standardize_attr, FrameType, conversion\nfrom swirlspy.qpf import sprog, dense_lucaskanade" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define 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" ] }, { "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 = 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 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 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\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" ] }, { "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" ] }, { "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\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 working_dir,\n \"../tests/outputs/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\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" ] }, { "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 working_dir,\n \"../tests/outputs/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" ] }, { "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 }