{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nRadar QPE (Quantitative Precipitation Estimate)\n===============================================\nThis example demonstrates how to generate radar QPE\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport xarray as xr\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom pyresample import utils\n\nfrom matplotlib.colors import BoundaryNorm\n\nfrom swirlspy.rad.iris import read_iris\nfrom swirlspy.qpe.utils import dbz2rr, rr2rf, timestamps_ending, locate_file\n\nplt.switch_backend('agg')\n\nFILE_LOC = os.getcwd()\nTHIS_DIR = FILE_LOC+'/../tests/int'\nos.chdir(THIS_DIR)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read IRIS RAW data files into xarray\n---------------------------------------------------\n\n1. Generate timestamps with timestamps_ending()\n2. Locate radar data files with locate_file()\n3. Define target grid as a pyresample AreaDefinition\n4. Read radar data files with read_iris()\n5. Convert reflectivity in dBZ to rainrate in mm/h with dbz2rr()\n6. Convert rainrate to rainfall in 6 mins with rr2rf()\n7. Concatenate all xarrays into one with xarray.concat()\n8. Compute the sum along the time dimension\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dir = THIS_DIR+\"/../samples/iris\"\nlocated_files = []\ntimestamps = timestamps_ending(pd.Timestamp(2018, 8, 22, 19, 00))\nfor timestamp in timestamps:\n    located_files.append(locate_file(dir, timestamp))\n\narea_id = \"hk1980_250km\"\ndescription = (\n    \"A 250 m resolution rectangular grid \"\n    \"centred at HKO and extending to 250 km \"\n    \"in each direction in HK1980 easting/northing coordinates\"\n    )\nproj_id = 'hk1980'\nprojection = (\n    '+proj=tmerc +lat_0=22.31213333333334 '\n    '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '\n    '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'\n    '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '\n    '+no_defs'\n    )\nx_size = 500\ny_size = 500\narea_extent = (587000, 569000, 1087000, 1069000)\narea_def = utils.get_area_def(\n    area_id, description, proj_id, projection, x_size, y_size, area_extent\n)\n\nrainfalls = []\nfor filename in located_files:\n    reflectivity = read_iris(\n        filename,\n        area_def=area_def,\n        coord_label=['easting', 'northing'],\n        radar_radius=250000\n    )\n    rainrate = dbz2rr(reflectivity, a=58.53, b=1.56)\n    rainfall = rr2rf(rainrate, time_shift=pd.Timedelta(minutes=6))\n    rainfalls.append(rainfall)\n\nrainfall_cat = xr.concat(rainfalls, dim='time')\nrainfall_hourly = rainfall_cat.sum(dim='time')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generate rainfall map\n---------------------------------------------------\nDefine colour scale and format of the rainfall map plot\nthen generate it with xarray.plot()\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "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\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, 10))\nplt.axis('equal')\nrainfall_hourly.plot(cmap=cmap, norm=norm)\nplt.savefig(\"radar_qpe.png\")"
      ]
    }
  ],
  "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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}