{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# PQPF (Hong Kong)\nThis example demonstrates how to perform\noperational PQPF 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, LinearSegmentedColormap\n\nfrom swirlspy.rad.irisref import read_irisref\nfrom swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending\nfrom swirlspy.qpe.utils import temporal_interp, multiple_acc\nfrom swirlspy.obs.rain import Rain\nfrom swirlspy.qpf.rover import rover\nfrom swirlspy.qpf.sla import sla\nfrom swirlspy.qpf.rover import rover\nfrom swirlspy.core.resample import grid_resample\n\nplt.switch_backend('agg')\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" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflectivity_list = [] # stores reflec from read_irisref()\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 three 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[-3]\nxarray2 = reproj_reflectivity_list[-2]\nxarray3 = reproj_reflectivity_list[-1]\n\ninitialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running ROVER and Semi-Lagrangian Advection\n\n1. Concatenate required reflectivity xarrays along time dimension.\n2. Run ROVER on all members, with the concatenated xarray as the input.\n3. Store motion field xarrays as a list.\n4. Perform Semi-Lagrangian Advection on all members using the motion fields\n from rover.\n5. Store forecasted reflectivities as list.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Combining two reflectivity DataArrays\n# the order of the coordinate keys is now ['y', 'x', 'time']\n# as opposed to ['time', 'x', 'y']\nreflec_concat_6min = xarray.concat([xarray2, xarray3], dim='time')\nreflec_concat_12min = xarray.concat([xarray1, xarray3], dim='time')\n\n# Running rover on 4 members\n# Mem-1\nu1, v1 = rover(\n reflec_concat_6min,\n start_level=1, max_level=7, sigma=2.5\n)\n# Mem-2\nu2, v2 = rover(\n reflec_concat_12min,\n start_level=2, max_level=7, sigma=2.5\n)\n# Rover-A\nu3, v3 = rover(\n reflec_concat_6min\n)\n# Mem-4\nu4, v4 = rover(\n reflec_concat_12min\n)\n\n# Storing motion fields for quiver plot\nmotion_list = [[u1, v1], [u2, v2], [u3, v3], [u4, v4]]\n\nrover_time = pd.Timestamp.now()\n\n# Running SLA on all members\n\nz1 = sla(\n reflec_concat_6min,\n u1, v1, steps=30\n)\n\nz2 = sla(\n reflec_concat_12min,\n u2, v2, steps=15\n)\n\nz3 = sla(\n reflec_concat_6min,\n u3, v3, steps=30\n)\n\nz4 = sla(\n reflec_concat_12min,\n u4, v4, steps=15\n)\n\n# appending all reflectivities to list\n\nz_sla_list = [z1, z2, z3, z4]\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.\n3. Concatenate reflectivities of different members along\n a fourth dimension.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "z_cat_list = []\nfor reflectivity in z_sla_list:\n z_all = reproj_reflectivity_list + [reflectivity[1:, ...]]\n z_cat = xarray.concat(z_all, dim='time')\n z_cat_list.append(z_cat)\n\n# Concatenating different members\nz_ens_6min = xarray.concat(\n [z_cat_list[0],\n z_cat_list[2]],\n xarray.IndexVariable('member', ['Mem-1', 'Mem-3'])\n)\n\nz_ens_12min = xarray.concat(\n [z_cat_list[1],\n z_cat_list[3]],\n xarray.IndexVariable('member', ['Mem-2', 'Mem-4'])\n)\nconcat_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating radar reflectivity maps\n\n1. Define the colour scale and the format of the\n reflectivity plots.\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. Set area extent.\n7. Plot GSHHS coastlines.\n8. Plot data using xarray.plot().\n9. Plot quiver using axes method.\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 levels\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# Plotting maps for forecasts with 6 minute tracking intervals\n\n# Defining the crs and coastline type\ncrs = area_def_tgt.to_cartopy_crs()\nhires = cfeature.GSHHSFeature(\n scale='h',\n levels=[1],\n edgecolor='black',\n facecolor='none'\n)\n\n# Define area extent\nzoom = (712000, 962000, 695000, 945000)\n\n# Generating a timelist for every hour\ntimelist = [\n (basetime + pd.Timedelta(minutes=60*i-6)).to_datetime64() for i in range(4)\n ]\n\nfor ensemble_mem in [z_ens_6min, z_ens_12min]:\n for idx, member in enumerate(ensemble_mem.coords['member'].values):\n images = []\n # Defining quiver for each member\n qx = motion_list[idx][0].coords['easting'].values[::5]\n qy = motion_list[idx][0].coords['northing'].values[::5]\n qu = motion_list[idx][0].values[::5, ::5]\n qv = motion_list[idx][1].values[::5, ::5]\n\n for time in timelist:\n # Figure\n fig = plt.figure(figsize=(28, 21))\n # Axes\n ax = plt.axes(projection=crs)\n # Zoom\n ax.set_extent(zoom, crs=crs)\n # Coastlines\n ax.add_feature(hires)\n # Plot data\n quadmesh = ensemble_mem.sel(\n member=member, time=time\n ).plot.pcolormesh(cmap=cmap, norm=norm)\n\n # Customize labels\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n ensemble_mem.attrs['long_name']+'[' +\n ensemble_mem.attrs['units']+']',\n fontsize=28\n )\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 # Plot quiver\n ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\n # Title\n plt.title(\n \"Quantitative Precipitation Forecast with Hourly Accumulated \"\n \"Rainfall\\n\"\n f\"Basetime: {basetime}H \"\n f\"Valid time: {pd.Timestamp(time)}H\\n\"\n f\"Member: {member}\",\n fontsize=32\n )\n savepath = THIS_DIR + \\\n (\"/../tests/outputs/rover-output-map-\"\n f\"{member}-{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png\")\n plt.savefig(savepath)\n images.append(imageio.imread(savepath))\n imageio.mimsave(\n THIS_DIR +\n f\"/../tests/outputs/rover-output-map-{member}.gif\",\n images, duration=1\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| Step 1: Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr().\n| Step 2: Replacing time coordinates with end times instead of start times.\n| Step 3: Interpolate rainrates temporally to ensure that rainrates are\n| equally spaced along the time axis.\n| Step 4: Convert rainrates to rainfalls in 6 and 12 minutes with rr2rf().\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 6 minute ensembles\nrainrates_ens_6min = dbz2rr(z_ens_6min, a=58.53, b=1.56)\nrainrates_ens_6min_endtime = rainrates_ens_6min.copy()\nrainrates_ens_6min_endtime.coords['time'] = \\\n [\n pd.Timestamp(t) + pd.Timedelta(minutes=6)\n for t in rainrates_ens_6min_endtime.coords['time'].values\n ]\nrainfalls_ens_6min = rr2rf(rainrates_ens_6min_endtime)\n\n# 12 minute ensembles\nrainrates_ens_12min = dbz2rr(z_ens_12min, a=58.53, b=1.56)\nrainrates_ens_12min_endtime = rainrates_ens_12min.copy()\nrainrates_ens_12min_endtime.coords['time'] = \\\n [\n pd.Timestamp(t) + pd.Timedelta(minutes=12)\n for t in rainrates_ens_12min_endtime.coords['time'].values\n ]\nrainrates_ens_12min_interp = temporal_interp(\n rainrates_ens_12min,\n pd.Timestamp('201902190700'),\n pd.Timestamp('201902191100')\n)\nrainfalls_ens_12min = rr2rf(rainrates_ens_12min_interp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Compute hourly accumulated rainfall every 30 minutes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "acc_rf_6min = multiple_acc(\n rainfalls_ens_6min, basetime, basetime+pd.Timedelta(hours=3)\n)\n\nacc_rf_12min = multiple_acc(\n rainfalls_ens_12min, basetime, basetime+pd.Timedelta(hours=3)\n)\n\nacc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculating Probability Exceeding Threshold\n\nIn this example, thresholds are 0.5mm, 5mm,\n10mm, 30mm, 50mm and 70 mm.\n\nProbabilities are 0%, 25%, 50%, 75% and 100%, as there are\nfour members.\n\nSteps:\n\n1. Define threshold.\n2. Concatenate rainfalls from 6 minute and 12 minute ensembles\n along the member dimension.\n3. Use a loop to calculate the probability of exceeding threshold\n for each gridcell and store results in list.\n4. Concatenate the contents of the list. Result is an xarray with\n dimensions (threshold, time, y, x).\n5. Plot the results using xarray.plot(). In this example,\n only the 0.5mm, 10mm and 50mm rainfall every hour will be plotted.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define threshold\nthreshold = [0.5, 5., 10., 30., 50., 70.]\n\n# Concatenating rainfalls_ens_6min and rainfalls_ens_12min\nacc_rf = xarray.concat(\n [acc_rf_6min, acc_rf_12min], dim='member'\n )\n\n# Define list to store probabilities of exceeding rainfall thresholds\n# List corresponds to threshold\nprob_list = []\n\n# Calculate probability\nfor th in threshold:\n bool_forecast = acc_rf >= th\n count = bool_forecast.sum(dim='member')\n prob = count/len(bool_forecast.coords['member'].values) * 100\n prob_list.append(prob)\n\n# Generate coordinate xarray for threshold\nth_xarray = xarray.IndexVariable(\n 'threshold', threshold\n)\n\n# concatenate\nprob_rainfall = xarray.concat(\n prob_list,\n dim=th_xarray\n)\n\nprob_rainfall.attrs['long_name'] = \"Probability of Exceeding Threshold\"\nprob_rainfall.attrs['units'] = \"%\"\n\n# Plot the results\ncmap = LinearSegmentedColormap.from_list(\n 'custom blue', ['#FFFFFF', '#000099']\n)\n\nfor th in prob_rainfall.coords['threshold'].values[::2]:\n images = []\n for time in prob_rainfall.coords['time'].values[2::2]:\n # Figure\n plt.figure(figsize=(28, 21))\n # Axes\n ax = plt.axes(projection=crs)\n # Zoom\n ax.set_extent(zoom, crs=crs)\n # Coastlines\n ax.add_feature(hires)\n\n # Plot data\n quadmesh = prob_rainfall.sel(threshold=th, time=time).plot(cmap=cmap)\n\n # Customising labels\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n prob_rainfall.attrs['long_name']+'[' +\n prob_rainfall.attrs['units']+']',\n fontsize=28\n )\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 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 \"Probability of Exceeding Threshold\\n\"\n f\"Basetime: {str(basetime)}H \"\n f\"Valid time: {pd.Timestamp(time)}H\\n\"\n f\"Start time: t{t_minus:+} minutes\"\n f\" End time: t{t_plus:+} minutes\\n\"\n f\"Threshold = {th}mm\",\n fontsize=32\n )\n savepath = THIS_DIR + \\\n (f\"/../tests/outputs/p_{pd.Timestamp(time).strftime('%Y%m%d%H%M')}\"\n f\"_threshold_{th}mm.png\")\n plt.savefig(savepath)\n images.append(imageio.imread(savepath))\n imageio.mimsave(\n THIS_DIR +\n f\"/../tests/outputs/p_{th}.gif\",\n images, duration=1\n )\n\nprob_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Rainfall percentiles\n\n1. Using xarray.DataArray.mean(),\n calculate the mean rainfall of all gridpoints.\n2. Using xarray.DataArray.min(),\n find the minimum rainfall of all gridpoints.\n3. Using xarray.DataArray.max(),\n find the maximum rainfall of all gridpoints.\n4. Using xarray.DataArray.quantile()\n find the 25th, 50th and 75th percentile rainfall\n of all gridpoints.\n5. Concatenate rainfall along percentile dimension.\n6. Plot results using xarray.plot(). In this example\n only the minimum, 50th and 75th percentile rainfall\n every hour will be plotted.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mean\nmean_rainfall = acc_rf.mean(dim='member', keep_attrs=True)\n# max\nmax_rainfall = acc_rf.max(dim='member', keep_attrs=True)\n# min\nmin_rainfall = acc_rf.min(dim='member', keep_attrs=True)\n# quartiles\nq_rainfall = acc_rf.quantile(\n [.25, .5, .75], dim='member',\n interpolation='nearest', keep_attrs=True)\n\n# generate index\npercentile = xarray.IndexVariable(\n 'percentile',\n ['Minimum',\n '25th Percentile',\n '50th Percentile',\n 'Mean',\n '75th Percentile',\n 'Maximum']\n )\n\n# concatenating\n\np_rainfall = xarray.concat(\n [min_rainfall,\n q_rainfall.sel(quantile=.25).squeeze().drop('quantile'),\n q_rainfall.sel(quantile=.50).squeeze().drop('quantile'),\n mean_rainfall,\n q_rainfall.sel(quantile=.75).squeeze().drop('quantile'),\n max_rainfall],\n dim=percentile\n)\n\np_rainfall.attrs['long_name'] = \"Hourly Accumulated Rainfall\"\np_rainfall.attrs['units'] = \"mm\"\n\n# Defining levels\n\n# 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# Plotting\nfor pos in p_rainfall.coords['percentile'].values[::2]:\n images = []\n for time in p_rainfall.coords['time'].values[2::2]:\n fig = plt.figure(figsize=(28, 21))\n ax = plt.axes(projection=crs)\n ax.add_feature(hires)\n ax.set_extent(zoom, crs=crs)\n\n quadmesh = p_rainfall.sel(percentile=pos, time=time).plot(\n cmap=cmap, norm=norm\n )\n cbar = quadmesh.colorbar\n cbar.ax.set_ylabel(\n p_rainfall.attrs['long_name']+'[' +\n p_rainfall.attrs['units']+']',\n fontsize=28\n )\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 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 f\"{pos} Rainfall Intensity\\n\"\n f\"Basetime: {str(basetime)}H \"\n f\"Valid time: {pd.Timestamp(time)}H\\n\"\n f\"Start time: t{t_minus:+} minutes \"\n f\"End time: t{t_plus:+} minutes\",\n fontsize=32\n )\n position = pos.split(\" \")[0]\n savepath = THIS_DIR + \\\n (\"/../tests/outputs/rainfall_\"\n f\"{pd.Timestamp(time).strftime('%Y%m%d%H%M')}\"\n f\"_{position}.png\")\n plt.savefig(savepath)\n images.append(imageio.imread(savepath))\n imageio.mimsave(\n THIS_DIR + f\"/../tests/outputs/rainfall_{position}.gif\",\n images, duration=1\n )\n\nptime = 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 timesteps (in this example, 30 minute intervals).\n3. Store rainfall values over time in an xarray.DataArray.\n4. Plot the time series of rainfall with boxplots at desired station.\n In this case, the 15th percentile member is plotted.\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 xarray.DataArray.\nstation_rf_list = []\nstation_name = []\nfor index, row in df.iterrows():\n station_rf_list.append(p_rainfall.sel(\n northing=row[1], easting=row[2],\n method='nearest'\n ).drop('northing').drop('easting'))\n station_name.append(row[0])\n\nstation_name_index = xarray.IndexVariable(\n 'ID', station_name\n )\nstation_rf = xarray.concat(\n station_rf_list,\n dim=station_name_index\n)\n\n# Extracting the 15th ranked station\nxr_15_percentile = station_rf.quantile(\n .15, dim='ID', interpolation='nearest').drop('quantile')\n\n# Plotting\n_, tax = plt.subplots(figsize=(20, 15))\nplt.plot(\n np.arange(1, len(xr_15_percentile.coords['time'].values) + 1),\n xr_15_percentile.loc['Mean'].values,\n 'ko-'\n) # plot line\n\n# Storing percentiles as dictionary to call ax.bxp for boxplot\nstats_list = []\nfor i in range(len(xr_15_percentile.coords['time'].values)):\n stats = {\n 'med': xr_15_percentile.loc['50th Percentile'].values[i],\n 'q1': xr_15_percentile.loc['25th Percentile'].values[i],\n 'q3': xr_15_percentile.loc['75th Percentile'].values[i],\n 'whislo': xr_15_percentile.loc['Maximum'].values[i],\n 'whishi': xr_15_percentile.loc['Minimum'].values[i]\n }\n stats_list.append(stats)\n\n# Plot boxplot\ntax.bxp(\n stats_list, showfliers=False\n)\n\n# Labels\nxcoords = xr_15_percentile.coords['time'].values\nxticklabels = [pd.to_datetime(str(t)).strftime(\"%-H:%M\") for t in xcoords]\ntax.set_xticklabels(xticklabels)\ntax.xaxis.set_tick_params(labelsize=20)\ntax.yaxis.set_tick_params(labelsize=20)\nplt.title('Time Series of Hourly Accumulated Rainfall', fontsize=25)\nplt.ylabel(\"Hourly Accumulated Rainfall [mm]\", fontsize=22)\nplt.xlabel(\"Time\", fontsize=18)\nplt.savefig(THIS_DIR+\"/../tests/outputs/pqpf_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(\n \"Calculate and plot probability exceeding threshold: \"\n f\"{prob_time}\"\n )\nprint(\n f\"Plotting rainfall maps: {ptime}\"\n)\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(\n \"Time to calculate and plot probability exceeding threshold: \"\n f\"{prob_time-acc_time}\"\n )\nprint(f\"Time to plot rainfall maps: {ptime-prob_time}\")\nprint(\n f\"Time to extract station data and plot time series: \"\n f\"{extract_time-ptime}\"\n )" ] } ], "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 }