{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# QPF (Hong Kong)\nThis example demonstrates how to perform\noperational deterministic QPF up to three hours from\nraingauge and radar data, using data from Hong Kong.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\n\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 cartopy.feature as cfeature\nimport cartopy.crs as ccrs\nimport cartopy.io.shapereader as shpreader\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\nfrom pyresample import utils\n\nfrom swirlspy.rad.iris import read_iris_grid\nfrom swirlspy.qpe.utils import locate_file, timestamps_ending\nfrom swirlspy.qpf import rover\nfrom swirlspy.qpf import sla\nfrom swirlspy.utils import standardize_attr, FrameType\nfrom swirlspy.utils.conversion import to_rainfall_depth, acc_rainfall_depth\nfrom swirlspy.core.resample import grid_resample\n\nplt.switch_backend('agg')\n\nTHIS_DIR = os.getcwd()\nos.chdir(THIS_DIR)\n\nstart_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialising\n\nThis section demonstrates extracting\nradar reflectivity data.\n\nStep 1: Define a basetime.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply basetime\nbasetime = pd.Timestamp('201902190800')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Using basetime, generate timestamps of desired radar files\ntimestamps_ending() and locate files using locate_file().\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Obtain radar files\ndir = THIS_DIR + '/../tests/samples/iris/ppi'\nlocated_files = []\nradar_ts = timestamps_ending(\n basetime,\n duration=pd.Timedelta(60, 'm')\n)\n\nfor timestamp in radar_ts:\n located_files.append(locate_file(dir, timestamp))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Read data from radar files into xarray.DataArray\nusing read_iris_grid().\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflectivity_list = [] # stores reflec from read_iris_grid()\nfor filename in located_files:\n reflec = read_iris_grid(\n filename\n )\n reflectivity_list.append(reflec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: Define the target grid as a pyresample AreaDefinition.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining target grid\narea_id = \"hk1980_250km\"\ndescription = (\"A 1km 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": [ "Step 5: Reproject the radar data from read_iris_grid() from Centered\nAzimuthal\n(source) projection to HK 1980 (target) projection.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Extracting the AreaDefinition of the source projection\narea_def_src = reflectivity_list[0].attrs['area_def']\n\n# Reprojecting\nreproj_reflectivity_list = []\nfor reflec in reflectivity_list:\n reproj_reflec = grid_resample(\n reflec, area_def_src, area_def_tgt,\n coord_label=['easting', 'northing']\n )\n reproj_reflectivity_list.append(reproj_reflec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: Assigning reflectivity xarrays at the last two timestamps to\nvariables for use during ROVER QPF.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "initialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running ROVER and Semi-Lagrangian Advection\n\n1. Concatenate two reflectivity xarrays along time dimension.\n2. Run ROVER, with the concatenated xarray as the input.\n3. Perform Semi-Lagrangian Advection using the motion fields from rover.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Combining the two reflectivity DataArrays\n# the order of the coordinate keys is now ['y', 'x', 'time']\n# as opposed to ['time', 'x', 'y']\nreflec_concat = xr.concat(reproj_reflectivity_list, dim='time')\nstandardize_attr(reflec_concat, frame_type=FrameType.dBZ, zero_value=9999.)\n\n# Rover\nmotion = rover(reflec_concat)\n\nrover_time = pd.Timestamp.now()\n\n# Semi Lagrangian Advection\nreflectivity = sla(\n reflec_concat, motion, 30\n)\n\nsla_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Concatenating observed and forecasted reflectivities\n\n1. Add forecasted reflectivity to reproj_reflectivity_list.\n2. Concatenate observed and forecasted reflectivity\n xarray.DataArrays along the time dimension.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflectivity = xr.concat([reflec_concat[:-1, ...], reflectivity], dim='time')\nreflectivity.attrs['long_name'] = 'Reflectivity'\nstandardize_attr(reflectivity)\n\nconcat_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\nDefine the color scale and format of the plots\nand plot using xarray.plot().\n\nIn this example, only hourly images will be plotted.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining colour scale and format\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 '#FFFFFF', '#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\n# Defining the crs\ncrs = area_def_tgt.to_cartopy_crs()\n\n# Generating a timelist for every hour\ntimelist = [\n (basetime + pd.Timedelta(60*i-6, 'm')) for i in range(4)\n]\n\n# Obtaining the slice of the xarray to be plotted\nda_plot = reflectivity.sel(time=timelist)\n\n# Defining motion quivers\nqx = motion.coords['easting'].values[::5]\nqy = motion.coords['northing'].values[::5]\nqu = motion.values[0, ::5, ::5]\nqv = motion.values[1, ::5, ::5]\n\n# Defining coastlines\nmap_shape_file = os.path.join(THIS_DIR, \"./../tests/samples/shape/rsmc\")\nocean_color = np.array([[[178, 208, 254]]], dtype=np.uint8)\nland_color = cfeature.COLORS['land']\ncoastline = cfeature.ShapelyFeature(\n list(shpreader.Reader(map_shape_file).geometries()),\n ccrs.PlateCarree()\n)\n\n# Plotting\np = da_plot.plot(\n col='time', col_wrap=2,\n subplot_kws={'projection': crs},\n cbar_kwargs={\n 'extend': 'max',\n 'ticks': levels[1:],\n 'format': '%.3g'\n },\n cmap=cmap,\n norm=norm\n)\nfor idx, ax in enumerate(p.axes.flat):\n # ocean\n ax.imshow(np.tile(ocean_color, [2, 2, 1]),\n origin='upper',\n transform=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180],\n zorder=-1)\n # coastline, color\n ax.add_feature(coastline,\n facecolor=land_color, edgecolor='none', zorder=0)\n # overlay coastline without color\n ax.add_feature(coastline, facecolor='none',\n edgecolor='gray', linewidth=0.5, zorder=3)\n ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\n ax.gridlines()\n ax.set_title(\n \"Reflectivity\\n\"\n f\"Based @ {basetime.strftime('%H:%MH')}\",\n loc='left',\n fontsize=9\n )\n ax.set_title(\n ''\n )\n ax.set_title(\n f\"{basetime.strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {timelist[idx].strftime('%H:%MH')} \",\n loc='right',\n fontsize=9\n )\nplt.savefig(\n THIS_DIR +\n f\"/../tests/outputs/rover-output-map-hk.png\",\n dpi=300\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+0min.\n\n#. Convert reflectivity in dBZ to rainfalls in 6 mins with to_rainfall_depth().\n#. Changing time coordinates of xarray from start time to endtime.\n#. Accumulate hourly rainfall every 30 minutes using multiple_acc().\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Convert reflectivity to rainrates\nrainfalls = to_rainfall_depth(reflectivity, a=58.53, b=1.56)\n\n# Converting the coordinates of xarray from start to endtime\nrainfalls.coords['time'] = [\n pd.Timestamp(t) + pd.Timedelta(6, 'm')\n for t in rainfalls.coords['time'].values\n]\n\n# Accumulate hourly rainfall every 30 minutes\nacc_rf = acc_rainfall_depth(\n rainfalls,\n basetime,\n basetime + pd.Timedelta(hours=3)\n)\nacc_rf.attrs['long_name'] = 'Rainfall accumulated over the past 60 minutes'\n\nacc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting rainfall maps\n\nDefine the colour scheme and format and plot using xarray.plot().\n\nIn this example, only hourly images will be plotted.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Defining the colour scheme\nlevels = [\n 0, 0.5, 2, 5, 10, 20,\n 30, 40, 50, 70, 100, 150,\n 200, 300, 400, 500, 600, 700\n]\n\ncmap = ListedColormap([\n '#ffffff', '#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\n# Defining projection\ncrs = area_def_tgt.to_cartopy_crs()\n# Defining zoom extent\nr = 64000\nproj_site = acc_rf.proj_site\nzoom = (\n proj_site[0]-r, proj_site[0]+r, proj_site[1]-r, proj_site[1]+r\n) # (x0, x1, y0, y1)\n\n# Defining times for plotting\ntimelist = [basetime + pd.Timedelta(i, 'h') for i in range(4)]\n\n# Obtaining xarray slice to be plotted\nda_plot = acc_rf.sel(\n easting=slice(zoom[0], zoom[1]),\n northing=slice(zoom[3], zoom[2]),\n time=timelist\n)\n\n# Plotting\np = da_plot.plot(\n col='time', col_wrap=2,\n subplot_kws={'projection': crs},\n cbar_kwargs={\n 'extend': 'max',\n 'ticks': levels,\n 'format': '%.3g'\n },\n cmap=cmap,\n norm=norm\n)\nfor idx, ax in enumerate(p.axes.flat):\n # ocean\n ax.imshow(np.tile(ocean_color, [2, 2, 1]),\n origin='upper',\n transform=ccrs.PlateCarree(),\n extent=[-180, 180, -180, 180],\n zorder=-1)\n # coastline, color\n ax.add_feature(coastline,\n facecolor=land_color, edgecolor='none', zorder=0)\n # overlay coastline without color\n ax.add_feature(coastline, facecolor='none',\n edgecolor='gray', linewidth=0.5, zorder=3)\n ax.gridlines()\n ax.set_xlim(zoom[0], zoom[1])\n ax.set_ylim(zoom[2], zoom[3])\n ax.set_title(\n \"Past Hour Rainfall\\n\"\n f\"Based @ {basetime.strftime('%H:%MH')}\",\n loc='left',\n fontsize=8\n )\n ax.set_title(\n ''\n )\n ax.set_title(\n f\"{basetime.strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {timelist[idx].strftime('%H:%MH')} \",\n loc='right',\n fontsize=8\n )\nplt.savefig(\n THIS_DIR +\n f\"/../tests/outputs/rainfall_hk.png\",\n dpi=300\n)\n\nrf_image_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extract the rainfall values at a specified location\n\nIn this example, the rainfall values at the location is assumed\nto be the same as the nearest gridpoint.\n\n1. Read information regarding the rain gauge stations\n into a pandas.DataFrame.\n2. Extract the rainfall values at the nearest gridpoint to location\n for given times (in this example, 30 minute intervals).\n3. Store rainfall values over time in a pandas.DataFrame.\n4. Plot the time series of rainfall at different stations.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Getting rain gauge station coordinates\ndf = pd.read_csv(\n os.path.join(THIS_DIR, \"../tests/samples/hk_raingauge.csv\"),\n usecols=[0, 1, 2, 3, 4]\n)\n\n# Extract rainfall values at gridpoint closest to the\n# location specified for given timesteps and storing it\n# in pandas.DataFrame.\n\nrf_time = []\nfor time in acc_rf.coords['time'].values:\n rf = []\n for index, row in df.iterrows():\n rf.append(acc_rf.sel(\n time=time, northing=row[1],\n easting=row[2],\n method='nearest'\n ).values)\n rf_time.append(rf)\n\nrf_time = np.array(rf_time)\n\nstation_rf = pd.DataFrame(\n data=rf_time,\n columns=df.iloc[:, 0],\n index=pd.Index(\n acc_rf.coords['time'].values,\n name='time'\n )\n)\nprint(station_rf)\n\n# Plotting time series graph\nax = station_rf.plot(title=\"Time Series of Hourly Accumulated Rainfall\",\n grid=True)\nax.set_ylabel(\"Hourly Accumulated Rainfall (mm)\")\nplt.savefig(THIS_DIR+\"/../tests/outputs/qpf_time_series.png\")\n\nextract_time = 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\"Rover time: {rover_time}\")\nprint(f\"SLA time: {sla_time}\")\nprint(f\"Concatenating time: {concat_time}\")\nprint(f\"Plotting radar image time: {radar_image_time}\")\nprint(f\"Accumulating rainfall time: {acc_time}\")\nprint(f\"Plotting rainfall map time: {rf_image_time}\")\nprint(f\"Extracting and plotting time series time: {extract_time}\")\n\nprint(f\"Time to initialise: {initialising_time-start_time}\")\nprint(f\"Time to run rover: {rover_time-initialising_time}\")\nprint(f\"Time to perform SLA: {sla_time-rover_time}\")\nprint(f\"Time to concatenate xarrays: {concat_time - sla_time}\")\nprint(f\"Time to plot radar image: {radar_image_time - concat_time}\")\nprint(f\"Time to accumulate rainfall: {acc_time - radar_image_time}\")\nprint(f\"Time to plot rainfall maps: {rf_image_time-acc_time}\")\nprint(f\"Time to extract and plot time series: {extract_time-rf_image_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 }