Delineate water bodies (lakes) from Landsat 8 using machine learning with Python and QGIS - Tutorial

The process of retrieving information from satellite imagery could be time consuming and faces challenges like image resolution, feature recognition and user criteria. With the use of machine learning tools we can preprocess images and combine them with standard QGIS tools to delineate objects in a much efficient way.

We have done a tutorial for the delimitation of water bodies as lakes from a panchromatic Landsat 8 image. The tutorial shows a mixed procedure that detects edges with Python and Scikit Image , traces paths with QGIS and the Trace Raster plugin and finally gets the lake extension as a Shapefile.

Input files are available under this Github repository:

https://github.com/SaulMontoya/delimitateWaterBodiesMachineLearningPythonandQGIS.git

and also from this link.

To run the hakuchik image from Docker just type this on the Powershell:

docker run -it --name hakuchik -p 8888:8888 -p 6543:5432 hatarilabs/hakuchik:a2d1cb82a605

Tutorial

Python code

This is the Python code that detects edges on the Landsat 8 image:

Import required packages

#python and matplotlib libraries
import numpy as np
import matplotlib.pyplot as plt

#geospatial libraries
import rasterio

#machine learning libraries
#!pip install scikit-image
from skimage import feature

Import raster and read band

imgRaster = rasterio.open('../Rst/landsatImage.tif')
im = imgRaster.read(1)
fig = plt.figure(figsize=(12,8))
plt.imshow(im)
plt.show()

Canny filter

# Compute the Canny filter for two values of sigma
edges1 = feature.canny(im)
edges2 = feature.canny(im, sigma=3)

# display results
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(18, 6),
                                    sharex=True, sharey=True)

ax1.imshow(im, cmap=plt.cm.gray)
ax1.axis('off')
ax1.set_title('noisy image', fontsize=20)

ax2.imshow(edges1, cmap=plt.cm.gray)
ax2.axis('off')
ax2.set_title(r'Canny filter, $\sigma=1$', fontsize=20)

ax3.imshow(edges2, cmap=plt.cm.gray)
ax3.axis('off')
ax3.set_title(r'Canny filter, $\sigma=3$', fontsize=20)

fig.tight_layout()

plt.show()

Export edges as geospatial TIFF

rasterOut = rasterio.open(
    '../Rst/landsatImage_edges2.tif',
    'w',
    driver='GTiff',
    height=edges2.shape[0],
    width=edges2.shape[1],
    count=3,
    dtype=im.dtype,
    crs=imgRaster.crs,
    transform=imgRaster.transform)

rasterOut.write(np.uint8(edges2*60),1)
rasterOut.write(np.uint8(edges2*120),2)
rasterOut.write(np.uint8(edges2*180),3)
rasterOut.close()