{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From NumPy to Leaflet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to display some raster geographic data in IPyLeaflet. The data is a NumPy array, which means that you have all the power of the Python scientific stack at your disposal to process it.\n",
    "\n",
    "The following libraries are needed:\n",
    "* `requests`\n",
    "* `tqdm`\n",
    "* `rasterio`\n",
    "* `numpy`\n",
    "* `scipy`\n",
    "* `pillow`\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 requests\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "import zipfile\n",
    "import rasterio\n",
    "from affine import Affine\n",
    "import numpy as np\n",
    "import scipy.ndimage\n",
    "from rasterio.warp import reproject, RESAMPLING\n",
    "import PIL\n",
    "import matplotlib.pyplot as plt\n",
    "from base64 import b64encode\n",
    "try:\n",
    "    from StringIO import StringIO\n",
    "    py3 = False\n",
    "except ImportError:\n",
    "    from io import StringIO, BytesIO\n",
    "    py3 = True\n",
    "from ipyleaflet import Map, ImageOverlay, basemap_to_tiles, basemaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download a raster file representing the flow accumulation for South America. This gives an idea of the river network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'http://earlywarning.usgs.gov/hydrodata/sa_30s_zip_grid/sa_acc_30s_grid.zip'\n",
    "filename = os.path.basename(url)\n",
    "name = filename[:filename.find('_grid')]\n",
    "adffile = name + '/' + name + '/w001001.adf'\n",
    "\n",
    "if not os.path.exists(adffile):\n",
    "    r = requests.get(url, stream=True)\n",
    "    with open(filename, 'wb') as f:\n",
    "        total_length = int(r.headers.get('content-length'))\n",
    "        for chunk in tqdm(r.iter_content(chunk_size=1024), total=(total_length/1024) + 1):\n",
    "            if chunk:\n",
    "                f.write(chunk)\n",
    "                f.flush()\n",
    "    zip = zipfile.ZipFile(filename)\n",
    "    zip.extractall('.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We transform the data a bit so that rivers appear thicker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = rasterio.open(adffile)\n",
    "acc_orig = dataset.read()[0]\n",
    "acc = np.where(acc_orig<0, 0, acc_orig)\n",
    "\n",
    "shrink = 1 # if you are out of RAM try increasing this number (should be a power of 2)\n",
    "radius = 5 # you can play with this number to change the width of the rivers\n",
    "circle = np.zeros((2*radius+1, 2*radius+1)).astype('uint8')\n",
    "y, x = np.ogrid[-radius:radius+1,-radius:radius+1]\n",
    "index = x**2 + y**2 <= radius**2\n",
    "circle[index] = 1\n",
    "acc = np.sqrt(acc)\n",
    "acc = scipy.ndimage.maximum_filter(acc, footprint=circle)\n",
    "acc[acc_orig<0] = np.nan\n",
    "acc = np.array(acc[::shrink, ::shrink])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original data is in the WGS 84 projection, but Leaflet uses Web Mercator, so we need to reproject."
   ]
  },
  {
   "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",
    "with rasterio.Env():\n",
    "    rows, cols = acc.shape\n",
    "    src_transform = list(dataset.affine)\n",
    "    src_transform[0] *= shrink\n",
    "    src_transform[4] *= shrink\n",
    "    src_transform = Affine(*src_transform[:6])\n",
    "    src_crs = {'init': 'EPSG:4326'}\n",
    "    source = acc\n",
    "\n",
    "    dst_crs = {'init': 'EPSG:3857'}\n",
    "    dst_transform, width, height = rasterio.warp.calculate_default_transform(src_crs, dst_crs, cols, rows, *dataset.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",
    "acc_web = destination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's convert our NumPy array to an image. For that we must specify a colormap (here `plt.cm.jet`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc_norm = acc_web - np.nanmin(acc_web)\n",
    "acc_norm = acc_norm / np.nanmax(acc_norm)\n",
    "acc_norm = np.where(np.isfinite(acc_web), acc_norm, 0)\n",
    "acc_im = PIL.Image.fromarray(np.uint8(plt.cm.jet(acc_norm)*255))\n",
    "acc_mask = np.where(np.isfinite(acc_web), 255, 0)\n",
    "mask = PIL.Image.fromarray(np.uint8(acc_mask), mode='L')\n",
    "im = PIL.Image.new('RGBA', acc_norm.shape[::-1], color=None)\n",
    "im.paste(acc_im, mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image is embedded in the URL as a PNG file, so that it can be sent to the browser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if py3:\n",
    "    f = BytesIO()\n",
    "else:\n",
    "    f = StringIO()\n",
    "im.save(f, 'png')\n",
    "data = b64encode(f.getvalue())\n",
    "if py3:\n",
    "    data = data.decode('ascii')\n",
    "imgurl = 'data:image/png;base64,' + data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can overlay our image and if everything went fine it should be exactly over South America."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = dataset.bounds\n",
    "bounds = [(b.bottom, b.left), (b.top, b.right)]\n",
    "io = ImageOverlay(url=imgurl, bounds=bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = [-10, -60]\n",
    "zoom = 2\n",
    "m = Map(center=center, zoom=zoom)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tile = basemap_to_tiles(basemaps.Esri.WorldStreetMap)\n",
    "m.add_layer(tile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can play with the opacity slider and check that rivers from our data file match the rivers on OpenStreetMap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_layer(io)\n",
    "io.interact(opacity=(0.0,1.0,0.01))"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}