Commit 3b122953 authored by David Youssefi's avatar David Youssefi

ENH: Remove old ipynb + NDWI threshold

parent f7c244cd
......@@ -111,11 +111,15 @@
"source": [
"### Compute ndwi : <b> Fill the <span style=\"color:black;background:yellow\">compute_ndwi</span> function </b> \n",
"\n",
"**Reminder** : NDWI2 is computed from green and nir bands (defined by McFeeters, 1996):\n",
"**Reminder:** NDWI2 is computed from green and nir bands (defined by McFeeters, 1996):\n",
"\n",
"\\begin{equation}\n",
"\\mbox{NDWI2}=\\frac{(Xgreen - Xnir)}{(Xgreen + Xnir)}\n",
"\\end{equation}"
"\\end{equation}\n",
"\n",
"**NB:** For this second variant of the NDWI, a threshold can also be found in https://www.mdpi.com/2072-4292/5/7/3544/htm (McFeeters, 2013):\n",
"- ```< 0.3 - Non-water```\n",
"- ```>= 0.3 - Water```"
]
},
{
......
{
"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",
"\n",
"## Step 2 : Compute Watermask\n",
"\n",
"The aim of this second exercise is to combine NDWI2 values to create a water mask\n",
"- what kind of function do we need to implement ?\n",
"- are all the pixels valid ?\n",
"- we want to maximize water extent\n",
"\n",
"We are also going to use OTB pipeline, to avoid writing intermediate file : \n",
"- the first step will be plugged as input of the second step\n",
"- only the final application will launch the pipeline, with ExecuteAndWriteOutput()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import otbApplication\n",
"import rasterio\n",
"\n",
"import display_api\n",
"import utils\n",
"\n",
"# Data directory\n",
"DATA_DIR = \"data\"\n",
"\n",
"# Output directory\n",
"OUTPUT_DIR = \"output\"\n",
"\n",
"# NDWI synthesis\n",
"NDWI_synthesis = os.path.join(OUTPUT_DIR,\"synthesis_NDWI2.tif\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute the NDWI2\n",
"\n",
"Do we need to adapt the formula to deal with invalid pixels ? \n",
"\n",
"What happens with clouds ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def app_radiometricindices(images, out_dir):\n",
" index = 0\n",
" for im in images:\n",
" # Here we create an application\n",
" app = otbApplication.Registry.CreateApplication(\"BandMath\")\n",
" # BandMath takes a list of images as input,\n",
" # so we have to give a Python list [ image ], or [ image1, image2, .., imageN]\n",
" app.SetParameterStringList(\"il\",[im])\n",
" out = os.path.join(out_dir, os.path.basename(im.replace(\"_RVBPIR.tif\",\"_NDWI2.tif\")))\n",
" app.SetParameterString(\"out\",out)\n",
" app.SetParameterString(\"exp\",\"(im1b2-im1b4)/(im1b2+im1b4)\")\n",
" # Plug the application in the pipeline \n",
" # does NOT launch the pipeline yet !\n",
" app.Execute()\n",
" # we save the application in a list\n",
" if (index == 0):\n",
" tab_App = [app]\n",
" else:\n",
" tab_App.append(app)\n",
" index += 1\n",
" # the list of applications is returned\n",
" return tab_App"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute watermask from a list of NDWI2 images\n",
"\n",
"We are going to use a list of images to compute a watermask (0 for land, 1 for water) : \n",
"- the inputs must be set in a list with app.AddImageToParameterInputImageList\n",
"- the expression will be use all the images im1b1, im2b1, etc."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def compute_watermask(list_app,outputfilename):\n",
" app = otbApplication.Registry.CreateApplication(\"BandMath\")\n",
"\n",
" index = 1\n",
" input_list = []\n",
" exp = \"?????\"\n",
" for appli in list_app:\n",
" app.AddImageToParameterInputImageList(\"il\",appli.GetParameterOutputImage(\"out\"))\n",
" if index==1:\n",
" exp = exp # + \"some function that involves im1b1\"\n",
" else:\n",
" exp = exp # +\", \"+\"im\"+str(index)+\"b1\"def compute_watermask(list_app,outputfilename):\n",
" index += 1\n",
"\n",
" print(exp)\n",
"\n",
" app.SetParameterString(\"out\",outputfilename)\n",
" app.SetParameterString(\"exp\",exp)\n",
" app.ExecuteAndWriteOutput()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute radiometric indices (NDVI & NDWI2) on a stack of images\n",
"list_app = app_radiometricindices(utils.list_images(DATA_DIR,\"_RVBPIR.tif\"), OUTPUT_DIR)\n",
"compute_watermask(list_app, NDWI_synthesis)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import rasterio\n",
"import display_api\n",
"\n",
"raster = rasterio.open(NDWI_synthesis)\n",
"m, dc = display_api.rasters_on_map([raster], OUTPUT_DIR, [\"NDWI synthesis\"])\n",
"m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Go to the [Step 3](./step3_segment_watermask_threshold_NDWI.ipynb)"
]
}
],
"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
}
......@@ -93,7 +93,7 @@
"source": [
"### Threshold ndwi : <b> Fill the <span style=\"color:black;background:yellow\">threshold_ndwi</span> function </b> \n",
"\n",
"**Tips:** The formula to threshold the NDWI images can be written using\n",
"**Tips:** The formula to threshold (th=0.3) the NDWI images can be written using\n",
"- binary operators:\n",
" - ‘+’ addition, ‘-‘ subtraction, ‘*’ multiplication, ‘/’ division\n",
" - ‘^’ raise x to the power of y\n",
......
{
"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 3 : Segment watermask thresholding NDWI"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Functions of the full OTB Pipeline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import otbApplication\n",
"\n",
"def app_radiometricindices(images):\n",
" index = 0\n",
" for im in images:\n",
" # Here we create an application\n",
" app = otbApplication.Registry.CreateApplication(\"BandMath\")\n",
" # BandMath takes a list of images as input,\n",
" # so we have to give a Python list [ image ], or [ image1, image2, .., imageN]\n",
" app.SetParameterStringList(\"il\",[im])\n",
" app.SetParameterString(\"out\",im.replace(\"RVBPIR.tif\",\"_indices.tif\"))\n",
" # the main parameter is the mathematical expression\n",
" # here, we compute NDVI : (nir - red) / (nir + red)\n",
" # The corresponding bands for NIR and Red are respectively the 4th and the 3rd bands (b4, b3)\n",
" # of the first image (im1)\n",
" app.SetParameterString(\"exp\",\"(im1b2+im1b4) > 0 ? (im1b2-im1b4)/(im1b2+im1b4) : -1 \")\n",
" app.Execute()\n",
" if (index == 0):\n",
" tab_App = [app]\n",
" else:\n",
" tab_App.append(app)\n",
" index += 1\n",
" return tab_App\n",
"\n",
"def compute_watermask(list_app,outputfilename):\n",
" app = otbApplication.Registry.CreateApplication(\"BandMath\")\n",
"\n",
" index = 1\n",
" input_list = []\n",
" exp = \"max(\"\n",
" for appli in list_app:\n",
" app.AddImageToParameterInputImageList(\"il\",appli.GetParameterOutputImage(\"out\"))\n",
" if index==1:\n",
" exp = exp+\"im\"+str(index)+\"b1\"\n",
" else:\n",
" exp = exp+\", \"+\"im\"+str(index)+\"b1\"\n",
" index += 1\n",
"\n",
" exp = exp +\")\"\n",
" print(exp)\n",
" print(input_list)\n",
"\n",
" app.SetParameterString(\"out\",outputfilename)\n",
" app.SetParameterString(\"exp\",exp)\n",
" app.Execute()\n",
" return app\n",
"\n",
"def thresholdNDWI(appWatermask,outputfilename,threshold):\n",
" app = otbApplication.Registry.CreateApplication(\"BandMath\")\n",
" app.AddImageToParameterInputImageList(\"il\", appWatermask.GetParameterOutputImage(\"out\"))\n",
" app.SetParameterString(\"exp\", \"im1b1 > \"+str(threshold)+\" ? 1 : 0\")\n",
" app.SetParameterString(\"out\", outputfilename)\n",
" app.ExecuteAndWriteOutput()\n",
" return app\n",
"\n",
"def smoothResult(thresholdNDWI, radius, outputfilename):\n",
" app = otbApplication.Registry.CreateApplication(\"BinaryMorphologicalOperation\")\n",
" app.SetParameterInputImage(\"in\",thresholdNDWI.GetParameterOutputImage(\"out\"))\n",
" app.SetParameterString(\"out\",outputfilename)\n",
" app.SetParameterString(\"structype\",\"ball\")\n",
" app.SetParameterFloat(\"structype.ball.xradius\", radius)\n",
" app.SetParameterFloat(\"structype.ball.yradius\", radius)\n",
" app.SetParameterString(\"filter\",\"closing\")\n",
" app.ExecuteAndWriteOutput()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Full OTB Pipeline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"\n",
"import utils\n",
"\n",
"# Threshold\n",
"threshold = 0.3\n",
"\n",
"# Data directory\n",
"DATA_DIR = \"data\"\n",
"\n",
"# Output directory\n",
"OUTPUT_DIR = \"output\"\n",
"\n",
"# Input / Output filenames\n",
"images_list = utils.list_images(DATA_DIR,\"_RVBPIR.tif\")\n",
"ndwi_max_fname = os.path.join(OUTPUT_DIR, \"max_NDWI2.tif\")\n",
"ndwi_thres_fname = os.path.join(OUTPUT_DIR, \"ndwi_threshold%d.tif\" % (threshold*100) )\n",
"\n",
"# Compute radiometric indices on a stack of images\n",
"list_app = app_radiometricindices(images_list)\n",
"\n",
"# Compute a water mask from NDWI2\n",
"wmask = compute_watermask(list_app, ndwi_max_fname)\n",
"\n",
"# Use a threshold to extract the mask\n",
"ndwi = thresholdNDWI(wmask, ndwi_thres_fname,threshold)\n",
"\n",
"# Smooth the result by applying a morphological operation\n",
"# smoothResult(ndwi,2.0,wdir+\"/result_smoothed.tif\")\n",
"\n",
"print (\"Full OTB Pipeline completed\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Display OTB result"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import rasterio\n",
"import display_api\n",
"\n",
"raster = rasterio.open(ndwi_thres_fname)\n",
"m, dc = display_api.rasters_on_map([raster], OUTPUT_DIR, [\"OTB result\"])\n",
"m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Go to the [Step 4](step4_filter_shapefile.ipynb)"
]
}
],
"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
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment