ENH: add step 4

display_api.py

@@ -17,7 +17,7 @@ from rasterio.warp import transform
 import warnings
 warnings.simplefilter(action='ignore', category=FutureWarning)
 
-from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl, LayersControl
+from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl, LayersControl, GeoJSON
 
 
 def normalize(array, quantile=2):
@@ -91,7 +91,7 @@ def write_quicklook(raster, filename, downfactor=4):
             dst.write(band, n+1)
 
 
-def rasters_on_map(rasters_list, out_dir, overlay_names_list):
+def rasters_on_map(rasters_list, out_dir, overlay_names_list, geojson_data=None):
     """
     displays a raster on a ipyleaflet map
     :param rasters_list: rasters to display (rasterio image)
@@ -121,6 +121,12 @@ def rasters_on_map(rasters_list, out_dir, overlay_names_list):
     m.add_control(LayersControl())
     m.add_control(FullScreenControl())
 
+    # - add geojson data
+    if geojson_data is not None:
+        geo_json = GeoJSON(data=geojson_data,
+                           style = {'color': 'green', 'opacity':1, 'weight':1.9, 'dashArray':'9', 'fillOpacity':0.1})
+        m.add_layer(geo_json)
+    
     # - add draw control
     dc = DrawControl()
     m.add_control(dc)

ogr_io.py

+"""
+Read/Write OGR functions
+"""
+from osgeo import ogr
+
+def openToRead(shapefile):
+    driver = ogr.GetDriverByName("SQLite")
+    if driver.Open(shapefile, 0):
+        dataSource = driver.Open(shapefile, 0)
+    else:
+        print("Not possible to open the file "+shapefile)
+        sys.exit(1)
+    return dataSource
+
+def openToReadWrite(shapefile):
+    driver = ogr.GetDriverByName("SQLite")
+    if driver.Open(shapefile, 1):
+        dataSource = driver.Open(shapefile, 1)
+    else:
+        print("Not possible to open the file "+shapefile)
+        sys.exit(1)
+    return dataSource

step4_filter_shapefile.ipynb

+{
+ "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": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Converting  output/ndwi_threshold30.tif  into a shapefile  output/watermask.sqlite\n"
+     ]
+    }
+   ],
+   "source": [
+    "import os\n",
+    "import subprocess\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",
+    "ndwi_thres_fname = os.path.join(OUTPUT_DIR, \"ndwi_threshold30.tif\")\n",
+    "watermask_fname = os.path.join(OUTPUT_DIR, \"watermask.sqlite\")\n",
+    "envelope_fname = os.path.join(OUTPUT_DIR, \"morbihan.sqlite\")\n",
+    "results_fname = os.path.join(OUTPUT_DIR, \"results.sqlite\")\n",
+    "geojson_fname = os.path.join(OUTPUT_DIR, \"geojson.json\")\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": 2,
+   "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": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Nb features (before filtering):  344\n",
+      "Nb islands (after filtering):  45\n"
+     ]
+    }
+   ],
+   "source": [
+    "from osgeo import ogr, osr\n",
+    "\n",
+    "# Save extent to a new Shapefile\n",
+    "outDriver = ogr.GetDriverByName(\"SQLite\")\n",
+    "\n",
+    "# Remove output shapefile if it already exists\n",
+    "if os.path.exists(results_fname):\n",
+    "    outDriver.DeleteDataSource(results_fname)\n",
+    "\n",
+    "# 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": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "17fd0284fe9f4b5c8e7a7861779afc87",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'attribution': 'Map data (c) <a href…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "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
+}