{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nBlending\n========================================================\nThis example demonstrates how to blend different \nreflectivity sources into one.\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 time calculations\nimport pandas as pd\n# Python package for numerical calculations\nimport numpy as np\n# Python package for projection\nimport cartopy.crs as ccrs\n# Python package for land/sea features\nimport cartopy.feature as cfeature\n# Python package for reading map shape file\nimport cartopy.io.shapereader as shpreader\n# Python package for creating plots\nfrom matplotlib import pyplot as plt\n# Python package for output grid format\nfrom matplotlib.gridspec import GridSpec\n# Python package for colorbars\nfrom matplotlib.colors import BoundaryNorm, ListedColormap\nfrom matplotlib.cm import ScalarMappable\n# Python package for projection description\nfrom pyresample import get_area_def\n\n# swirlspy regrid function\nfrom swirlspy.core.resample import grid_resample\n# swirlspy iris parser function\nfrom swirlspy.rad.iris import read_iris_grid\n# swirlspy h8 parser function\nfrom swirlspy.sat.h8 import read_h8_data\n# swirlspy blending function\nfrom swirlspy.blending import comp_qpe, Raw\n# directory constants\nfrom swirlspy.tests.samples import DATA_DIR\nfrom swirlspy.tests.outputs import OUTPUT_DIR\n\nwarnings.filterwarnings(\"ignore\")\nplt.switch_backend('agg')\n\nstart_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialising\n---------------------------------------------------\n\nThis section demonstrates parsing\nHimawari-8 data.\n\nStep 1: Define necessary parameter.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Define base time\nbase_time = pd.Timestamp(\"2019-07-31T07:00\")\n\n# Define data boundary in WGS84 (latitude)\nlatitude_from = 30.\nlatitude_to = 16.\nlongitude_from = 105.\nlongitude_to = 122.\n\narea = (\n    latitude_from, latitude_to,\n    longitude_from, longitude_to\n)\n\n# Define grid size, use negative value for descending range\ngrid_size = (-.025, .025)\n\n# list of source data\nsources = []\n\ninitialising_time = pd.Timestamp.now()\n\n# Load map shape\nmap_shape_file = os.path.abspath(os.path.join(DATA_DIR, 'shape/se_asia'))\n\n# coastline and province\nmap_with_province = cfeature.ShapelyFeature(\n    list(shpreader.Reader(map_shape_file).geometries()),\n    ccrs.PlateCarree()\n)\n\n\ndef plot_base(ax: plt.Axes, extents: list, crs: ccrs.Projection):\n    \"\"\"\n    base map function\n    \"\"\"\n    ax.set_extent(extents, crs=crs)\n    # fake the ocean color\n    ax.imshow(\n        np.tile(np.array([[[178, 208, 254]]], dtype=np.uint8), [2, 2, 1]),\n        origin='upper', transform=ccrs.PlateCarree(),\n        extent=[-180, 180, -180, 180], zorder=-1\n    )\n    # coastline, province and state, color\n    ax.add_feature(\n        map_with_province, facecolor=cfeature.COLORS['land'],\n        edgecolor='none', zorder=0\n    )\n    # overlay coastline, province and state without color\n    ax.add_feature(\n        map_with_province, facecolor='none', edgecolor='gray', linewidth=0.5\n    )\n    ax.set_title('')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 2: Read data from radar files into xarray.DataArray\nusing read_iris_grid().\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radar = read_iris_grid(os.path.join(DATA_DIR, \"iris/RAD190731150000.REF2256\"))\n\nradar_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 3: Define the target grid as a pyresample AreaDefinition.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Defining target grid\narea_id = \"WGS84\"\ndescription = 'World Geodetic System 1984'\nproj_id = 'WGS84'\nprojection = '+proj=longlat +datum=WGS84 +no_defs'\nx_size = (longitude_to - longitude_from) / grid_size[1] + 1\ny_size = (latitude_to - latitude_from) / grid_size[0] + 1\narea_extent = (longitude_from, latitude_from, longitude_to, latitude_to)\nradar_area_def = get_area_def(\n    area_id, description, proj_id, projection, x_size, y_size, area_extent\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 5: Reproject the radar data from read_iris_grid() from Centered\nAzimuthal (source) projection to World Geodetic System 1984 projection.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Extracting the AreaDefinition of the source projection\narea_def_src = radar.attrs['area_def']\n\n# Reprojecting\nreproj_radar = grid_resample(\n    radar, area_def_src, radar_area_def,\n    coord_label=['x', 'y']\n).sortby(\n    ['y'], ascending=False\n)\n\n# fix floating point issue\ny_coords = np.linspace(\n    latitude_from,\n    latitude_to,\n    reproj_radar.data.shape[1],\n    dtype=np.float32\n)\nx_coords = np.linspace(\n    longitude_from,\n    longitude_to,\n    reproj_radar.data.shape[2],\n    dtype=np.float32\n)\nreproj_radar.coords['y'] = np.array(y_coords)\nreproj_radar.coords['x'] = np.array(x_coords)\nreproj_radar = reproj_radar.sel(time=reproj_radar.coords['time'].values[0])\n\n\nradar_site = (\n    reproj_radar.attrs['proj_site'][1],\n    reproj_radar.attrs['proj_site'][0],\n    1.8,  # radius\n    0.76  # weight sigma\n)\n\nsources.append(Raw(\n    reproj_radar,\n    [radar_site],  # sites configuration, list of available sites useful for mosaic data\n    0.1  # data weight\n))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 6: Define data directory\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Supply data directory.\n# Please make sure H8 data filename is follow the naming pattern -\n# HS_H08_{date}_{time}_B{channel:02}_FLDK_R{rsol:02}_S{seg:02}10.DAT\n# example:\n#   base time = 2019-07-31 07:00 UTC\n#   channel = 4\n#   resolution = 10\n#   segment = 2\n#   ========================================\n#   filename: HS_H08_20190731_0700_B04_FLDK_R10_S0410.DAT\ndata_dir = os.path.join(DATA_DIR, \"h8\")\n\nsat_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 7: Parse data into reflectivity as xarray.DataArray\nusing read_h8_data().\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sat = read_h8_data(\n    data_dir,\n    base_time,\n    area,\n    grid_size\n)\n# remove time axis\nsat = sat.sel(time=sat.coords['time'].values[0])\n\n# no site data used, treat all points of data with same weight\nsources.append(Raw(\n    sat,\n    weight=0.01  # data weight\n))\n\nblend_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 8: Blend all data together.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reflec = comp_qpe(\n    grid_size,\n    area,\n    sources\n)\n\n\npost_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Step 9: Remove invalid data if needed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reflec.values[reflec.values < 13.] = reflec.attrs['zero_value']\n\n# update sat data for plotting\nsat.values[sat.values < 13.] = sat.attrs['zero_value']"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Generating radar reflectivity maps\n-----------------------------------\n\nDefine the color scale and format of the plots\nand plot using xarray.plot().\n\nIn this example, only hourly images will be plotted.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Defining colour scale and format\nlevels = [\n    -32768,\n    10, 15, 20, 24, 28, 32,\n    34, 38, 41, 44, 47, 50,\n    53, 56, 58, 60, 62\n]\ncmap = ListedColormap([\n    '#FFFFFF00', '#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# colorbar map\nmappable = ScalarMappable(cmap=cmap, norm=norm)\nmappable.set_array([])\n\n# Defining the crs\ncrs = ccrs.PlateCarree()\nextents = (longitude_from, longitude_to, latitude_from, latitude_to)\n\n# Plotting\nfig: plt.Figure = plt.figure(figsize=(24, 8), frameon=False)\ngs = GridSpec(\n    1, 3, figure=fig, wspace=0.03, hspace=-0.25, top=0.95,\n    bottom=0.05, left=0.17, right=0.845\n)\n\n# plot radar\nax = fig.add_subplot(gs[0, 0], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(reproj_radar.values, cmap=cmap, norm=norm, interpolation='nearest',\n               extent=extents)\n\nax.set_title(\n    \"RADAR\\n\"\n    f\"Based @ {base_time.strftime('%H:%MH')}\",\n    loc='left',\n    fontsize=9\n)\nax.set_title(\n    ''\n)\nax.set_title(\n    f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n    f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n    loc='right',\n    fontsize=9\n)\n\n# plot H8\nax = fig.add_subplot(gs[0, 1], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(sat.values, cmap=cmap, norm=norm, interpolation='nearest',\n               extent=extents)\n\nax.set_title(\n    \"H8\\n\"\n    f\"Based @ {base_time.strftime('%H:%MH')}\",\n    loc='left',\n    fontsize=9\n)\nax.set_title(\n    ''\n)\nax.set_title(\n    f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n    f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n    loc='right',\n    fontsize=9\n)\n\n# plot blended\nax = fig.add_subplot(gs[0, 2], projection=crs)\n\nplot_base(ax, extents, crs)\n\nim = ax.imshow(reflec.values, cmap=cmap, norm=norm, interpolation='nearest',\n               extent=extents)\n\nax.set_title(\n    \"Blended\\n\"\n    f\"Based @ {base_time.strftime('%H:%MH')}\",\n    loc='left',\n    fontsize=9\n)\nax.set_title(\n    ''\n)\nax.set_title(\n    f\"{base_time.strftime('%Y-%m-%d')} \\n\"\n    f\"Valid @ {(base_time + pd.Timedelta(minutes=10)).strftime('%H:%MH')} \",\n    loc='right',\n    fontsize=9\n)\n\ncbar_ax = fig.add_axes([0.875, 0.125, 0.03, 0.75])\ncbar = fig.colorbar(\n    mappable, cax=cbar_ax, ticks=levels[1:], extend='max', format='%.3g')\ncbar.ax.set_ylabel('Reflectivity', rotation=90)\n\nfig.savefig(\n    os.path.join(OUTPUT_DIR, \"blending.png\"),\n    bbox_inches='tight'\n)\n\nimage_time = pd.Timestamp.now()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Checking run time of each component\n--------------------------------------------------------------------\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\"Read radar time: {radar_time}\")\nprint(f\"Parse H8 data time: {sat_time}\")\nprint(f\"Blending time: {blend_time}\")\nprint(f\"Post blending time: {post_time}\")\nprint(f\"Plotting blended image time: {image_time}\")\n\nprint(f\"Time to initialise: {initialising_time - start_time}\")\nprint(f\"Time to run read radar: {radar_time - initialising_time}\")\nprint(f\"Time to run data parsing: {sat_time - radar_time}\")\nprint(f\"Time to run blending: {blend_time - sat_time}\")\nprint(f\"Time to perform post process: {post_time - blend_time}\")\nprint(f\"Time to plot reflectivity image: {image_time - post_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
}