{
"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 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 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_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\nhires = cfeature.GSHHSFeature(\n levels=[1],\n scale='h',\n edgecolor='k'\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 ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\n ax.add_feature(hires) # coastlines\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 ax.add_feature(hires) # using GSHHS coastlines defined previously\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
}