{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nQPF (Thailand)\n===================================================\nThis example demonstrates how to perform\noperational deterministic QPF up to three hours\nfrom hourly composite rainfall data, using data from Thailand.\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\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 projection conversion\nfrom utm.conversion import from_latlon\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\n# Python package for area definition\nfrom pyresample import get_area_def\n\n# Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection\nfrom swirlspy.qpf import rover, sla\n# swirlspy Thailand netcdf file parser function\nfrom swirlspy.rad import read_netcdf_th_refl\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_depth, acc_rainfall_depth\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 a basetime.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Supply basetime\nbaseT = \"202010161000\"\nbasetime = pd.Timestamp(baseT)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 2: Using basetime, generate timestamps of desired radar files.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Obtain radar files\nlocated_files = []\ninterval = pd.Timedelta(15, 'm')\nduration = pd.Timedelta(60, 'm')\n\nwhile duration >= pd.Timedelta(0):\n located_files.append(\n os.path.abspath(os.path.join(\n DATA_DIR,\n 'netcdf_tmd',\n (basetime - duration).strftime('RADAR_Thailand2_%Y%m%dT%H%M%S000Z.nc')\n )\n ))\n duration -= interval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3: Read data from radar files into xarray.DataArray\nusing read_netcdf_th_refl().\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# defining the a and b value of the Z-R relationship\nZRa = 300\nZRb = 1.4\n\n\nreflect_list = [] # stores reflectivity from read_netcdf_th-refl()\nfor filename in located_files:\n reflec = read_netcdf_th_refl(\n filename, ZRa, ZRb\n )\n reflect_list.append(reflec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4 : Define the target grid\n\nDefining target grid\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "area_id = \"Thai1975\"\ndescription = (\"Covering whole territory of Thailand\")\nproj_id = 'Thai'\nprojection = \"+proj=utm +zone=47 +a=6377276.345 +b=6356075.41314024 +towgs84=210,814,289,0,0,0,0 +units=m +no_defs \"\nx_size = len(reflect_list[0][0][0])\ny_size = len(reflect_list[0][0][:])\n\nlon = reflect_list[0][0][0].coords['lon']\nlat = reflect_list[0][0].coords['lat']\n\n# Using utm to calculate the coordinate system\nll = from_latlon(float(lat[-1]), float(lon[0]), 47)\nur = from_latlon(float(lat[0]), float(lon[-1]), 47)\n\narea_extent = (ll[0], ll[1], ur[0], ur[1])\n\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_netcdf_th_refl()\n(source) projection to Thai (target) projection.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Extracting the AreaDefinition of the source projection\n\narea_def_src = reflect_list[0].attrs['area_def']\n\nreproj_reflect_list = []\nfor reflect in reflect_list:\n reproj_reflect = grid_resample(\n reflect, area_def_src, area_def_tgt,\n coord_label=['easting', 'northing']\n )\n reproj_reflect_list.append(reproj_reflect)" ] }, { "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-------------------------------------------\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']\n\nreflect_concat = xr.concat(reproj_reflect_list, dim='time')\nstandardize_attr(reflect_concat, frame_type=FrameType.dBZ, zero_value=9999.)\n\n\n# Rover\nmotion = rover(reflect_concat)\n\nrover_time = pd.Timestamp.now()\n\n\n# Semi Lagrangian Advection\nreflectivity = sla(reflect_concat, motion, 12)\n\nsla_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Concatenating observed and forecasted reflectivities\n-------------------------------------------------------\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([reflect_concat[:-1, ...], reflectivity], dim='time')\nreflectivity.attrs['long_name'] = '2 km radar composite reflectivity'\nstandardize_attr(reflectivity)\n\nconcat_time = pd.Timestamp.now()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating radar reflectivity maps\n-----------------------------------\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\n\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n# Defining the crs\ncrs = area_def_tgt.to_cartopy_crs()\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\n# Generating a timelist for every hour\ntimelist = [\n (basetime + pd.Timedelta(minutes=60*i-15)) for i in range(4)\n]\n\n\n# Obtaining the slice of the xarray to be plotted\nda_plot = reflectivity.sel(time=timelist)\n\n\n# Defining motion quivers\nqx = motion.coords['easting'].values[::50]\nqy = motion.coords['northing'].values[::50]\nqu = motion.values[0, ::50, ::50]\nqv = motion.values[1, ::50, ::50]\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)\n\n\nfor 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', linestyle=':', linewidth=0.5)\n # precipitation\n ax.quiver(qx, qy, qu, qv, pivot='mid', scale=20, scale_units='inches')\n ax.gridlines()\n ax.set_title(\n \"Precipitation\\n\"\n f\"Based @ {basetime.strftime('%H:%MZ')}\",\n loc='left',\n fontsize=7\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:%MZ')} \",\n loc='right',\n fontsize=7\n )\n\n\nplt.savefig(\n os.path.join(OUTPUT_DIR, f\"rover-output-map-mn{baseT}.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------------------------------------------------\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=ZRa, b=ZRb)\n\n# Converting the coordinates of xarray from start to endtime\nrainfalls.coords['time'] = [\n pd.Timestamp(t) + pd.Timedelta(15, 'm')\n for t in rainfalls.coords['time'].values\n]\n\n\nintervalstep = 30\nresult_step_size = pd.Timedelta(minutes=intervalstep)\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 result_step_size=result_step_size\n)\nacc_rf.attrs['long_name'] = 'Rainfall accumulated over the past 60 minutes'\n\n\nacc_time = pd.Timestamp.now()\n\n\n# Defining times for plotting\ntimelist = [basetime + pd.Timedelta(i, 'h') for i in range(4)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting rainfall maps\n---------------------------------------\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\nlevels2 = [\n 0., 0.5, 2., 5., 10., 20.,\n 30., 40., 50., 70., 100., 150.,\n 200., 300., 400., 500., 600., 700.\n]\n\ncmap2 = ListedColormap([\n '#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646',\n '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063',\n '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd'\n])\n\nnorm2 = BoundaryNorm(levels2, ncolors=cmap2.N, clip=True)\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 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': levels2,\n 'format': '%.3g'\n },\n cmap=cmap2,\n norm=norm2\n)\nfor 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', linestyle=':', linewidth=0.5)\n ax.gridlines()\n ax.set_title(\n \"Past Hour Rainfall\\n\"\n f\"Based @ {basetime.strftime('%H:%MZ')}\",\n loc='left',\n fontsize=7\n )\n ax.set_title(\n ''\n )\n ax.set_title(\n f\"{timelist[idx].strftime('%Y-%m-%d')} \\n\"\n f\"Valid @ {timelist[idx].strftime('%H:%MZ')} \",\n loc='right',\n fontsize=7\n )\n\n\nplt.savefig(\n os.path.join(OUTPUT_DIR, f\"rainfall_{baseT}.png\"),\n dpi=300\n)\n\nrf_image_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 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}\")\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\"Total Running time: {rf_image_time - start_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 }