{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Quantitative Precipitation Estimate (QPE)\n\nThis example demonstrates how to perform QPE,\nusing raingauge data from Manila.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport pandas\nimport numpy as np\nimport pandas as pd\nfrom pyresample import utils\nimport matplotlib\nfrom matplotlib.colors import BoundaryNorm\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\n\n\nimport swirlspy.qpe.rfmap as rfmap\nfrom swirlspy.obs.rain import Rain\nfrom swirlspy.qpe.utils import locate_file, timestamps_ending\n\nTHIS_DIR = os.getcwd()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialising\n\n1. Generate timestamps with timestamps_ending()\n2. Locate rainfall data files with locate_file()\n3. Define target grid as a pyresample AreaDefinition\n4. Construct a Rain object containing rainfall data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dir = THIS_DIR+\"/../tests/samples\"\nlocated_files = []\ntimestamps = timestamps_ending(\n    pd.Timestamp(2018, 8, 11, 00, 40),\n    duration=pd.Timedelta('1 hour'),\n    interval=pd.Timedelta('10 minutes'),\n    exclude_end=False\n    )\n\nfor timestamp in timestamps:\n    located_files.append(locate_file(dir, timestamp))\n\n# Writing AreaDefinition of target grid\narea_id = \"WGS84\"\ndescription = \"Plate Carree\"\nproj_id = \"WGS84\"\nprojection = '+proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees'\nx_size = 500\ny_size = 500\narea_extent = (\n    119.8, 13.6, 122.0, 16.5\n    )\narea_def = utils.get_area_def(\n    area_id, description, proj_id, projection, x_size, y_size, area_extent\n)\n\n# Constructing a Rain object using data from the text file\nrg_object = Rain(\n    located_files,\n    proj=\"Plate Carree\",\n    duration=pd.Timedelta('10 minutes'),\n    NAN=3276.7\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Interpolation\n\nPerform RBF interpolation. Interpolated rainfall is returned\nalong with metadata as an xarray.DataArray.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Interpolate into target grid specified by AreaDefinition\ninterpolated_rf = rfmap.rg_interpolate(\n    rg_object, area_def, 'rbf',\n    coord_label=['Longitude', 'Latitude']\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Generate rainfall map\n\nDefine the colour scale and format of the reflectivity plots.\nGenerate the image using xarray.DataArray.plot().\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Plotting\nlevels = [\n    0, 0.5, 2, 5, 10, 20,\n    30, 40, 50, 70, 100, 150,\n    200, 300, 400, 500, 600, 700\n]\n\ncmap = 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\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\nplt.figure(figsize=(15, 15))\nplt.axis(\"equal\")\n\n# AreaDefinition.to_cartopy_crs() does not work for\n# Geographic Coordinate Systems\nax = plt.axes(projection=ccrs.PlateCarree())\nax.set_extent(\n    [area_extent[0], area_extent[2],\n        area_extent[1], area_extent[3]]\n)\nax.coastlines(resolution='10m')  # plot coastlines\n# Plot data\ninterpolated_rf.plot(\n    cmap=cmap, norm=norm, extend='neither'\n    )\n# Show stations\nrfmap.show_raingauge(rg_object)\n# Plot title\nplt.title(\n    \"Accumulated rainfall from \"\n    f\"{rg_object.start_time} to {rg_object.end_time}\"\n)\n# timestring\nstart_timestr = rg_object.start_time.strftime(\"%Y%m%d%H%M\")\nend_timestr = rg_object.end_time.strftime(\"%Y%m%d%H%M\")\nplt.savefig(\n    os.path.join(\n        THIS_DIR,\n        \"../tests/outputs\",\n        f\"RBF_{start_timestr}_{end_timestr}.png\"\n    )\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
}