step4_filter_shapefile.ipynb 6.14 KB
Newer Older
David Youssefi's avatar
David Youssefi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
{
 "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",
David Youssefi's avatar
David Youssefi committed
28
   "execution_count": null,
David Youssefi's avatar
David Youssefi committed
29
   "metadata": {},
David Youssefi's avatar
David Youssefi committed
30
   "outputs": [],
David Youssefi's avatar
David Youssefi committed
31 32 33
   "source": [
    "import os\n",
    "import subprocess\n",
David Youssefi's avatar
David Youssefi committed
34
    "from osgeo import ogr, osr\n",
David Youssefi's avatar
David Youssefi committed
35 36 37 38 39 40 41 42 43 44
    "\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",
David Youssefi's avatar
David Youssefi committed
45
    "envelope_fname = os.path.join(DATA_DIR, \"morbihan.sqlite\")\n",
David Youssefi's avatar
David Youssefi committed
46 47 48 49 50
    "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",
David Youssefi's avatar
David Youssefi committed
51 52 53 54 55 56 57 58 59 60 61
    "# 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",
David Youssefi's avatar
David Youssefi committed
62 63 64 65 66 67 68 69 70 71 72
    "# 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",
David Youssefi's avatar
David Youssefi committed
73
   "execution_count": null,
David Youssefi's avatar
David Youssefi committed
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
   "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",
David Youssefi's avatar
David Youssefi committed
93
   "execution_count": null,
David Youssefi's avatar
David Youssefi committed
94
   "metadata": {},
David Youssefi's avatar
David Youssefi committed
95
   "outputs": [],
David Youssefi's avatar
David Youssefi committed
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
   "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",
David Youssefi's avatar
David Youssefi committed
144
   "execution_count": null,
David Youssefi's avatar
David Youssefi committed
145
   "metadata": {},
David Youssefi's avatar
David Youssefi committed
146
   "outputs": [],
David Youssefi's avatar
David Youssefi committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
   "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"
   ]
David Youssefi's avatar
David Youssefi committed
169 170 171 172 173 174 175
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
David Youssefi's avatar
David Youssefi committed
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
  }
 ],
 "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
}