Skip to content

GitLab

O
OTB Guided Tour
Commit 34282070 authored by David Youssefi's avatar David Youssefi
Browse files
Options

ENH: add step1 + display_api.py and utils.py

parent 368858e3
Hide whitespace changes
Inline Side-by-side
Showing with 517 additions and 102 deletions
README.md
View file @ 34282070
......@@ -9,7 +9,7 @@
4. Create Virtual Machine (2048Mo, VMDK, 15Gio) and install VirtualBox Guest Additions
5. Install librairies and softwares
- Miniconda : ./Miniconda3-latest-Linux-x86_64.sh
- Create otb env : conda create -n otb python=3.5 numpy rasterio jupyter matplotlib
- Create otb env : conda create -n otb python=3.5 numpy rasterio jupyter matplotlib gdal
- Install ipyleaflet : pip install --user ipyleaflet +
- OTB : ./OTB-6.6.1-Linux64.run --target /home/otb/OTB
- Jupyter :
......
display_api.py 0 → 100644
View file @ 34282070
"""
Some functions for Jupyter displaying
- quicklook creation
- map displaying (imshow, ipyleaflet)
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.warp import transform_bounds
from rasterio.warp import transform
# ignore rasterio FurtureWarnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl, LayersControl
def normalize(array, quantile=2):
"""
normalizes bands into 0-255 scale
:param array: numpy array to normalize
:param quantile: quantile for ignoring outliers
"""
array = np.nan_to_num(array)
array_min, array_max = np.percentile(array, quantile), np.percentile(array, 100-quantile)
normalized = 255*((array - array_min)/(array_max - array_min))
normalized[normalized>255] = 255
normalized[normalized<0] = 0
return normalized.astype(np.uint8)
def imshow_RGBPIR(raster, colors):
"""
show R,G,B, PIR image
:param array: raster (rasterio image)
:param colors: colors names list
"""
print ("Quicklook by channel")
multiband = []
fig, ax = plt.subplots(1,len(colors)+1, figsize=(21,7))
for nfig, color in enumerate(colors):
band = raster.read(nfig+1, out_shape=(int(raster.height/16), int(raster.width/16)))
band = normalize(band)
if nfig < 3: multiband.append(band)
try:
ax[nfig].imshow(band, cmap=color+'s')
except ValueError:
ax[nfig].imshow(band, cmap='Oranges')
ax[nfig].set_title(color+' channel')
ax[len(colors)].imshow(np.dstack(multiband))
ax[len(colors)].set_title(colors[0]+colors[1]+colors[2])
def write_quicklook(raster, filename, downfactor=4):
"""
write a JPG preview
:param raster: raster (rasterio image)
:param filename: path to the output preview
:param downfactor: downsampling factor
"""
profile = raster.profile
# update size in profile
newwidth = int(raster.width/downfactor)
newheight = int(raster.height/downfactor)
aff = raster.affine
newaffine = rasterio.Affine(aff.a/downfactor, aff.b, aff.c,
aff.d, aff.e/downfactor, aff.f)
profile.update(dtype=rasterio.uint8, count=3, compress='lzw', driver='JPEG',
width=newwidth, height=newheight, transform=newaffine, affine=newaffine)
# write raster
with rasterio.open(filename, 'w', **profile) as dst:
for n in range(3):
if raster.count == 1:
band = raster.read(1, out_shape=(int(raster.height/downfactor), int(raster.width/downfactor)))
else:
band = raster.read(n+1, out_shape=(int(raster.height/downfactor), int(raster.width/downfactor)))
band = normalize(band)
dst.write(band, n+1)
def rasters_on_map(rasters_list, out_dir, overlay_names_list):
"""
displays a raster on a ipyleaflet map
:param rasters_list: rasters to display (rasterio image)
:param out_dir: path to the output directory (preview writing)
:param overlay_names_list: name of the overlays for the map
"""
# - get bounding box
raster = rasters_list[0]
epsg4326 = {'init': 'EPSG:4326'}
bounds = transform_bounds(raster.crs, epsg4326, *raster.bounds)
center = [(bounds[0]+bounds[2])/2, (bounds[1]+bounds[3])/2]
# - get centered map
m = Map(center=(center[-1], center[0]), zoom=10)
# - plot quicklook
for raster, overlay_name in zip(rasters_list, overlay_names_list):
bounds = transform_bounds(raster.crs, epsg4326, *raster.bounds)
quicklook_url = os.path.join(out_dir, os.path.splitext(overlay_name.replace("/", "_"))[0] + "_PREVIEW.JPG")
write_quicklook(raster, quicklook_url)
quicklook = ImageOverlay(
url=quicklook_url,
bounds=((bounds[1], bounds[0]),(bounds[3], bounds[2])),
name=overlay_name
)
m.add_layer(quicklook)
m.add_control(LayersControl())
m.add_control(FullScreenControl())
# - add draw control
dc = DrawControl()
m.add_control(dc)
return m, dc
def get_bounding_box_from_draw(raster, dc):
"""
returns bounding box in a image
:param raster: displayed raster (rasterio image)
:param dc: drawcontrol ipyleaflet
"""
try:
# Get last draw from the map
coordinates = np.array(dc.last_draw['geometry']['coordinates'][0])
# lon, lat to Sentinel-2 coordinate reference system
lon_min_max = [np.amin(coordinates[:,0]), np.amax(coordinates[:,0])]
lat_min_max = [np.amin(coordinates[:,1]), np.amax(coordinates[:,1])]
lons, lats = np.meshgrid(lon_min_max, lat_min_max)
xs, ys = transform({'init': 'EPSG:4326'}, raster.crs, lons.flatten(), lats.flatten())
# Get the region of interest in the image
rows, cols = [], []
for x, y in zip(xs, ys):
row, col = ~raster.affine * (x, y)
rows.append(row)
cols.append(col)
row_min, row_max = np.amin(rows), np.amax(rows)
col_min, col_max = np.amin(cols), np.amax(cols)
startx, starty = map(lambda v:int(np.floor(v)), [col_min, row_min])
endx, endy = map(lambda v:int(np.ceil(v)), [col_max, row_max])
print ("Bounding box computed")
except:
print ("Draw a polygon in the displayed map")
startx, starty, endx, endy = None, None, None, None
return startx, starty, endx, endy
otb_workshop.ipynb
View file @ 34282070
......@@ -6,8 +6,8 @@
"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",
"# OTB Guided Tour - JURSE 2019 Vannes France - May 21th\n",
"## Yannick TANGUY and David YOUSSEFI (CNES, French Space Agency)\n",
"\n",
"<br>\n",
"\n",
......@@ -16,7 +16,25 @@
"This tutorial will present the ORFEO Toolbox and showcase available applications for processing and manipulating satellite imagery.\n",
"\n",
"As a first step, we will learn how to manipulate OTB-Applications through their different interfaces.\n",
"Then, we will follow a guided tour of two main processing frameworks (segmentation & classification): we will use data from different sensors (Pléiades or Spot 6/7 with a Sentinel 2 time series) to make a value-added classification map of the region of Rennes (Brittany).\n"
"Then, we will follow a guided tour of two main processing frameworks (segmentation & classification): we will use data from different sensors (Pléiades or Spot 6/7 with a Sentinel 2 time series) to make a value-added classification map of the region of Rennes (Brittany).\n",
"\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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sentinel 2 Dataset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Information about dataset: directory (./data), band order, mgrs, level, dates"
]
},
{
......@@ -27,50 +45,66 @@
"source": [
"from glob import glob\n",
"import rasterio\n",
"from rasterio.plot import show\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import os\n",
"import display_api\n",
"import utils\n",
"\n",
"# Data directory\n",
"DATA_DIR = \"data\"\n",
"\n",
"# Product information\n",
"MGRS = \"30TWT\"\n",
"LEVEL = \"L2A\"\n",
"DATE = \"20180805\"\n",
"# Output directory\n",
"OUTPUT_DIR = \"output\"\n",
"\n",
"# Color bands\n",
"COLORS_BANDS = ['Blue', 'Green', 'Red', 'NIR']\n",
"\n",
"# Get S2 product\n",
"S2_product = glob(os.path.join(DATA_DIR, \"SENTINEL2?_{0}-*_{1}_T{2}_*\".format(DATE, LEVEL, MGRS)))[0]\n",
"COLOR_BANDS = ['Red', 'Green', 'Blue', 'NIR']\n",
"\n",
"# Normalize bands into 0.0 - 1.0 scale\n",
"def normalize(array, quantile=2):\n",
" array_min, array_max = np.percentile(array, quantile), np.percentile(array, 100-quantile)\n",
" normalized = ((array - array_min)/(array_max - array_min))\n",
" normalized[normalized>1] = 1\n",
" normalized[normalized<0] = 0\n",
" return normalized\n",
"# Product information\n",
"MGRS = \"30TWT\"\n",
"LEVEL = \"L2A\"\n",
"\n",
"# Date list\n",
"print (list(map(lambda x:os.path.basename(x)[11:11+8],glob(os.path.join(DATA_DIR, \"SENTINEL2*.tif\")))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Choose your dataset (by date)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"DATE = \"20180805\"\n",
"raster = rasterio.open(glob(os.path.join(DATA_DIR, \"SENTINEL2?_{0}-*_{1}_T{2}_*.tif\".format(DATE, LEVEL, MGRS)))[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Display with imshow"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot bands from SENTINEL-2 product\n",
"raster = rasterio.open(S2_product)\n",
"quicklooks = {}\n",
"\n",
"fig, ax = plt.subplots(1,len(COLORS_BANDS)+1, figsize=(21,7))\n",
"for nfig, color in enumerate(COLORS_BANDS):\n",
" quicklooks[color] = normalize(raster.read(nfig+1, out_shape=(int(raster.height / 8), int(raster.width / 8))))\n",
" try:\n",
" ax[nfig].imshow(quicklooks[color], cmap=color+'s')\n",
" except ValueError:\n",
" ax[nfig].imshow(quicklooks[color], cmap='Oranges')\n",
" ax[nfig].set_title(color+' channel')\n",
"\n",
"ax[len(COLORS_BANDS)].imshow(np.dstack((quicklooks['Red'], quicklooks['Green'], quicklooks['Blue'])))\n",
"ax[len(COLORS_BANDS)].set_title('RGB')\n",
"\n",
"print (\"Quicklook by channel\")"
"display_api.imshow_RGBPIR(raster, colors=COLOR_BANDS)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Display on a map"
]
},
{
......@@ -81,41 +115,12 @@
},
"outputs": [],
"source": [
"from rasterio.warp import transform_bounds\n",
"from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl, LayersControl\n",
"\n",
"# Plot footprint on a map\n",
"\n",
"# - get bounding box\n",
"epsg4326 = {'init': 'EPSG:4326'}\n",
"bounds = transform_bounds(raster.crs, epsg4326, *raster.bounds)\n",
"center = [(bounds[0]+bounds[2])/2, (bounds[1]+bounds[3])/2]\n",
"\n",
"# - get centered map\n",
"m = Map(center=(center[-1], center[0]), zoom=10)\n",
"\n",
"# - create quicklook\n",
"import otbApplication\n",
"quicklook_url = os.path.splitext(S2_product)[0] + \"_PREVIEW.JPG\"\n",
"app = otbApplication.Registry.CreateApplication(\"DynamicConvert\")\n",
"app.SetParameterString('in', S2_product)\n",
"app.SetParameterString('out', quicklook_url)\n",
"app.ExecuteAndWriteOutput()\n",
"\n",
"# - plot quicklook\n",
"quicklook = ImageOverlay(\n",
" url=quicklook_url,\n",
" bounds=((bounds[1], bounds[0]),(bounds[3], bounds[2])),\n",
" name=\"QUICKLOOK S2A\"\n",
")\n",
"\n",
"m.add_layer(quicklook)\n",
"m.add_control(LayersControl())\n",
"m.add_control(FullScreenControl())\n",
"\n",
"dc = DrawControl()\n",
"m.add_control(dc)\n",
"m\n"
"# displays a raster on a ipyleaflet map (\n",
"# rasters_list: rasters to display (rasterio image),\n",
"# out_dir: path to the output directory (preview writing)\n",
"# overlay_names_list: name of the overlays for the map)\n",
"m, dc = display_api.rasters_on_map([raster], OUTPUT_DIR, [DATE])\n",
"m"
]
},
{
......@@ -124,40 +129,16 @@
"metadata": {},
"outputs": [],
"source": [
"from rasterio.warp import transform\n",
"\n",
"def get_bounding_box_from_draw(dc=dc):\n",
" # Get last draw from the map\n",
" coordinates = np.array(dc.last_draw['geometry']['coordinates'][0])\n",
" \n",
" # lon, lat to Sentinel-2 coordinate reference system\n",
" lon_min_max = [np.amin(coordinates[:,0]), np.amax(coordinates[:,0])]\n",
" lat_min_max = [np.amin(coordinates[:,1]), np.amax(coordinates[:,1])]\n",
" lons, lats = np.meshgrid(lon_min_max, lat_min_max)\n",
" xs, ys = transform({'init': 'EPSG:4326'}, raster.crs, lons.flatten(), lats.flatten())\n",
" \n",
" # Get the region of interest in the image\n",
" rows, cols = [], []\n",
" for x, y in zip(xs, ys):\n",
" row, col = ~raster.affine * (x, y)\n",
" rows.append(row)\n",
" cols.append(col)\n",
" row_min, row_max = np.amin(rows), np.amax(rows)\n",
" col_min, col_max = np.amin(cols), np.amax(cols)\n",
" startx, starty = map(lambda v:int(np.floor(v)), [col_min, row_min])\n",
" endx, endy = map(lambda v:int(np.ceil(v)), [col_max, row_max]) \n",
" return startx, starty, endx, endy\n",
"\n",
"startx, starty, endx, endy = get_bounding_box_from_draw()\n",
"startx, starty, endx, endy = display_api.get_bounding_box_from_draw(raster, dc)\n",
"print (startx, starty, endx, endy)"
]
},
{
"cell_type": "code",
"execution_count": null,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": []
"source": [
"### Go to the [Step 1](http://localhost:8888/notebooks/otb-guided-tour/step1_compute_radiometric_indices.ipynb)"
]
}
],
"metadata": {
......
step1_compute_radiometric_indices.ipynb 0 → 100644
View file @ 34282070
{
"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 1 : Compute radiometric indices\n",
"\n",
"The aim of this first exercise is to compute radiometric indices from Sentinel2 images\n",
"- The NDVI is already computed\n",
"- you have to compute the NDWI2 and save it in a new file, suffixed by _NDWI2.tif\n",
"\n",
"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}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Given example: compute ndvi"
]
},
{
"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",
"print(utils.list_images(DATA_DIR,\"_RVBPIR.tif\"))\n",
"\n",
"def compute_ndvi(images, out_dir):\n",
" \"\"\"\n",
" Compute NDVI from images with Red, Green, Blue and NIR bands\n",
" :param images: list of images. Their filenames match the pattern [...]_RVBPIR.tif\n",
" :param out_dir: path to the output directory\n",
" \"\"\"\n",
" out_list = []\n",
" for im in images:\n",
" out = os.path.join(out_dir, os.path.basename(im.replace(\"_RVBPIR.tif\",\"_NDVI.tif\")))\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\", out)\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\",\"(im1b4-im1b1)/(im1b4+im1b1)\")\n",
" app.ExecuteAndWriteOutput()\n",
" out_list.append(out)\n",
" \n",
" return out_list"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute radiometric indices (NDVI & NDWI2) on a stack of images\n",
"ndvi_filenames_list = compute_ndvi(utils.list_images(DATA_DIR,\"_RVBPIR.tif\"), OUTPUT_DIR)\n",
"print (ndvi_filenames_list)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Display on a map\n",
"rasters = list(map(rasterio.open, ndvi_filenames_list))\n",
"overnames = list(map(os.path.basename, ndvi_filenames_list))\n",
"m, dc = display_api.rasters_on_map(rasters, OUTPUT_DIR, overnames)\n",
"m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute ndwi : <b> Fill the <span style=\"color:black;background:yellow\">compute_ndwi</span> function </b> "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def compute_ndwi(images, out_dir):\n",
" \"\"\"\n",
" Compute NDWI2 from images with Red, Green, Blue and NIR bands\n",
" :param images: list of images. Their filenames match the pattern [...]_RVBPIR.tif\n",
" :param out_dir: path to the output directory\n",
" \"\"\"\n",
" # TODO: fill the compute_ndwi_function\n",
" return None"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute radiometric indices (NDVI & NDWI2) on a stack of images\n",
"ndwi_filenames_list = compute_ndwi(utils.list_images(DATA_DIR,\"_RVBPIR.tif\"), OUTPUT_DIR)\n",
"if ndwi_filenames_list is not None:\n",
" print (ndwi_filenames_list)\n",
"else:\n",
" print (\"TODO: fill the compute_ndwi_function\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if ndwi_filenames_list is not None:\n",
" # Display on a map\n",
" rasters = list(map(rasterio.open, ndwi_filenames_list))\n",
" overnames = list(map(os.path.basename, ndwi_filenames_list))\n",
" m, dc = display_api.rasters_on_map(rasters, OUTPUT_DIR, overnames)\n",
" m\n",
"else:\n",
" print (\"TODO: fill the compute_ndwi_function\")"
]
}
],
"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
}
step2_compute_watermask.ipynb 0 → 100644
View file @ 34282070
{
"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 : \n",
"\n",
"The aim of this second exercise is to compute a water mask, from NDWI2 values\n",
"- what kind of function do we need to implement ?\n",
"- are all the pixels valid ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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\",\"_ndwi2.tif\"))\n",
" app.SetParameterString(\"exp\",\"(im1b2-im1b4)/(im1b2+im1b4)\")\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"
]
}
],
"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
}
utils.py 0 → 100644
View file @ 34282070
"""
Some high-level file operations
"""
import os
def list_images(data_dir, suffix):
"""
returns a list of file names, contained in a directory, with a name that matches a suffix