{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nQPE (Manila)\n====================================================\n\nThis example demonstrates how to perform QPE,\nusing raingauge data from Manila and radar data from Subic.\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 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\nimport pyproj\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\n\n# swirlspy raingauge data object\nfrom swirlspy.obs import Rain\n# swirlspy Philippine UF file parser function\nfrom swirlspy.rad.uf_ph import read_uf_ph\n# swirlspy raingauge data interpolate and blending\nfrom swirlspy.qpe.rfmap import rg_interpolate, comp_qpe, show_raingauge\n# swirlspy test data source locat utils function\nfrom swirlspy.qpe.utils import timestamps_ending, locate_file\n# swirlspy standardize data function\nfrom swirlspy.utils import FrameType, standardize_attr, conversion\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')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initialising\n-----------------------------------------------------------------\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 = os.path.join(DATA_DIR, 'uf_ph/sub')\nrg_dir = os.path.join(DATA_DIR, '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 = 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 = xr.concat(reflec_list, dim='time')\nstandardize_attr(reflectivity, frame_type=FrameType.dBZ)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Accumulating and interpolating rainfall\n-----------------------------------------------------------------\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 = conversion.to_rainfall_rate(reflectivity, False, a=300, b=1.4)\n\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\nstandardize_attr(rainrates_time_endtime, frame_type=FrameType.mmph)\nrainfalls = conversion.to_rainfall_depth(rainrates_time_endtime)" ] }, { "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 = conversion.acc_rainfall_depth(\n rainfalls,\n rg_object.end_time,\n rg_object.end_time,\n pd.Timedelta(minutes=60)\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compositing rainfall\n-----------------------------------------------------------------\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---------------------------------------------------------------\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(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\n# Plotting function for neatness.\n\n\ndef plot_map(\n da, rg_object, acctime, area_def,\n based='raingauge and radar',\n savepath='', area_extent=None\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 = 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\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=os.path.join(OUTPUT_DIR, f'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=os.path.join(OUTPUT_DIR, f'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=os.path.join(OUTPUT_DIR, f'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 }