{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nQPF (Laos)\n========================================================\nThis example demonstrates how to perform\noperational deterministic QPF up to three hours using \nreflectivity data.\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 timestamp\nimport pandas as pd\n# Python package for xarrays to read and handle netcdf data\nimport xarray as xr\n# Python package for numerical calculations\nimport numpy as np\n# Python com-swirls package to standardize attributes\nfrom swirlspy.utils import standardize_attr, FrameType\n# Python com-swirls package to calculate motion field (rover) and semi-lagrangian advection\nfrom swirlspy.qpf import rover, sla\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")\n\n# data files\ndata_paths = [\n    os.path.abspath(os.path.join(DATA_DIR, 'laos_h8/laos_20190731070000.nc')),\n    os.path.abspath(os.path.join(DATA_DIR, 'laos_h8/laos_20190731071000.nc')),\n]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialising\n---------------------------------------------------\n\nThis section demonstrates preprocessing\nreflectivity data.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read data from saved files\nreflec_datas = []\nfor file_path in data_paths:\n    d = xr.open_dataarray(file_path)\n    reflec_datas.append(d)\n\n# Concatenate list on time axis\ndata_frames = xr.concat(reflec_datas, dim='time')\n\n# Make sure data is ordered by time\ndata_frames = data_frames.sortby('time', ascending=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Nowcast (SWIRLS-Radar-Advection)\n---------------------------------------------------\nThe swirls radar advection was performed using the observed radar data\nFirstly, some attributes necessary for com-swirls input variable is added\nSecondly, rover function is invoked to compute the motion field\nThirdly, semi-lagrangian advection is performed to advect the radar data using the rover motion field\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Adding in some attributes that is step_size <10 mins in pandas.Timedelta>, zero_value <numpy.nan> and frame_type <FrameType.dBZ>\nstandardize_attr(data_frames)\n\n# Rover motion field computation\nmotion = rover(data_frames)\n\nrover_time = pd.Timestamp.now()\n\n# Semi-Lagrangian Advection\nswirls = sla(data_frames, motion, 17)  # Radar time goes from earliest to latest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting result\n---------------------------------------------------\nStep 1: Import plotting library and necessary library\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Python package for reading map shape file\nimport cartopy.io.shapereader as shpreader\n# Python package for land/sea features\nimport cartopy.feature as cfeature\n# Python package for projection\nimport cartopy.crs as ccrs\n# Python package for creating plots\nfrom matplotlib import pyplot as plt\n# Python package for output import grid\nfrom matplotlib.gridspec import GridSpec\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\n# Python package for scalar data to RGBA mapping\nfrom matplotlib.cm import ScalarMappable\n\nplt.switch_backend('agg')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 2: Defining the dBZ levels, colorbar parameters and projection\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# levels of colorbar (dBZ)\nlevels = [-32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44,\n          47, 50, 53, 56, 58, 60, 62]\n# hko colormap for dBZ at each levels\ncmap = ListedColormap([\n    '#FFFFFF00', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n    '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n    '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n])\n# boundary\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n# scalar data to RGBA mapping\nscalar_map = ScalarMappable(cmap=cmap, norm=norm)\nscalar_map.set_array([])\n# Defining plot parameters\nmap_shape_file = os.path.abspath(os.path.join(\n    DATA_DIR,\n    'shape/se_asia'\n))\n# coastline and province\nse_asia = cfeature.ShapelyFeature(\n    list(shpreader.Reader(map_shape_file).geometries()),\n    ccrs.PlateCarree()\n)\n# output area\nextents = [99.5, 108., 13.75, 22.75]\n\n# base_map plotting function\ndef plot_base(ax: plt.Axes):\n    ax.set_extent(extents, crs=ccrs.PlateCarree())\n\n    # fake the ocean color\n    ax.imshow(np.tile(np.array([[[178, 208, 254]]],\n                               dtype=np.uint8), [2, 2, 1]),\n              origin='upper',\n              transform=ccrs.PlateCarree(),\n              extent=[-180, 180, -180, 180],\n              zorder=-1)\n    # coastline, state, color\n    ax.add_feature(se_asia,\n                   facecolor=cfeature.COLORS['land'], edgecolor='none', zorder=0)\n    # overlay coastline, state without color\n    ax.add_feature(se_asia, facecolor='none',\n                   edgecolor='gray', linewidth=0.5)\n\n    ax.set_title('')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 3: Filtering\nvalues <= 15dbZ are not plotted\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "swirls = swirls.where(swirls > 15, np.nan)\n\n# Defining motion quivers\nqx = motion.coords['x'].values[::5]\nqy = motion.coords['y'].values[::5]\nqu = motion.values[0, ::5, ::5]\nqv = motion.values[1, ::5, ::5]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 4: Plotting the swirls-radar-advection, nwp-bias-corrected, blended 3 hours ahead\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig: plt.Figure = plt.figure(\n    figsize=(3 * 3.5 + 1, 3),\n    frameon=False\n)\ngs = GridSpec(\n    1, 3, figure=fig,\n    wspace=0, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nbasetime = pd.Timestamp(swirls.time.values[0])\ninterval = pd.Timedelta(swirls.attrs['step_size'])\n\nfor col in range(3):\n    time_index = (col + 1) * 6 - 1\n    timelabel = basetime + pd.Timedelta(interval * (time_index + 1), 'm')\n\n    ax: plt.Axes = fig.add_subplot(\n        gs[0, col],\n        projection=ccrs.PlateCarree()\n    )\n    z = swirls[time_index].values\n    y = swirls[time_index].y\n    x = swirls[time_index].x\n\n    # plot base map\n    plot_base(ax)\n\n    # plot reflectivity\n    ax.contourf(\n        x, y, z, 60,\n        transform=ccrs.PlateCarree(),\n        cmap=cmap, norm=norm, levels=levels\n    )\n\n    # plot motion field\n    ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)\n\n    ax.set_title(\n        f\"Nowcast\\n\" +\n        f\"Initial @ {basetime.strftime('%H:%MZ')}\",\n        loc='left', fontsize=8.75\n    )\n    ax.set_title('')\n    ax.set_title(\n        f\"Initial {basetime.strftime('%Y-%m-%d')} \\n\" +\n        f\"Valid @ {timelabel.strftime('%H:%MZ')} \",\n        loc='right', fontsize=8.75\n    )\n\ncbar_ax = fig.add_axes([0.85, 0.105, 0.01, 0.845])\ncbar = fig.colorbar(\n    scalar_map, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g'\n)\ncbar.ax.set_ylabel('Reflectivity (dBZ)', rotation=90)\n\nfig.savefig(\n    os.path.join(\n        OUTPUT_DIR,\n        \"qpf_laos.png\"\n    ),\n    dpi=450,\n    bbox_inches=\"tight\",\n    pad_inches=0.25\n)\n\nradar_image_time = pd.Timestamp.now()"
      ]
    }
  ],
  "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
}