{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# QPF (Malaysia)\nThis example demonstrates how to perform\noperational deterministic QPF up to three hours using\nnational radar data.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Definitions\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 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\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\nstart_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialising\n\nThis section demonstrates extracting\nradar reflectivity data.\n\nStep 1: Define your input data directory and output directory\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Supply the directory of radar and nwp data\ndata_dir = os.path.abspath(\n    os.path.join(DATA_DIR, 'netcdf_ms')\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 2: Define a basetime\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Supply basetime\nbasetime = pd.Timestamp('201908090900')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 3: Read data files from the radar data using xarray()\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Radar data listed from the basetime[0] --> 3 hours before the basetime[17] (descending time)\ninterval = 10             # Interval of radar data\nradar_datas = []\nfor i in range(0, 2):\n    t = basetime - pd.Timedelta(minutes=i * interval)\n    # Radar data nomenclature\n    filename = os.path.join(\n        data_dir,\n        t.strftime(\"radar_d03_%Y-%m-%d_%H_%M_00.rapids.nc\")\n    )\n    reflec = xr.open_dataset(filename)\n    radar_datas.append(reflec)\n\n# Concatenate list by time\nreflec_concat = xr.concat(radar_datas, dim='time')\n\n# Extracting the radar data: The radar dBZ variable is named 'Zradar', therefore, we extract 'Zradar'\nradar = reflec_concat['Zradar']\n\n# Reversing such that time goes from earliest to latest; 3 hours before basetime[0] --> basetime[17]\nradar = radar.sortby('time', ascending=True)\n\n# Filtering\nradar = radar.where(radar > 15, np.nan)\n\ninitialising_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Nowcast (SWIRLS-Radar-Advection)\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 <9999.> frame_type <FrameType.dBZ>\nstandardize_attr(radar, frame_type=FrameType.dBZ, zero_value=np.nan)\n\n# Rover motion field computation\nmotion = rover(radar)\nrover_time = pd.Timestamp.now()\n\n# Semi-Lagrangian Advection\nswirls = sla(radar, motion, 18)  # Radar time goes from earliest to latest\nsla_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plotting result\nStep 1: Defining the dBZ levels, colorbar parameters and projection\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    '#FFFFFF', '#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, 120, 0.5, 7.25]\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 2: 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=(5 + 1, 3 * 2),\n    frameon=False\n)\ngs = GridSpec(\n    3, 1, figure=fig,\n    wspace=0.03, hspace=0, top=0.95, bottom=0.05, left=0.17, right=0.845\n)\n\nfor row in range(3):\n    time_index = (row + 1) * 6\n    timelabel = basetime + pd.Timedelta(interval * (time_index), 'm')\n\n    ax: plt.Axes = fig.add_subplot(\n        gs[row, 0],\n        projection=ccrs.PlateCarree()\n    )\n    z = swirls[time_index].values\n    lats = swirls[time_index].latitude\n    lons = swirls[time_index].longitude\n    title = 'SWIRLS Reflectivity'\n\n    # plot base map\n    plot_base(ax)\n\n    # plot reflectivity\n    ax.contourf(\n        lons, lats, z, 60,\n        transform=ccrs.PlateCarree(),\n        cmap=cmap, norm=norm, levels=levels\n    )\n\n    ax.set_title(\n        f\"{title}\\n\" +\n        f\"Initial @ {basetime.strftime('%H:%MZ')}\",\n        loc='left', fontsize=9\n    )\n    ax.set_title('')\n    ax.set_title(\n        f\"Initial {basetime.strftime('%Y-%m-%d')} \\n\" +\n        f\"Forecast Valid @ {timelabel.strftime('%H:%MZ')} \",\n        loc='right', fontsize=9\n    )\n\ncbar_ax = fig.add_axes([0.9, 0.105, 0.04, 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        \"swirls_ms_fcs.png\"\n    ),\n    dpi=450,\n    bbox_inches=\"tight\",\n    pad_inches=0.1\n)\n\nradar_image_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\"SLA time: {sla_time}\")\nprint(f\"Plotting radar image time: {radar_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 plot radar image: {radar_image_time - sla_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
}