{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left; margin:0px 15px 15px 0px; width:120px\" src=\"https://www.orfeo-toolbox.org/wp-content/uploads/2016/03/logo-orfeo-toolbox.png\">\n",
    "\n",
    "# OTB Guided Tour - JURSE 2019 Vannes France - May 21th\n",
    "## Yannick TANGUY and David YOUSSEFI (CNES, French Space Agency)\n",
    "\n",
    "<br>\n",
    "\n",
    "<b> Press <span style=\"color:black;background:yellow\">SHIFT+ENTER</span> to execute the notebook interactively cell by cell </b></div>\n",
    "\n",
    "## Step 4 : Filter shapefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polygonize with GDAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "from osgeo import ogr, osr\n",
    "\n",
    "import ogr_io\n",
    "\n",
    "# Data directory\n",
    "DATA_DIR = \"data\"\n",
    "\n",
    "# Output directory\n",
    "OUTPUT_DIR = \"output\"\n",
    "\n",
    "# Input / Output filenames\n",
    "envelope_fname = os.path.join(DATA_DIR, \"morbihan.sqlite\")\n",
    "ndwi_thres_fname = os.path.join(OUTPUT_DIR, \"ndwi_threshold30.tif\")\n",
    "watermask_fname = os.path.join(OUTPUT_DIR, \"watermask.sqlite\")\n",
    "results_fname = os.path.join(OUTPUT_DIR, \"results.sqlite\")\n",
    "geojson_fname = os.path.join(OUTPUT_DIR, \"geojson.json\")\n",
    "\n",
    "# Get SQLite driver\n",
    "outDriver = ogr.GetDriverByName(\"SQLite\")\n",
    "\n",
    "# Remove results if it already exists\n",
    "if os.path.exists(results_fname):\n",
    "    outDriver.DeleteDataSource(results_fname)\n",
    "    \n",
    "# Remove watermask if it already exists\n",
    "if os.path.exists(watermask_fname):\n",
    "    outDriver.DeleteDataSource(watermask_fname)\n",
    "\n",
    "# Convert the TIF file in a shapefile with polygons\n",
    "try:\n",
    "    subprocess.call([\"gdal_polygonize.py\", ndwi_thres_fname,\"-f\", \"SQLite\", watermask_fname])\n",
    "    print (\"Converting \",ndwi_thres_fname, \" into a shapefile \", watermask_fname)\n",
    "except:\n",
    "    print(\"Error during conversion\")\n",
    "    exit(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read shapefile with ogr\n",
    "watermask = ogr_io.openToRead(watermask_fname)\n",
    "envelope = ogr_io.openToRead(envelope_fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Crop Islands\n",
    "\n",
    "Crop the \"watermask\" shapefile with envelope shape (~ Morbihan gulf) and also filter the smallest areas and the two biggest (main land and sea)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the spatial reference, WGS84\n",
    "srs = osr.SpatialReference()\n",
    "srs.ImportFromEPSG(32630)\n",
    "\n",
    "# Create the output shapefile\n",
    "outDataSource = outDriver.CreateDataSource(results_fname)\n",
    "outLayer = outDataSource.CreateLayer(\"islands\", srs, geom_type=ogr.wkbPolygon)\n",
    "\n",
    "# Collect all islands inside the envelope\n",
    "geomcol = ogr.Geometry(ogr.wkbGeometryCollection)\n",
    "inLayer = watermask.GetLayer()\n",
    "layer_envelope = envelope.GetLayer()\n",
    "\n",
    "for shape in layer_envelope:\n",
    "    inLayer.SetSpatialFilter(shape.GetGeometryRef())\n",
    "    print(\"Nb features (before filtering): \",inLayer.GetFeatureCount())\n",
    "    for feature in inLayer:\n",
    "        outLayer.CreateFeature(feature)\n",
    "\n",
    "featureDefn = outLayer.GetLayerDefn()\n",
    "\n",
    "areaMax = 0\n",
    "index = 1\n",
    "biggestFeature = 0\n",
    "for feature in outLayer:\n",
    "    area = feature.GetGeometryRef().GetArea()\n",
    "    # Our features contain : \n",
    "    # - some land (because the cropping shape is not very precise\n",
    "    # - the sea !\n",
    "    # - some very small submarine dark areas : these are supposely rocks\n",
    "    # --> we eliminate every feature larger than 4 km² (\"ile aux moines\" area) \n",
    "    # and also smaller than 10000 m² (1 Ha)\n",
    "    if (area > 4000000.0):\n",
    "        #print(\"Del too large feature area = \",str(area))\n",
    "        outLayer.DeleteFeature(feature.GetFID())\n",
    "    if (area < 10000.0):\n",
    "        #print(\"Del too small feature area = \",str(area))\n",
    "        outLayer.DeleteFeature(feature.GetFID())\n",
    "\n",
    "print (\"Nb islands (after filtering): \", outLayer.GetFeatureCount())\n",
    "    \n",
    "# unload the driver -> writes the file\n",
    "outDriver = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from glob import glob\n",
    "import display_api\n",
    "import json\n",
    "\n",
    "# Convert sqlite to geojson EPSG:4326\n",
    "try:\n",
    "    subprocess.call([\"ogr2ogr\", \"-f\", \"GeoJSON\", geojson_fname, results_fname, \"-t_srs\", \"EPSG:4326\"])\n",
    "    print (\"Converting sqlite to geojson EPSG:4326\")\n",
    "except:\n",
    "    print(\"Error during conversion\")\n",
    "    exit(-1)\n",
    "    \n",
    "with open(geojson_fname, 'r') as f:\n",
    "    fc = json.load(f)\n",
    "\n",
    "DATE = \"20180805\"\n",
    "raster = rasterio.open(glob(os.path.join(DATA_DIR, \"SENTINEL2?_{0}-*.tif\".format(DATE)))[0])\n",
    "m, dc = display_api.rasters_on_map([raster], OUTPUT_DIR, [DATE], geojson_data=fc)\n",
    "m"
   ]
  },
  {
   "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.5.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}