notebooks/LargeScaleMosaic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distributed processing using Dynamic Tiler + mosaicJSON\n",
    "\n",
    "\n",
    "### Goal\n",
    "Create Mosaic of Landsat 8, using rio-tiler-mosaic pixel selection methods\n",
    "\n",
    "\n",
    "### Requirement\n",
    "```\n",
    "!pip install tdqm cogeo-mosaic\n",
    "```\n",
    "\n",
    "![](https://user-images.githubusercontent.com/10407788/74463971-b47ae000-4e60-11ea-8413-b058e00a2b09.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import time\n",
    "import random\n",
    "import requests\n",
    "import datetime\n",
    "import urllib.parse\n",
    "from io import BytesIO\n",
    "from concurrent import futures\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import mercantile\n",
    "from supermercado.burntiles import tile_extrema\n",
    "from affine import Affine\n",
    "\n",
    "import rasterio\n",
    "from rasterio.io import MemoryFile\n",
    "from rasterio.windows import Window\n",
    "\n",
    "from rio_cogeo.utils import _meters_per_pixel\n",
    "from rio_cogeo.cogeo import cog_translate\n",
    "from rio_cogeo.profiles import cog_profiles\n",
    "\n",
    "from rio_tiler_mosaic.mosaic import _filter_futures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# See https://github.com/developmentseed/awspds-mosaic to deploy your own \n",
    "endpoint = \"{YOUR ENDPOINT HERE}\"\n",
    "\n",
    "# Define AOI here\n",
    "bounds = [-5.833740234375, 46.15700496290803, -0.9118652343749999, 49.23912083246698]\n",
    "\n",
    "# Date Filters\n",
    "start = datetime.datetime.strptime(\"2019-01-01\", \"%Y-%m-%d\").strftime(\"%Y-%m-%dT00:00:00Z\")\n",
    "end = datetime.datetime.strptime(\"2019-12-31\", \"%Y-%m-%d\").strftime(\"%Y-%m-%dT23:59:59Z\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SAT-API QUERY\n",
    "# see http://sat-utils.github.io/sat-api/\n",
    "query = {\n",
    "    \"bbox\": bounds,\n",
    "    \"time\": f\"{start}/{end}\",\n",
    "    \"query\": {\n",
    "        \"eo:sun_elevation\": {\"gt\": 0},\n",
    "        \"landsat:tier\": {\"eq\": \"T1\"},\n",
    "        \"collection\": {\"eq\": \"landsat-8-l1\"},\n",
    "        \"eo:cloud_cover\": {\"gte\": 0, \"lt\": 3}  # Cloud Filters\n",
    "    },\n",
    "    \"sort\": [{\n",
    "      \"field\": \"eo:cloud_cover\",\n",
    "      \"direction\": \"asc\"\n",
    "    }],\n",
    "}\n",
    "\n",
    "tilescale = 2  # Each tile will be 512x512px\n",
    "tilesize = 256 * tilescale\n",
    "\n",
    "# We post the query to the mosaic endpoint with some optional parameters\n",
    "results = requests.post(\n",
    "    f\"{endpoint}/mosaic/create\", \n",
    "    json=query,\n",
    "    params={\n",
    "        # Landsat covers zoom 7 to 12 but \n",
    "        # because we don't want to use the mosaic for visualization \n",
    "        # purpose we don't really care about lower zoom level. \n",
    "        # We could use zoom 11 but don't want to make the mosaicJSON file too big.\n",
    "        # Minzoom define the quadkey zoom and thus the number of quadkey list\n",
    "        # See https://github.com/developmentseed/mosaicjson-spec/tree/master/0.0.2\n",
    "        \"minzoom\": 9, \n",
    "        # We filter the season to have greenest data \n",
    "        \"seasons\": \"spring,summer\",\n",
    "        # Here we can also pass some tile option to be added to the tile url.\n",
    "        \"tile_format\":\"tif\",  # We use GeoTIFF output from the tiler\n",
    "        \"tile_scale\": tilescale,\n",
    "    },\n",
    ").json()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Mercator tiles covering the AOI\n",
    "# For Landsat the maxzoom is 12 (~38m resolution, in Web Mercator projection)\n",
    "# because we use 512x512 pixel tiles we can \n",
    "# fetch tiles at zoom 11 instead of 12.\n",
    "zoom = results[\"maxzoom\"] - (tilescale - 1)\n",
    "\n",
    "extrema = tile_extrema(bounds, zoom)\n",
    "tiles = []\n",
    "for x in range(extrema[\"x\"][\"min\"], extrema[\"x\"][\"max\"]):\n",
    "    for y in range(extrema[\"y\"][\"min\"], extrema[\"y\"][\"max\"]):\n",
    "        tiles.append(mercantile.Tile(z=zoom, x=x, y=y))\n",
    "\n",
    "random.shuffle(tiles)\n",
    "print(len(tiles))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create tile URL\n",
    "# See https://landsatlive.live/tiles/docs for the list of options you can pass to the tile url\n",
    "query_params = dict(\n",
    "    bands=\"4,3,2\", # True Color RGB\n",
    "    color_ops=\"gamma RGB 3.5, saturation 1.7, sigmoidal RGB 15 0.35\", # Looks nice\n",
    "    pixel_selection=\"median\",  # See https://github.com/cogeotiff/rio-tiler-mosaic/blob/master/rio_tiler_mosaic/methods/defaults.py#L86\n",
    ")\n",
    "tiles_url = results[\"tiles\"][0] + urllib.parse.urlencode(query_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Output COG parameters\n",
    "width = (extrema[\"x\"][\"max\"] - extrema[\"x\"][\"min\"]) * tilesize\n",
    "height = (extrema[\"y\"][\"max\"] - extrema[\"y\"][\"min\"]) * tilesize\n",
    "w, n = mercantile.xy(\n",
    "    *mercantile.ul(extrema[\"x\"][\"min\"], extrema[\"y\"][\"min\"], zoom)\n",
    ")\n",
    "res = _meters_per_pixel(zoom, 0, tilesize=tilesize)\n",
    "\n",
    "params = dict(\n",
    "    driver=\"GTiff\",\n",
    "    dtype=\"uint8\",\n",
    "    count=3,\n",
    "    width=width,\n",
    "    height=height,\n",
    "    crs=\"epsg:3857\",\n",
    "    transform=Affine(res, 0, w, 0, -res, n),\n",
    "    nodata=0,\n",
    "    tiled=True,\n",
    "    blockxsize=tilesize,\n",
    "    blockysize=tilesize,\n",
    ")\n",
    "output_profile = cog_profiles.get(\"deflate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _worker(tile, tile_url, extrema, tilesize=256, retry=0):\n",
    "    \"\"\"\n",
    "    Tile Worker.\n",
    "    \n",
    "    Call the tile endpoint for each mercator tile.\n",
    "    \n",
    "    \"\"\"\n",
    "    url = tile_url.format(z=tile.z, x=tile.x, y=tile.y)\n",
    "    img = requests.get(url)\n",
    "    if not img.status_code == 200:\n",
    "        time.sleep(3)\n",
    "        if retry == 3:\n",
    "            raise Exception(\"Empty\")\n",
    "        return _worker(tile, retry=retry+1)\n",
    "\n",
    "    row = (tile.y - extrema[\"y\"][\"min\"]) * tilesize\n",
    "    col = (tile.x - extrema[\"x\"][\"min\"]) * tilesize\n",
    "    window = Window(col_off=col, row_off=row, width=tilesize, height=tilesize)\n",
    "\n",
    "    return window, img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Distribute the processes using multithreading\n",
    "# and then re-assemble the tiles in one COG\n",
    "with rasterio.Env():\n",
    "    with MemoryFile() as memfile:\n",
    "        with memfile.open(**params) as mem:\n",
    "            with futures.ThreadPoolExecutor(max_workers=50) as executor:\n",
    "                future_work = [\n",
    "                    executor.submit(_worker, tile, tiles_url, extrema, tilesize) for tile in tiles\n",
    "                ]\n",
    "                \n",
    "                for f in tqdm(futures.as_completed(future_work), total=len(future_work)):               \n",
    "                    pass\n",
    "\n",
    "            for f in _filter_futures(future_work):\n",
    "                window, img = f\n",
    "                with rasterio.open(BytesIO(img.content)) as src_dst:\n",
    "                    mem.write(src_dst.read(indexes=(1,2,3)), window=window)\n",
    "\n",
    "            cog_translate(\n",
    "                mem,\n",
    "                f\"Bretagne_medianPixel_2019.tif\",\n",
    "                output_profile,\n",
    "                config=dict(GDAL_NUM_THREADS=\"ALL_CPUS\", GDAL_TIFF_OVR_BLOCKSIZE=\"128\"),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Capture d’écran, le 2020-02-13 à 12 51 36](https://user-images.githubusercontent.com/10407788/74463338-b2fce800-4e5f-11ea-9216-76ffb0dea2b4.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Using rio-tiler Expression and Custom Pixel_Selection methods\n",
    "\n",
    "### Create Mosaic and visualize True Color RGB with highest NDVI value\n",
    "\n",
    "The idea is to use rio-tiler expression method to create a 4th band on which we'll base the pixel selection method.\n",
    "\n",
    "The `lastband` method will:\n",
    "- analyse the stack of tiles\n",
    "- find where the value for the 4th band is the highest\n",
    "- construct a tile with the index obtained\n",
    "- automatically remove the last band of the tile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create tile URL\n",
    "# See https://landsatlive.live/tiles/docs for the list of options you can pass to the tile url\n",
    "query_params = dict(\n",
    "    expr=\"b4,b3,b2,(b5-b4)/(b5+b4)\",  # Red,Green,Blue,NDVI\n",
    "    color_ops=\"gamma RGB 3.5, saturation 1.7, sigmoidal RGB 15 0.35\", # Looks nice\n",
    "    pixel_selection=\"lastband\",  # See https://github.com/developmentseed/awspds-mosaic/blob/master/awspds_mosaic/pixel_methods.py#L30\n",
    ")\n",
    "tiles_url = results[\"tiles\"][0] + urllib.parse.urlencode(query_params)\n",
    "\n",
    "\n",
    "with rasterio.Env():\n",
    "    with MemoryFile() as memfile:\n",
    "        with memfile.open(**params) as mem:\n",
    "            with futures.ThreadPoolExecutor(max_workers=50) as executor:\n",
    "                future_work = [\n",
    "                    executor.submit(_worker, tile, tiles_url, extrema, tilesize) for tile in tiles\n",
    "                ]\n",
    "                \n",
    "                for f in tqdm(futures.as_completed(future_work), total=len(future_work)):               \n",
    "                    pass\n",
    "\n",
    "            for f in _filter_futures(future_work):\n",
    "                window, img = f\n",
    "                with rasterio.open(BytesIO(img.content)) as src_dst:\n",
    "                    mem.write(src_dst.read(indexes=(1,2,3)), window=window)\n",
    "\n",
    "            cog_translate(\n",
    "                mem,\n",
    "                f\"Bretagne_highNDVI_2019.tif\",\n",
    "                output_profile,\n",
    "                config=dict(GDAL_NUM_THREADS=\"ALL_CPUS\", GDAL_TIFF_OVR_BLOCKSIZE=\"128\"),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://user-images.githubusercontent.com/10407788/74463335-b2fce800-4e5f-11ea-80e2-9ac94d9132f5.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## At Scale\n",
    "\n",
    "Let's try at country scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# La France\n",
    "bounds = [-6.8115234375, 42.35854391749705, 9.580078125, 52.07950600379697]\n",
    "\n",
    "# Date Filters\n",
    "start = datetime.datetime.strptime(\"2013-01-01\", \"%Y-%m-%d\").strftime(\"%Y-%m-%dT00:00:00Z\")\n",
    "end = datetime.datetime.strptime(\"2019-12-31\", \"%Y-%m-%d\").strftime(\"%Y-%m-%dT23:59:59Z\")\n",
    "\n",
    "# SAT-API QUERY\n",
    "# see http://sat-utils.github.io/sat-api/\n",
    "query = {\n",
    "    \"bbox\": bounds,\n",
    "    \"time\": f\"{start}/{end}\",\n",
    "    \"query\": {\n",
    "        \"eo:sun_elevation\": {\"gt\": 0},\n",
    "        \"landsat:tier\": {\"eq\": \"T1\"},\n",
    "        \"collection\": {\"eq\": \"landsat-8-l1\"},\n",
    "        \"eo:cloud_cover\": {\"gte\": 0, \"lt\": 5}  # Cloud Filters\n",
    "    }\n",
    "}\n",
    "\n",
    "tilescale = 2  # Each tile will be 512x512px\n",
    "tilesize = 256 * tilescale\n",
    "\n",
    "results = requests.post(\n",
    "    f\"{endpoint}/mosaic/create\", \n",
    "    json=query,\n",
    "    params={\n",
    "        \"minzoom\": 9, \n",
    "        \"seasons\": \"summer\",\n",
    "        \"tile_format\":\"tif\",\n",
    "        \"tile_scale\": tilescale,\n",
    "    },\n",
    ").json()\n",
    "\n",
    "mosaicid = results[\"name\"]\n",
    "\n",
    "# List of tiles\n",
    "zoom = results[\"maxzoom\"] - (tilescale - 1)\n",
    "extrema = tile_extrema(bounds, zoom)\n",
    "tiles = []\n",
    "for x in range(extrema[\"x\"][\"min\"], extrema[\"x\"][\"max\"]):\n",
    "    for y in range(extrema[\"y\"][\"min\"], extrema[\"y\"][\"max\"]):\n",
    "        tiles.append(mercantile.Tile(z=zoom, x=x, y=y))\n",
    "random.shuffle(tiles)\n",
    "\n",
    "# Define Output COG parameters\n",
    "width = (extrema[\"x\"][\"max\"] - extrema[\"x\"][\"min\"]) * tilesize\n",
    "height = (extrema[\"y\"][\"max\"] - extrema[\"y\"][\"min\"]) * tilesize\n",
    "w, n = mercantile.xy(\n",
    "    *mercantile.ul(extrema[\"x\"][\"min\"], extrema[\"y\"][\"min\"], zoom)\n",
    ")\n",
    "res = _meters_per_pixel(zoom, 0, tilesize=tilesize)\n",
    "\n",
    "params = dict(\n",
    "    driver=\"GTiff\",\n",
    "    dtype=\"uint8\",\n",
    "    count=3,\n",
    "    width=width,\n",
    "    height=height,\n",
    "    crs=\"epsg:3857\",\n",
    "    transform=Affine(res, 0, w, 0, -res, n),\n",
    "    nodata=0,\n",
    "    tiled=True,\n",
    "    blockxsize=tilesize,\n",
    "    blockysize=tilesize,\n",
    ")\n",
    "output_profile = cog_profiles.get(\"deflate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create tile URL\n",
    "# See https://landsatlive.live/tiles/docs for the list of options you can pass to the tile url\n",
    "query_params = dict(\n",
    "    bands=\"4,3,2\", # True Color RGB\n",
    "    color_ops=\"gamma RGB 3.5, saturation 1.7, sigmoidal RGB 15 0.35\", # Looks nice\n",
    "    pixel_selection=\"median\",\n",
    ")\n",
    "\n",
    "# For this large scale process we are going to use another endpoint which can run up to 30seconds per tile.\n",
    "tiles_url = f\"{endpoint}/batch/{mosaicid}/{{z}}/{{x}}/{{y}}@{tilescale}x.tif?\" + urllib.parse.urlencode(query_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distribute the processes using multithreading\n",
    "# and then re-assemble the tiles in one COG\n",
    "with rasterio.Env():\n",
    "    with MemoryFile() as memfile:\n",
    "        with memfile.open(**params) as mem:\n",
    "            with futures.ThreadPoolExecutor(max_workers=100) as executor:\n",
    "                future_work = [\n",
    "                    executor.submit(_worker, tile, tiles_url, extrema, tilesize) for tile in tiles\n",
    "                ]\n",
    "                \n",
    "                for f in tqdm(futures.as_completed(future_work), total=len(future_work)):               \n",
    "                    pass\n",
    "\n",
    "            for f in _filter_futures(future_work):\n",
    "                window, img = f\n",
    "                with rasterio.open(BytesIO(img.content)) as src_dst:\n",
    "                    mem.write(src_dst.read(indexes=(1,2,3)), window=window)\n",
    "\n",
    "            cog_translate(\n",
    "                mem,\n",
    "                f\"France_Median_2013-2019_Z12.tif\",\n",
    "                output_profile,\n",
    "                config=dict(GDAL_NUM_THREADS=\"ALL_CPUS\", GDAL_TIFF_OVR_BLOCKSIZE=\"128\"),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://user-images.githubusercontent.com/10407788/74540814-70deaf80-4f0e-11ea-879f-f51954d545d0.jpg)"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}