examples/introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First install the following packages:\n",
    "\n",
    "- `requests`\n",
    "- `tqdm`\n",
    "- `xarray-leaflet`\n",
    "- `dask`\n",
    "\n",
    "The recommended way is: `conda install -c conda-forge requests tqdm xarray_leaflet dask`.\n",
    "\n",
    "We start by importing the modules that we will need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import os\n",
    "import warnings\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import zipfile\n",
    "import rioxarray\n",
    "import xarray_leaflet\n",
    "from ipyleaflet import Map, basemaps, LayersControl, WidgetControl\n",
    "from ipywidgets import FloatSlider\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will display the DEM (digital elevation model) for Asia from the [HydroSHEDS](https://hydrosheds.org) dataset, which represents the terrain. Let's first download the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/hydrosheds/sa_30s_zip_grid/as_dem_30s_grid.zip'\n",
    "filename = os.path.basename(url)\n",
    "name = filename[:filename.find('_grid')]\n",
    "adffile = os.path.join(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": [
    "It is a dataset that [Rasterio](https://rasterio.readthedocs.io) can open, but in order to get a nice `DataArray` with all the metadata taken into account, we open it with rioxarray:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = rioxarray.open_rasterio(adffile, masked=True)\n",
    "da"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The projection is `EPSG:4326` (aka `WGS84`). Here the coordinate `x` corresponds to longitudes, and `y` to latitudes (in degrees). There is only one band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = da.sel(band=1)\n",
    "da.name = 'DEM'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset can be too big to hold in memory, so we will chunk it into smaller pieces. That will also improve performances as the generation of a tile can be done in parallel using Dask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da = da.chunk((1000, 1000))\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just need to create a map before passing it to our `DataArray` extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Map(center=[40, 115], zoom=3, basemap=basemaps.CartoDB.DarkMatter, interpolation='nearest')\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To show our data on the map, we call `leaflet.plot()` on our `DataArray`, and pass the map as parameter. We get back a layer, that we can further control with e.g. a slider to set the opacity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "l = da.leaflet.plot(m, colormap=plt.cm.terrain)\n",
    "# l.interact(opacity=(0., 1.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have a slider outside of the map, but we are going to insert it in the map. We will also insert a layer control."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layers_control = LayersControl(position='topright')\n",
    "m.add_control(layers_control)\n",
    "\n",
    "opacity_slider = FloatSlider(description='Opacity:', min=0, max=1, value=1)\n",
    "\n",
    "def set_opacity(change):\n",
    "    l.opacity = change['new']\n",
    "\n",
    "opacity_slider.observe(set_opacity, names='value')\n",
    "slider_control = WidgetControl(widget=opacity_slider, position='bottomleft')\n",
    "m.add_control(slider_control)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `select()` method allows to select a region by clicking and dragging a box on the map (first click on the █ button)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da.leaflet.select()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then get the selected data array back, and e.g. plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_selected = da.leaflet.get_selection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if da_selected is not None:\n",
    "    da_selected.plot.imshow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da.leaflet.unselect()  # remove the selection buttons"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}