{ "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\nfrom pyresample import utils\nimport xarray\nimport cartopy.feature as cfeature\nimport cartopy.crs as ccrs\nimport matplotlib\nimport imageio\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import BoundaryNorm\n\nfrom swirlspy.rad.irisref import read_irisref\nfrom swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending\nfrom swirlspy.qpe.utils import multiple_acc\nfrom swirlspy.obs.rain import Rain\nfrom swirlspy.qpf.rover import rover\nfrom swirlspy.qpf.sla import sla\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(minutes=60)\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_irisref().\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflectivity_list = [] # stores reflec from read_iris()\nfor filename in located_files:\n reflec = read_irisref(\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 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": [ "Step 5: Reproject the radar data from read_irisref() from Centered Azimuthal\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": [ "xarray1 = reproj_reflectivity_list[-2]\nxarray2 = reproj_reflectivity_list[-1]\n\ninitialising_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 = xarray.concat([xarray1, xarray2], dim='time')\n\n# Rover\nmotion_u, motion_v = rover(\n reflec_concat\n)\n\nrover_time = pd.Timestamp.now()\n\n# Semi Lagrangian Advection\nreflectivity = sla(\n reflec_concat, motion_u, motion_v, steps=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": [ "reproj_reflectivity_list.append(reflectivity[1:, ...])\nreflectivity = xarray.concat(reproj_reflectivity_list, dim='time')\n\nconcat_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\n#. Define the colour scale and the format of the\n reflectivity plots.\n#. Initialise the figure and the cartopy GeoAxes. Use the\n to_cartopy_crs() method of the AreaDefinition\n to obtain the projection system of the axes.\n#. Add GSHHS coastlines using the add_feature() method\n of cartopy GeoAxes.\n#. Plot data using xarray.plot() and make desired adjustments\n to labels.\n#. Plot quiver using the quiver() method of cartopy GeoAxes.\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 = matplotlib.colors.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# generate timesteps for plotting\ntimesteps = [basetime+pd.Timedelta(minutes=60*i-6) for i in range(4)]\n\n# Defining quiver parameters\nqx = motion_u.coords['easting'].values[::5]\nqy = motion_u.coords['northing'].values[::5]\nqu = motion_u.values[::5, ::5]\nqv = motion_v.values[::5, ::5]\n\n# Defining coastlines\nhires = cfeature.GSHHSFeature(\n scale='h',\n levels=[1],\n edgecolor='black',\n facecolor='none'\n )\n\n# Plotting\nimages = []\nfor time in timesteps:\n # Figure\n plt.figure(figsize=(28, 21))\n plt.axis(\"equal\")\n # Axes\n crs = area_def_tgt.to_cartopy_crs()\n ax = plt.axes(projection=crs)\n # Coastlines\n ax.add_feature(hires)\n\n # Plot reflectivity data\n quadmesh = reflectivity.sel(time=time).plot(\n cmap=cmap,\n norm=norm,\n extend='neither'\n )\n\n # Customising labels\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n reflectivity.attrs['long_name']+'['+reflectivity.attrs['units']+']',\n fontsize=28\n )\n cbar.ax.tick_params(labelsize=24)\n ax.xaxis.set_visible(True)\n ax.yaxis.set_visible(True)\n ax.xaxis.set_tick_params(labelsize=24)\n ax.yaxis.set_tick_params(labelsize=24)\n ax.xaxis.label.set_size(28)\n ax.yaxis.label.set_size(28)\n\n # Quiver\n ax.quiver(\n qx, qy, qu, qv, pivot='mid', regrid_shape=20\n )\n\n # Title\n plt.title(\n \"Quantitative Precipitation Forecast with Reflectivity Fields\\n\"\n f\"Base time: {basetime}H\\n\"\n f\"Valid time {pd.Timestamp(time)}H\",\n fontsize=32\n )\n savepath = os.path.join(\n THIS_DIR,\n (\"../tests/outputs/rover-output-map-\"\n f\"{time.strftime('%Y%m%d%H%M')}.png\")\n )\n plt.savefig(savepath)\n images.append(imageio.imread(savepath))\n\nimageio.mimsave(\n THIS_DIR +\n \"/../tests/outputs/rover-output-map.gif\",\n images, duration=1\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 rainrates in mm/h with dbz2rr().\n#. Changing time coordinates of xarray from start time to endtime.\n#. Convert rainrates to rainfalls in 6 mins with rr2rf().\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\nrainrates = dbz2rr(reflectivity, a=58.53, b=1.56)\n\n# Converting the coordinates of xarray from start to endtime\nrainrates_endtime = rainrates.copy()\nrainrates_endtime.coords['time'] = \\\n [\n pd.Timestamp(t) + pd.Timedelta(minutes=6)\n for t in rainrates_endtime.coords['time'].values\n ]\n\n# Convert rainrates to accumulated rainfalls every 6 minutes with rr2rf().\nrainfalls = rr2rf(rainrates_endtime)\n\n# Accumulate hourly rainfall every 30 minutes\nacc_rf = multiple_acc(\n rainfalls, basetime, basetime+pd.Timedelta(hours=3)\n)\n\nacc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting rainfall maps\n\n1. Define the colour scheme.\n2. Obtaining the crs of the projection system using to_cartopy_crs()\n method of AreaDefinition.\n3. Defining the zoom or area extent. Tuple order is (x0, x1, y0, y1)\n as opposed to pyresample (x0, y0, x1, y1).\n4. Initialise figure.\n5. Initialise cartopy GeoAxes.\n6. Adding area extent defined in (3) to axes.\n7. Adding coastlines (this example uses GSHHS)\n8. Plot xarray using xarray.plot() and make desired adjustments to labels.\n9. Adding title (if required).\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 = matplotlib.colors.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 = xarray1.proj_site\nzoom = (\n proj_site[0]-r, proj_site[0]+r, proj_site[1]-r, proj_site[1]+r\n)\n\nimages = []\nfor time in acc_rf.coords['time'].values:\n # Define figure\n plt.figure(figsize=(28, 21))\n plt.axis('equal')\n # Plotting axes\n ax = plt.axes(projection=crs)\n # setting extent\n ax.set_extent(zoom, crs=crs)\n # Adding coastlines\n ax.add_feature(hires)\n\n # Plot data\n quadmesh = acc_rf.sel(time=time).plot(\n cmap=cmap,\n norm=norm,\n extend='neither'\n )\n\n # Customizing labels\n # Customising labels\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n acc_rf.attrs['long_name']+'['+acc_rf.attrs['units']+']',\n fontsize=28\n )\n cbar.ax.tick_params(labelsize=24)\n ax.xaxis.set_visible(True)\n ax.yaxis.set_visible(True)\n ax.xaxis.set_tick_params(labelsize=24)\n ax.yaxis.set_tick_params(labelsize=24)\n ax.xaxis.label.set_size(28)\n ax.yaxis.label.set_size(28)\n\n # Add title\n t_minus = pd.Timestamp(time) - \\\n basetime-pd.Timedelta(minutes=60)\n t_minus = round(t_minus.total_seconds()) // 60\n t_plus = pd.Timestamp(time) - basetime\n t_plus = round(t_plus.total_seconds()) // 60\n plt.title(\n \"Quantitative Precipitation Forecast with Hourly\"\n \" Accumulated Rainfall.\\n\"\n f\"Basetime: {basetime} Valid time: {pd.Timestamp(time)}\\n\"\n \"Start time: t\"\n f\"{t_minus:+} minutes End time: t{t_plus:+} minutes\",\n fontsize=32\n )\n savepath = THIS_DIR + \\\n (\"/../tests/outputs/forecast-\" +\n f\"{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png\")\n plt.savefig(savepath)\n images.append(imageio.imread(savepath))\n\nimageio.mimsave(\n THIS_DIR +\n \"/../tests/outputs/forecast.gif\",\n images, duration=1\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 radar stations 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 radar 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 }