{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image overlay\n",
    "\n",
    "This notebook shows how you can overlay an image on a Leaflet map, and use a slider to walk through a time series. This is similar to the Video overlay example, but offers more control as you can do what you like with the slider, including going backwards.\n",
    "\n",
    "The following libraries are needed:\n",
    "* `tqdm`\n",
    "* `rasterio`\n",
    "* `matplotlib`\n",
    "* `pandas`\n",
    "* `ipyleaflet`\n",
    "\n",
    "The recommended way is to `conda install -c conda-forge` them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ftplib\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "import rasterio\n",
    "from rasterio.warp import reproject, Resampling\n",
    "from affine import Affine\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import subprocess\n",
    "import datetime\n",
    "import asyncio\n",
    "from time import time\n",
    "from base64 import b64encode\n",
    "from ipyleaflet import Map, ImageOverlay, basemaps\n",
    "from ipywidgets import VBox, Checkbox, SelectionRangeSlider\n",
    "try:\n",
    "    from StringIO import StringIO\n",
    "    py3 = False\n",
    "except ImportError:\n",
    "    from io import StringIO, BytesIO\n",
    "    py3 = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we download some satellite rainfall estimates ([NASA's Global Precipitation Measurement](https://www.nasa.gov/mission_pages/GPM/main/index.html)). They come in the WGS84 projection and cover the entire globe between latitudes 60N and 60S. They have a spatial resolution of 0.1° and a temporal resolution of 30 minutes. We first download each individual files for the day 2017/01/01 (48 files)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs('tif', exist_ok=True)\n",
    "ftp = ftplib.FTP('arthurhou.pps.eosdis.nasa.gov')\n",
    "passwd = 'ebwje/cspdibsuAhnbjm/dpn'\n",
    "login = ''\n",
    "for c in passwd:\n",
    "    login += chr(ord(c) - 1)\n",
    "ftp.login(login, login)\n",
    "ftp.cwd('gpmdata/2017/01/01/gis')\n",
    "lst = [i for i in ftp.nlst() if i.startswith('3B-HHR-GIS.MS.MRG.3IMERG.') and i.endswith('.V06B.tif')]\n",
    "for filename in tqdm(lst):\n",
    "    if not os.path.exists('tif/' + filename):\n",
    "        ftp.retrbinary(\"RETR \" + filename, open('tif/' + filename, 'wb').write)\n",
    "ftp.quit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we convert them into images, we will need to apply a color map to the data. In order to scale this color map, we need to extract the maximum value present in the whole data set. Since the maximum value will appear very rarely, the visual rendering will be better if we saturate the images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_vmax = 0\n",
    "for f in os.listdir('tif'):\n",
    "    dataset = rasterio.open('tif/' + f)\n",
    "    data = dataset.read()[0][300:1500]\n",
    "    data = np.where(data==9999, np.nan, data) / 10 # in mm/h\n",
    "    image_vmax = max(image_vmax, np.nanmax(data))\n",
    "image_vmax *= 0.1 # saturate images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We reproject our data (originally in WGS84, also known as EPSG:4326) to Web Mercator (also known as EPSG:3857), which is the projection used by Leaflet. After applying a color map (here `plt.cm.jet`), we save each image to a PNG file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# At this point if GDAL complains about not being able to open EPSG support file gcs.csv, try in the terminal:\n",
    "# export GDAL_DATA=`gdal-config --datadir`\n",
    "\n",
    "def to_web_mercator(data):\n",
    "    with rasterio.Env():\n",
    "        rows, cols = 1200, 3600\n",
    "        src_transform = Affine(0.1, 0, -180, 0, -0.1, 60)\n",
    "        src_crs = {'init': 'EPSG:4326'}\n",
    "        source = np.where((data==9999) | (~np.isfinite(data)), 0, data) / 10 # in mm/h\n",
    "        dst_crs = {'init': 'EPSG:3857'}\n",
    "        bounds = [-180, -60, 180, 60]\n",
    "        dst_transform, width, height = rasterio.warp.calculate_default_transform(src_crs, dst_crs, cols, rows, *bounds)\n",
    "        dst_shape = height, width\n",
    "\n",
    "        destination = np.zeros(dst_shape)\n",
    "\n",
    "        reproject(\n",
    "            source,\n",
    "            destination,\n",
    "            src_transform=src_transform,\n",
    "            src_crs=src_crs,\n",
    "            dst_transform=dst_transform,\n",
    "            dst_crs=dst_crs,\n",
    "            resampling=Resampling.nearest)\n",
    "\n",
    "    return destination\n",
    "\n",
    "def to_png(data, vmax, fname):\n",
    "    fig, ax = plt.subplots(1, figsize=(36, 12))\n",
    "    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)\n",
    "    ax.imshow(data, vmin=0, vmax=vmax, interpolation='nearest', cmap=plt.cm.jet)\n",
    "    ax.axis('tight')\n",
    "    ax.axis('off')\n",
    "    plt.savefig(fname)\n",
    "    plt.close()\n",
    "\n",
    "os.makedirs('png', exist_ok=True)\n",
    "\n",
    "for fname in tqdm(os.listdir('tif')):\n",
    "    png_name = fname[:-3] + 'png'\n",
    "    if not os.path.exists('png/' + png_name):\n",
    "        dataset = rasterio.open('tif/' + fname)\n",
    "        data = dataset.read()[0][300:1500]\n",
    "        data_web = to_web_mercator(data)\n",
    "        to_png(data_web, image_vmax, 'png/' + png_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create a map, a slider and a check box. The check box allows to walk through the images as snapshots (the left and right boundaries of the slider are tied together). When unchecked, you can move the left and right boundaries as you want and the mean precipitation in that range will be displayed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = [0, 0]\n",
    "zoom = 2\n",
    "m = Map(center=center, zoom=zoom, interpolation='nearest', basemap=basemaps.CartoDB.DarkMatter)\n",
    "\n",
    "bounds = [(-60, -180), (60, 180)]\n",
    "start_date = datetime.datetime(2017, 1, 1)\n",
    "end_date = datetime.datetime(2017, 1, 2)\n",
    "\n",
    "dates = pd.date_range(start_date, end_date, freq='30min')\n",
    "options = [(date.strftime(' %d:%H:%M '), date) for date in dates]\n",
    "index = (0, 1)\n",
    "date_slider = SelectionRangeSlider(\n",
    "    options=options,\n",
    "    index=index,\n",
    "    description='Time',\n",
    "    orientation='horizontal',\n",
    "    layout={'width': '900px'}\n",
    ")\n",
    "\n",
    "snapshot = Checkbox(\n",
    "    value=True,\n",
    "    description='Snapshot',\n",
    "    disabled=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computing of the mean precipitation can take some time, so we don't want to trigger it every time the slider moves. The following code implements a debouncer which ensures that a new computation is triggered only when the user stops moving the slider. There is also a throttler which ensures that snapshot images are not displayed too often, which gives a smoother animation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Timer:\n",
    "    def __init__(self, timeout, callback):\n",
    "        self._timeout = timeout\n",
    "        self._callback = callback\n",
    "        self._task = asyncio.ensure_future(self._job())\n",
    "\n",
    "    async def _job(self):\n",
    "        await asyncio.sleep(self._timeout)\n",
    "        self._callback()\n",
    "\n",
    "    def cancel(self):\n",
    "        self._task.cancel()\n",
    "\n",
    "def debounce(wait):\n",
    "    \"\"\" Decorator that will postpone a function's\n",
    "        execution until after `wait` seconds\n",
    "        have elapsed since the last time it was invoked. \"\"\"\n",
    "    def decorator(fn):\n",
    "        timer = None\n",
    "        def debounced(*args, **kwargs):\n",
    "            nonlocal timer\n",
    "            def call_it():\n",
    "                fn(*args, **kwargs)\n",
    "            if timer is not None:\n",
    "                timer.cancel()\n",
    "            timer = Timer(wait, call_it)\n",
    "        return debounced\n",
    "    return decorator\n",
    "\n",
    "def throttle(wait):\n",
    "    \"\"\" Decorator that prevents a function from being called\n",
    "        more than once every wait period. \"\"\"\n",
    "    def decorator(fn):\n",
    "        time_of_last_call = 0\n",
    "        scheduled = False\n",
    "        new_args, new_kwargs = None, None\n",
    "        def throttled(*args, **kwargs):\n",
    "            nonlocal new_args, new_kwargs, time_of_last_call, scheduled\n",
    "            def call_it():\n",
    "                nonlocal new_args, new_kwargs, time_of_last_call, scheduled\n",
    "                time_of_last_call = time()\n",
    "                fn(*new_args, **new_kwargs)\n",
    "                scheduled = False\n",
    "            time_since_last_call = time() - time_of_last_call\n",
    "            new_args = args\n",
    "            new_kwargs = kwargs\n",
    "            if not scheduled:\n",
    "                new_wait = max(0, wait - time_since_last_call)\n",
    "                Timer(new_wait, call_it)\n",
    "                scheduled = True\n",
    "        return throttled\n",
    "    return decorator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the logic for displaying a new image when the slider moves. The images are sent to the browser by embedding the data into the URL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snapshot_fname = None\n",
    "io = None\n",
    "\n",
    "def get_tif_name(date):\n",
    "    hour = str(date.hour).zfill(2)\n",
    "    minu1 = str(date.minute).zfill(2)\n",
    "    minu2 = str(date.minute + 29).zfill(2)\n",
    "    minu3 = str(date.hour * 60 + date.minute).zfill(4)\n",
    "    fname = '3B-HHR-GIS.MS.MRG.3IMERG.20170101-S' + hour + minu1 + '00-E' + hour + minu2 + '59.' + minu3 + '.V06B.tif'\n",
    "    return fname\n",
    "\n",
    "def display_png(fname):\n",
    "    global io\n",
    "    with open(fname, 'rb') as f:\n",
    "        data = b64encode(f.read())\n",
    "    if py3:\n",
    "        data = data.decode('ascii')\n",
    "    imgurl = 'data:image/png;base64,' + data\n",
    "    if io is None:\n",
    "        io = ImageOverlay(url=imgurl, bounds=bounds, opacity=0.4)\n",
    "        m.add_layer(io)\n",
    "    else:\n",
    "        io.url = imgurl\n",
    "\n",
    "@debounce(0.2)\n",
    "def show_mean():\n",
    "    mean_data = None\n",
    "    mean_vmax = 0\n",
    "    nb = 0\n",
    "    date = date_slider.value[0]\n",
    "    while date < date_slider.value[1]:\n",
    "        nb += 1\n",
    "        fname = get_tif_name(date)\n",
    "        dataset = rasterio.open('tif/' + fname)\n",
    "        data = dataset.read()[0][300:1500]\n",
    "        data = np.where(data==9999, np.nan, data) / 10 # in mm/h\n",
    "        mean_vmax = max(mean_vmax, np.nanmax(data))\n",
    "        if mean_data is None:\n",
    "            mean_data = data\n",
    "        else:\n",
    "            mean_data = mean_data + data\n",
    "        date += datetime.timedelta(minutes=30)\n",
    "    mean_data = mean_data / nb\n",
    "    mean_vmax *= 0.1 / nb # saturate image\n",
    "    data_web = to_web_mercator(mean_data)\n",
    "    to_png(data_web, mean_vmax, 'png/mean.png')\n",
    "    display_png('png/mean.png')\n",
    "\n",
    "@throttle(0.1)\n",
    "def show_snapshot(change):\n",
    "    global snapshot_fname\n",
    "    if change['new'][0] == end_date:\n",
    "        date_slider.value = [end_date - datetime.timedelta(minutes=30), end_date]\n",
    "    elif change['new'][1] == start_date:\n",
    "        date_slider.value = [start_date, start_date + datetime.timedelta(minutes=30)]\n",
    "    elif change['new'][0] != change['old'][0]:\n",
    "        # left has changed\n",
    "        date_slider.value = [change['new'][0], change['new'][0] + datetime.timedelta(minutes=30)]\n",
    "    else:\n",
    "        # right has changed    \n",
    "        date_slider.value = [change['new'][1] - datetime.timedelta(minutes=30), change['new'][1]]\n",
    "    new_fname = get_tif_name(date_slider.value[0])[:-3] + 'png'\n",
    "    if new_fname != snapshot_fname:\n",
    "        display_png('png/' + new_fname)\n",
    "        snapshot_fname = new_fname\n",
    "\n",
    "def date_changed(change):\n",
    "    if snapshot.value:\n",
    "        show_snapshot(change)\n",
    "        \n",
    "    else:\n",
    "        show_mean()\n",
    "\n",
    "date_slider.observe(date_changed, 'value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "VBox([date_slider, snapshot, m])"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}