{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Video overlay\n",
    "\n",
    "This notebook shows how you can overlay a video on a Leaflet map. This video can come off-the-shelf from an URL, but we will also see how to create your own video from NumPy arrays.\n",
    "\n",
    "The following libraries are needed:\n",
    "* `tqdm`\n",
    "* `rasterio`\n",
    "* `numpy`\n",
    "* `matplotlib`\n",
    "* `ipyleaflet`\n",
    "\n",
    "The recommended way is to try to `conda install` them first, and if they are not found then `pip install`."
   ]
  },
  {
   "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 subprocess\n",
    "from base64 import b64encode\n",
    "from ipyleaflet import Map, VideoOverlay\n",
    "try:\n",
    "    from StringIO import StringIO\n",
    "    py3 = False\n",
    "except ImportError:\n",
    "    from io import StringIO, BytesIO\n",
    "    py3 = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = [25, -130]\n",
    "zoom = 4\n",
    "m = Map(center=center, zoom=zoom)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you happen to find a pre-built video for which you know the geographic bounds, overlaying it is as simple as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = [(13, -130), (32, -100)]\n",
    "io = VideoOverlay(url='https://www.mapbox.com/bites/00188/patricia_nasa.webm', bounds=bounds)\n",
    "m.add_layer(io)\n",
    "io.interact(opacity=(0.0,1.0,0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However if you zoom in and play with the opacity, you will see that the coastline in the video does not overlay perfectly with the Leaflet map. This is because Leaflet uses a projection known as Web Mercator, whereas the video comes from satellite imagery. The view of the Earth from a geostationary satellite is called a disk, which is the part of the Earth that the satellite sees. Usually the disk is reprojected to e.g. WGS84, which is slightly different from Web Mercator. On a relatively small area, like in this video, the difference will be acceptable, but if you want to show some data over the whole map, it can be a problem. In that case you will need to reproject, as we will see.\n",
    "\n",
    "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('.V05B.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": [
    "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",
    "    vmax = max(vmax, np.nanmax(data))\n",
    "vmax *= 0.5 # saturate a little bit"
   ]
  },
  {
   "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",
    "os.makedirs('png', exist_ok=True)\n",
    "\n",
    "for f in tqdm(os.listdir('tif')):\n",
    "    png_name = f[:-3] + 'png'\n",
    "    if not os.path.exists('png/' + png_name):\n",
    "        dataset = rasterio.open('tif/' + f)\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",
    "            data = dataset.read()[0][300:1500]\n",
    "            source = np.where((data==9999) | (~np.isfinite(data)), 0, data) / 10 # in mm/h\n",
    "\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",
    "        data_web = destination\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_web, vmin=0, vmax=vmax, interpolation='nearest', cmap=plt.cm.jet)\n",
    "        ax.axis('tight')\n",
    "        ax.axis('off')\n",
    "        plt.savefig('png/' + png_name)\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use `ffmpeg` to create a video from our individual images ([ffmpeg](https://www.ffmpeg.org) can be installed on many systems). This utility needs our files to be named with a sequential number, which is why we rename them. The `ffmpeg` commands are pretty obscure but they do the job, and we finally get a `rain.mp4` video!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_files = os.listdir('png')\n",
    "png_files.sort()\n",
    "for i, f in enumerate(png_files):\n",
    "    if f.startswith('3B-HHR-GIS.MS.MRG.3IMERG.'):\n",
    "        os.rename('png/' + f, 'png/f' + str(i).zfill(2) + '.png')\n",
    "\n",
    "if not os.path.exists('mp4/rain.mp4'):\n",
    "    os.makedirs('mp4', exist_ok=True)\n",
    "    bitrate = '4000k'\n",
    "    framerate = '12'\n",
    "    cmd = 'ffmpeg -r $framerate -y -f image2 -pattern_type glob -i \"png/*.png\" -c:v libx264 -preset slow -b:v $bitrate -pass 1 -c:a libfdk_aac -b:a 0k -f mp4 -r $framerate -profile:v high -level 4.2 -pix_fmt yuv420p -movflags +faststart -vf \"scale=trunc(iw/2)*2:trunc(ih/2)*2\" /dev/null && \\\n",
    "    ffmpeg -r $framerate -f image2 -pattern_type glob -i \"png/*.png\" -c:v libx264 -preset slow -b:v $bitrate -pass 2 -c:a libfdk_aac -b:a 0k -r $framerate -profile:v high -level 4.2 -pix_fmt yuv420p -movflags +faststart -vf \"scale=trunc(iw/2)*2:trunc(ih/2)*2\" mp4/rain.mp4'\n",
    "    cmd = cmd.replace('$framerate', framerate).replace('$bitrate', bitrate)\n",
    "    subprocess.check_output(cmd, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The video can be sent to the browser by embedding the data into the URL, et voilà!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "center = [0, -70]\n",
    "zoom = 3\n",
    "m = Map(center=center, zoom=zoom, interpolation='nearest')\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if py3:\n",
    "    f = BytesIO()\n",
    "else:\n",
    "    f = StringIO()\n",
    "with open('mp4/rain.mp4', 'rb') as f:\n",
    "    data = b64encode(f.read())\n",
    "if py3:\n",
    "    data = data.decode('ascii')\n",
    "videourl = 'data:video/mp4;base64,' + data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = [(-60, -180), (60, 180)]\n",
    "io = VideoOverlay(url=videourl, bounds=bounds)\n",
    "m.add_layer(io)\n",
    "io.interact(opacity=(0.0,1.0,0.01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}