{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nPQPF (Hong Kong)\n========================================================\nThis example demonstrates how to perform\noperational PQPF up to three hours from\nraingauge and radar data, using data from Hong Kong.\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 ordered dictionary\nfrom collections import OrderedDict\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 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 colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap, LinearSegmentedColormap\n\n# Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection\nfrom swirlspy.qpf import rover, sla\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 standardize_attr, FrameType\n# swirlspy data convertion function\nfrom swirlspy.utils.conversion import to_rainfall_rate, to_rainfall_depth, acc_rainfall_depth, temporal_interpolate\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")\nplt.switch_backend('agg')\n\nstart_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initialising\n-----------------------------------------------------------\n\nThis section demonstrates extracting\nradar reflectivity data.\n\nStep 1: Define the basetime.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define 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 = os.path.join(DATA_DIR, '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" ] }, { "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 = 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 (source) projection to HK 1980 (target) projection.\n\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)\n\ninitialising_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: Concatenate the reflectivity xarrays in the list along the time\ndimension.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Concatenating all the (observed) reflectivity\n# xarray.DataArrays in reproj_reflectivity_list\nreflectivityObs = xr.concat(reproj_reflectivity_list, dim='time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Running ROVER and Semi-Lagrangian Advection\n---------------------------------------------------------\n\nThis section demonstrates how to run ROVER and\nSemi-Lagrangian Advection for multiple members.\n\nFirst, write a function to run these steps:\n\n#. Selecting the two relevant xarrays to be used as the ROVER input.\n#. Generate motion field using ROVER, with the selected xarray as the input.\n#. Perform Semi-Lagrangian Advection using the motion fields from ROVER.\n\nWriting the function makes it easier to implement multiprocessing if\nrequired.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ensemble_rover_sla(\n reflectivityObs,\n basetime,\n interval,\n duration,\n member,\n start_level,\n max_level,\n rho, alpha, sigma,\n track_interval\n):\n \"\"\"\n A function to run ROVER and SLA.\n\n Parameters\n -------------\n reflectivityObs: xarray.DataArray\n The xarray containing observed reflectivity.\n basetime: pandas.Timestamp\n The basetime of the forecast.\n Equivalent to the end time of the latest radar scan.\n interval: pandas.Timedelta\n The interval between the radar scans.\n duration: pandas.Timedelta\n The duration of the forecast.\n member: str\n The name of the member.\n start_level, max_level, rho, alpha, sigma, track_interval: floats\n Parameters of the members.\n\n Returns\n ------------\n reflectivityFcst: xarray.DataArray\n Contains foreacsted reflectivity.\n \"\"\"\n\n # Generate timestamps of relevant radar scans\n latestRadarScanTime = basetime - interval\n\n roverTimestrings = timestamps_ending(\n latestRadarScanTime, duration=track_interval,\n interval=track_interval,\n format='%Y%m%d%H%M',\n exclude_end=False\n )\n\n roverTimestamps = [pd.Timestamp(t) for t in roverTimestrings]\n\n qpf_xarray = reflectivityObs.sel(time=roverTimestamps)\n standardize_attr(qpf_xarray, frame_type=FrameType.dBZ, zero_value=9999.0)\n\n # Generate motion field using ROVER\n motion = rover(\n qpf_xarray,\n start_level=start_level,\n max_level=max_level,\n rho=rho,\n sigma=sigma,\n alpha=alpha\n )\n\n # Perform SLA\n steps = int(duration / interval)\n reflectivityFcst = sla(qpf_xarray, motion, steps)\n\n reflectivityFcst = reflectivityFcst.expand_dims(\n dim=OrderedDict({'member': [member]})\n ).copy()\n\n return reflectivityFcst" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, run the different members by calling the function above.\nIn this example. only four members are run.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Define the member parameters\n# Ensure that values at the same position across\n# the tuples correspond to the same member\nstartLevelTuple = (1, 2, 1, 1)\nmaxLevelTuple = (7, 7, 7, 7)\nsigmaTuple = (2.5, 2.5, 1.5, 1.5)\nrhoTuple = (9., 9., 9., 9.)\nalphaTuple = (2000., 2000., 2000., 2000.)\nmemberTuple = ('Mem-1', 'Mem-2', 'Mem-3', 'Mem-4')\ntrackIntervalList = [pd.Timedelta(i, 'm') for i in [6, 12, 6, 12]]\n\nprocesses = 4\n\nreflectivityFcstList = []\nfor i in range(processes):\n args = (\n reflectivityObs,\n basetime,\n pd.Timedelta(6, 'm'),\n pd.Timedelta(3, 'h'),\n memberTuple[i],\n startLevelTuple[i],\n maxLevelTuple[i],\n rhoTuple[i],\n alphaTuple[i],\n sigmaTuple[i],\n trackIntervalList[i]\n )\n reflectivityFcst = ensemble_rover_sla(*args)\n reflectivityFcstList.append(reflectivityFcst)\n\nrover_sla_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Concatenating observed and forecasted reflectivities\n--------------------------------------------------------\n\nFrom the results from all members:\n\n#. Change the timestamps of observed reflectivity\n and forecasted reflectivity of all members\n from start time to end time.\n#. Add forecasted reflectivity to reproj_reflectivity_list.\n#. 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": [ "# List to store combined reflectivity arrays of each member\nreflectivityConcatList = []\n\n# Changing the time coordinates from start time to end time\nreflectivityObs.coords['time'] = [\n pd.Timestamp(t) + pd.Timedelta(6, 'm')\n for t in reflectivityObs.time.values\n]\n\nfor reflectivityFcst in reflectivityFcstList:\n # Identify the member\n member = reflectivityFcst.member.values[0]\n # Add member dimension to observed reflectivity\n reflectivityObsExpanded = reflectivityObs.expand_dims(\n dim=OrderedDict({'member': [member]})\n ).copy()\n\n # Checking for the track interval of the forecasted\n # array\n timeIntervalsSet = set(np.diff(reflectivityFcst.time.values))\n if pd.Timedelta(12, 'm').to_timedelta64() in timeIntervalsSet:\n reflectivityFcst.coords['time'] = [\n pd.Timestamp(t) + pd.Timedelta(12, 'm')\n for t in reflectivityFcst.time.values\n ]\n else:\n reflectivityFcst.coords['time'] = [\n pd.Timestamp(t) + pd.Timedelta(6, 'm')\n for t in reflectivityFcst.time.values\n ]\n\n # Combining\n reflectivityList = [\n reflectivityObsExpanded,\n reflectivityFcst.isel(time=slice(1, None))\n ]\n reflectivity = xr.concat(reflectivityList, dim='time')\n standardize_attr(reflectivity, frame_type=FrameType.dBZ,\n zero_value=reflectivityFcst.attrs['zero_value'])\n reflectivityConcatList.append(reflectivity)\n\nconcat_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+0min.\n\n#. Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr().\n#. For rainrates every 12 minutes, interpolate rainfall\n temporally to once every 6 minutes to ensure that rainrates are\n equally spaced along the time axis.\n#. Convert rainrates to rainfalls with rr2rf().\n#. Compute hourly accumulated rainfall every 30 minutes.\n#. Concatenate accumulated rainfall of different members along the\n member dimension.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "accList = []\nfor reflectivity in reflectivityConcatList:\n # Converting reflectivity to rainrates\n rainrates = to_rainfall_rate(\n reflectivity, logarithmic=False, a=58.53, b=1.56)\n\n # Interpolate rainfall to every 12 minutes if\n # the member has a track interval of 12 minutes\n if pd.Timedelta(12, 'm').to_timedelta64() in timeIntervalsSet:\n rainrates = temporal_interpolate(rainrates,\n basetime - pd.Timedelta(1, 'h'),\n basetime + pd.Timedelta(3, 'h'),\n pd.Timedelta(6, 'm'),\n 'linear')\n # Convert rainrates to rainfalls every 6 minutes\n rainfalls = to_rainfall_depth(rainrates)\n\n # Accumulate rainfalls every 30 minutes\n accr = acc_rainfall_depth(rainfalls,\n basetime,\n basetime + pd.Timedelta(3, 'h'),\n pd.Timedelta(30, 'm'),\n pd.Timedelta(1, 'h'))\n accList.append(accr)\n\n# Concatenating along the mmeber dimension\nacc_rf = xr.concat(accList, dim='member')\n\nacc_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculating Probability Exceeding Threshold\n----------------------------------------------\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\n#. Define threshold.\n#. Use a loop to calculate the probability of exceeding threshold\n for each gridcell and store results in list.\n#. Concatenate the contents of the list. Result is an xarray with\n dimensions (threshold, time, y, x).\n#. Plot the results using xarray.plot(). In this example,\n only the probability exceeding 0.5mm rainfall\n 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# 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 = xr.IndexVariable('threshold', threshold)\n\n# concatenate\nprob_rainfall = xr.concat(prob_list, dim=th_xarray)\n\nprob_rainfall.attrs['long_name'] = \"Probability of Exceeding Threshold\"\nprob_rainfall.attrs['units'] = \"%\"\n\n# Plot the results\n\n# Extracting the crs\ncrs = area_def_tgt.to_cartopy_crs()\n\n# Defining cmap\ncmap = LinearSegmentedColormap.from_list(\n 'custom blue', ['#FFFFFF', '#000099']\n)\n\n# Obtaining xarray slices to be plotted\nthreshold_list = [0.5]\ntimelist = [basetime + pd.Timedelta(hours=i) for i in range(4)]\nda_plot = prob_rainfall.sel(threshold=threshold_list, time=timelist)\n\n# Defining coastlines\nmap_shape_file = os.path.join(DATA_DIR, \"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\nfor th in threshold_list:\n # Plotting\n p = da_plot.sel(threshold=th).plot(\n col='time', col_wrap=2,\n subplot_kws={'projection': crs},\n cbar_kwargs={'ticks': [0, 25, 50, 75, 100],\n 'format': '%.3g'},\n cmap=cmap\n )\n\n for 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() # gridlines\n ax.set_title(\n f\"% Exceeding {th}mm 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 )\n plt.savefig(\n os.path.join(OUTPUT_DIR, f\"p_{th}.png\"),\n dpi=300\n )\n\nprob_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rainfall percentiles\n----------------------------\n\n#. Using xarray.DataArray.mean(),\n calculate the mean rainfall of all gridpoints.\n#. Using xarray.DataArray.min(),\n find the minimum rainfall of all gridpoints.\n#. Using xarray.DataArray.max(),\n find the maximum rainfall of all gridpoints.\n#. Using xarray.DataArray.quantile()\n find the 25th, 50th and 75th percentile rainfall\n of all gridpoints.\n#. Concatenate rainfall along percentile dimension.\n#. Plot results using xarray.plot(). In this example\n only the maximum 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 = xr.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 = xr.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 = 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# Obtaining xarray slices to be plotted\nplot_positions = ['Maximum']\ntimelist = [basetime + pd.Timedelta(hours=i) for i in range(4)]\nda_plot = p_rainfall.sel(percentile=plot_positions, time=timelist)\n\nfor pos in plot_positions:\n # Plotting\n p = da_plot.sel(percentile=pos).plot(\n col='time', col_wrap=2,\n subplot_kws={'projection': crs},\n cbar_kwargs={'ticks': levels,\n 'format': '%.3g'},\n cmap=cmap,\n norm=norm\n )\n\n for 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_title(\n f\"{pos} Radar-Based 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 )\n position = pos.split(\" \")[0]\n plt.savefig(\n os.path.join(OUTPUT_DIR, f\"rainfall_{position}.png\"),\n dpi=300\n )\nptime = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract the rainfall values at a specified location\n------------------------------------------------------------------\n\nIn this example, the rainfall values at the location is assumed\nto be the same as the nearest gridpoint.\n\n#. Read information regarding the radar stations into a pandas.DataFrame.\n#. Extract the rainfall values at the nearest gridpoint to location\n for given timesteps (in this example, 30 minute intervals).\n#. Store rainfall values over time in an xarray.DataArray.\n#. 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(DATA_DIR, \"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 = xr.IndexVariable(\n 'ID', station_name\n)\nstation_rf = xr.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(\n os.path.join(OUTPUT_DIR, \"pqpf_time_series.png\")\n)\n\nextract_time = 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\"Rover and SLA time: {rover_sla_time}\")\nprint(f\"Concatenating time: {concat_time}\")\nprint(f\"Accumulating rainfall time: {acc_time}\")\nprint(\"Calculate and plot probability exceeding threshold: \"\n f\"{prob_time}\")\nprint(f\"Plotting rainfall maps: {ptime}\")\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 and SLA: {rover_sla_time-initialising_time}\")\nprint(f\"Time to concatenate xarrays: {concat_time - rover_sla_time}\")\nprint(f\"Time to accumulate rainfall: {acc_time - concat_time}\")\nprint(\"Time to calculate and plot probability exceeding threshold: \"\n f\"{prob_time-acc_time}\")\nprint(f\"Time to plot rainfall maps: {ptime-prob_time}\")\nprint(f\"Time to extract station data and plot time series: \"\n f\"{extract_time-ptime}\")" ] } ], "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 }