{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nQuantitative Precipitation Forecast (QPF)\n====================================================\nThis example demonstrates how to perform QPF using radar data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import xarray as xr\nimport os.path\nimport numpy as np\nimport os\nfrom pyresample import utils\nfrom matplotlib.colors import BoundaryNorm\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nfrom swirlspy.qpf.rover import rover\nfrom swirlspy.rad.iris import read_iris\n\ntry:\n    THIS_DIR = os.path.dirname(os.path.abspath(__file__))\nexcept NameError:\n    import sys\n    THIS_DIR = os.path.dirname(os.path.abspath(sys.argv[0]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialising\n---------------------------------------------------\n1. Define target grid as a pyresample AreaDefinition.\n2. Read radar data files with read_iris().\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Defining target grid as a pyresample AreaDefinition.\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\n# Read data from files (in this example, the files are\n# in iris format) and storing data as xarrays.\n\nvalid_time = '201805150624'\nsteps = 6\ntime_step_size = 6\n\nfile_name1 = \"../tests/samples/TMS170524110001.RAW2SGJ\"\nfile_name2 = \"../tests/samples/TMS170524110619.RAW2SJY\"\n\n# Reading file data, storing as xarray\nxarray1 = read_iris(\n    file_name1,\n    area_def=area_def,\n    coord_label=['easting', 'northing'],\n    radar_radius=250000\n)\nxarray2 = read_iris(\n    file_name2,\n    area_def=area_def,\n    coord_label=['easting', 'northing'],\n    radar_radius=250000\n)\nqpf_xarray = xr.concat([xarray1, xarray2], dim='time')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Extrapolation\n------------------------------\nPerform rover extrapolation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Calling rover\nmotion_u, motion_v, reflectivity = rover(\n    qpf_xarray, '201705241100', 6, 6\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generating labelled radar colour images\n------------------------\n\nDefine the colour scale and the format of the\nreflectivity plots. Then generate using xarray.plot().\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "levels = [\n    -32768,\n    10, 15, 20, 24, 28, 32,\n    34, 38, 41, 44, 47, 50,\n    53, 56, 58, 60, 62\n]\ncmap = matplotlib.colors.ListedColormap([\n    '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',\n    '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',\n    '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'\n])\n\nnorm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n\n# Plot radar colour images\nfor i in range(1, len(reflectivity)):\n    plt.figure(figsize=(15, 10))\n    plt.axis(\"equal\")\n    reflectivity[i].plot(cmap=cmap, norm=norm)\n    plt.savefig(\"rover-output-\"+\"{:02d}\".format(i)+\".png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generating radar reflectivity maps\n----------\nUsing the colour scale generated in the previous section,\ngenerate radar reflectivity maps with coastlines.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for j in range(1, len(reflectivity)):\n    plt.figure(figsize=(15, 10))\n    crs = area_def.to_cartopy_crs()\n    ax = plt.axes(projection=crs)\n    ax.coastlines(resolution='10m')\n    plt.axis(\"equal\")\n    reflectivity[j].plot(cmap=cmap, norm=norm)\n    plt.savefig(\"rover-output-map-\"+\"{:02d}\".format(j)+\".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
}