{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# QPE (Manila)\n\nThis example demonstrates how to perform QPE,\nusing raingauge data from Manila and radar data from Subic.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definitions\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport time\nimport numpy as np\nimport pandas as pd\nimport copy\nimport xarray\nimport scipy\nimport pyproj\nimport matplotlib\nfrom matplotlib.colors import BoundaryNorm\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nfrom pyresample import utils\nfrom cartopy.io import shapereader\nimport cartopy.feature as cfeature\nimport cartopy.crs as ccrs\nimport cartopy.io.shapereader as shpreader\n\nfrom swirlspy.obs import Rain\nfrom swirlspy.rad.uf_ph import read_uf_ph\nfrom swirlspy.core.resample import grid_resample\nfrom swirlspy.qpe.rfmap import rg_interpolate, comp_qpe, show_raingauge\nfrom swirlspy.qpe.utils import timestamps_ending, locate_file, dbz2rr, rr2rf, \\\n temporal_interp, acc\n\nplt.switch_backend('agg')\nTHIS_DIR = os.getcwd()\nos.chdir(THIS_DIR)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialising\n\nThis section demonstrates extracting raingauge and radar data.\n\nStep 1: Defining an end-time for accumulating rainfall.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "acctime = pd.Timestamp('20180811112000').floor('min')\nacctime_str = acctime.strftime('%Y%m%d%H%M')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Setting up directories for raingauge and radar files.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rad_dir = THIS_DIR + '/../tests/samples/uf_ph/sub/'\nrg_dir = THIS_DIR + '/../tests/samples/rfmap/'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Generating timestamps and pattern for both radar and\nraingauge files.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Timestamps of raingauges\nrg_timestrings = timestamps_ending(\n acctime + pd.Timedelta(minutes=10),\n duration=pd.Timedelta(hours=1),\n interval=pd.Timedelta(minutes=10)\n)\n\n# Raingauge pattern\nrg_pattern = ['rf60m_20'+ts for ts in rg_timestrings]\n\n# Finding time nearest radar file\n# to accumulation end time\nminute = acctime.minute\nnearest_6min = acctime.minute // 10 * 10\n\nnearest_rad_timestamp = pd.Timestamp(\n acctime_str[:-2]+f'{nearest_6min:02}'\n)\n\nrad_timestrings = timestamps_ending(\n nearest_rad_timestamp,\n duration=pd.Timedelta(hours=1),\n interval=pd.Timedelta(minutes=10)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: Extracting raingauge and radar files from\ntheir respective directories.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "located_rg_files = []\nfor pat in rg_pattern:\n located_rg_files.append(locate_file(rg_dir, pat))\n\n\nlocated_radar_files = []\nfor ts in rad_timestrings:\n located_radar_files.append(locate_file(rad_dir, ts))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 5: Read data from raingauge files into a Rain object.\nCoordinates are geodetic, following that in the files.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rg_object_geodetic = Rain(\n located_rg_files,\n 'WGS84',\n duration=pd.Timedelta(minutes=5),\n NAN=[3276.7, 32767]\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 6: 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 = \"epsg3123_240km\"\ndescription = (\"A 240 m resolution rectangular grid \"\n \"centred at Subic RADAR and extending to 240 km \"\n \"in each direction\")\nproj_id = 'epsg3123'\nprojection = ('+proj=tmerc +lat_0=0 '\n '+lon_0=121 +k=0.99995 +x_0=500000 '\n '+y_0=0 +ellps=clrk66 +towgs84=-127.62,-67.24,'\n '-47.04,-3.068,4.903,1.578,-1.06 +units=m '\n '+no_defs')\nx_size = 1000\ny_size = 1000\narea_extent = (191376.04113, 1399386.68659, 671376.04113, 1879386.68659)\narea_def = utils.get_area_def(\n area_id, description, proj_id, projection, x_size, y_size, area_extent\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 7: Convert coordinates of raingauge object to desired projection.\nIn this example, the desired projection is PRS92.\nThis can be achieved by the .reproject() method of the Rain object.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "inProj = pyproj.Proj(init=\"epsg:4326\")\noutProj = pyproj.Proj(area_def.proj_str)\n\nrg_object = rg_object_geodetic.reproject(inProj, outProj, \"PRS92\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 8: Read radar files into xarray.DataArrays using read_uf_ph().\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reflec_list = []\nfor file in located_radar_files:\n reflec = read_uf_ph(\n file, area_def=area_def,\n coord_label=['easting', 'northing'],\n indicator='deg_ppi', elevation=0.5\n )\n\n reflec_list.append(reflec)\n\nreflectivity = xarray.concat(reflec_list, dim='time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accumulating and interpolating rainfall\n\nInterpolate rainfall recorded by raingauges into the user-defined grid\nand accumulate radar rainfall over an hour\nafter making the necessary\nadjustments.\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 1: Interpolate Rain object to user-defined grid.\nIn this example, a multiquadric Radial Basis Function\nis used.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "interpolated_rg = rg_interpolate(\n rg_object, area_def, 'rbf',\n coord_label=['easting', 'northing']\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Convert to radar reflectivity to rainrates,\nconvert rainrates to times of raingauges,\nand accumulate rainfalls every 10 minutes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rainrates = dbz2rr(reflectivity, a=300, b=1.4)\n# Convert time coordinates of rainrates from start time\n# to end time\nrainrates_time_endtime = rainrates.copy()\nrainrates_time_endtime.coords['time'] = \\\n [\n pd.Timestamp(t) + pd.Timedelta(minutes=10)\n for t in rainrates.coords['time'].values\n ]\n\nrainfalls = rr2rf(rainrates_time_endtime, scan_duration=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Accumulate rainfall over an hour.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "acc_rf = acc(\n rainfalls,\n rg_object.end_time,\n acc_period=pd.Timedelta(minutes=60)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compositing rainfall\n\nPerform compositing on radar and raingauge derived rainfall\nto obtain a composite QPE.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "comprf = comp_qpe(\n area_def,\n rg_object=rg_object,\n rg_interp=interpolated_rg,\n rad_rf=acc_rf\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting\n\nPlot composited radar and raingauge rainfall.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 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 function for neatness.\n\ndef plot_map(\n da, rg_object, acctime, area_def,\n based='raingauge and radar',\n savepath='', area_extent=None\n):\n\n \"\"\"\n A custom function for plotting a map.\n\n Parameters\n --------------\n da: xarray.DataArray\n Contains data to be plotted.\n rg_object: Rain\n Contains raingauge data.\n acctime: pd.Timestamp\n Contains the endtime of the accumulation\n period.\n area_def: pyresample.geometry.AreaDefinition\n AreaDefinition of the grid.\n based: str\n Type of data plotted in the map.\n savepath: str\n Path to save the image to.\n area_extent: tuple\n Area extent (x0, x1, y0, y1) to be plotted.\n Defaults to None.\n\n \"\"\"\n # Defining the colour scheme\n levels = [\n 0, 0.5, 2, 5, 10, 20,\n 30, 40, 50, 70, 100, 150,\n 200, 300, 400, 500, 600, 700\n ]\n\n cmap = 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\n norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n # Plotting axes\n plt.figure(figsize=(24, 21))\n crs = area_def.to_cartopy_crs()\n ax = plt.axes(projection=crs)\n\n if area_extent is not None:\n ax.set_extent(area_extent, crs=crs)\n\n # Plot data\n quadmesh = da.plot(\n cmap=cmap,\n norm=norm,\n extend='max',\n cbar_kwargs={'ticks': levels, 'format': '%.3g'}\n )\n\n # Adjusting size of colorbar\n cb = quadmesh.colorbar\n cb.ax.set_ylabel(\n da.attrs['long_name']+'['+da.attrs['units']+']',\n fontsize=28\n )\n\n cb.ax.tick_params(labelsize=24)\n\n # Setting labels\n ax.xaxis.set_visible(True)\n ax.yaxis.set_visible(True)\n\n for tick in ax.xaxis.get_major_ticks():\n tick.label.set_fontsize(24)\n for tick in ax.yaxis.get_major_ticks():\n tick.label.set_fontsize(24)\n\n ax.xaxis.label.set_size(28)\n ax.yaxis.label.set_size(28)\n\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\n # Show raingauge\n show_raingauge(\n rg_object, ax, show_value=True, color='red', markersize=20,\n fontsize=20\n )\n # Show title\n plt.title(\n (f\"Last Hour Rainfall at {acctime.strftime('%H:%MH %d-%b-%Y')}\\n\"\n f\"(based on {based} data)\"),\n fontsize=32\n )\n\n plt.savefig(savepath)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Plotting maps\n\nr = 64000\nproj_site = reflectivity.proj_site\narea_extent = (\n proj_site[0]-r,\n proj_site[0]+r,\n proj_site[1]-r,\n proj_site[1]+r\n)\n\n# Raingauge only\nplot_map(\n interpolated_rg, rg_object,\n acctime, area_def,\n based='raingauge',\n savepath=THIS_DIR+f'/../tests/outputs/raingauge_{acctime_str}.png',\n area_extent=area_extent\n)\n\n# Radar only\nplot_map(\n acc_rf, rg_object,\n acctime, area_def,\n based='radar',\n savepath=THIS_DIR+f'/../tests/outputs/radar_{acctime_str}.png',\n area_extent=area_extent\n)\n\n# Composite raingauge and radar\nplot_map(\n comprf, rg_object,\n acctime, area_def,\n savepath=THIS_DIR+f'/../tests/outputs/comp_{acctime_str}.png',\n area_extent=area_extent\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 }