Detect water with several sensors

Let’s use EOReader for water detection, using several sensors over the same area.


# Imports
import os
import numpy as np
import xarray as xr

Core functions

  • A function loading the Normalized Water Index (NDWI)

  • A function extracting water through a threshold

def load_ndwi(prod, res):
    Load NDWI index (and rename the array)
    # Read NDWI index
    ndwi = prod.load(NDWI)[NDWI]
    ndwi_name = f"{ndwi.attrs['sensor']} NDWI"
    return ndwi.rename(ndwi_name)

def extract_water(ndwi):
    Extract water from NDWI index (and rename the array)
    # Assert water bodies when NDWI index > 0.2
    water = xr.where(ndwi > 0.2, 1, 0)

    # Set nodata where ndwi is nan.
    # WARNING: the function xr.DataArray.where sets by default np.nan where the condition is false !
    # See here:
    water = water.where(~np.isnan(ndwi))

    # Plot a subsampled version
    water_name = f"{ndwi.attrs['sensor']} WATER"
    water = water.rename(water_name)
    water.attrs["long_name"] = "Water detection"
    return water.rename(water_name)

Create logger

import logging

logger = logging.getLogger("eoreader")

# create console handler and set level to debug
ch = logging.StreamHandler()

# create formatter
formatter = logging.Formatter('%(message)s')

# add formatter to ch

# add ch to logger

Linking some data

First of all, we need some satellite data.
Let’s take 3 products covering approximately the same area (over DAX in France):

  • One Landsat-8 OLCI collection 2

  • One Landsat-5 TM collection 2

  • One Sentinel-2 L2A

prod_folder = os.path.join("/home", "data", "DS3", "CI", "extracteo", "water")
paths = [  
    # Landsat-8 OLCI collection 2    
    os.path.join(prod_folder, "LC08_L1TP_200030_20201220_20210310_02_T1.tar"),
    # Landsat-5 TM collection 2    
    os.path.join(prod_folder, "LT05_L1TP_200030_20111110_20200820_02_T1.tar"),
    # Sentinel-2 L2A
    os.path.join(prod_folder, ""),

Process these products

from eoreader.reader import Reader
from eoreader.bands import *

# Create the reader
eoreader = Reader()

# Loop on all the products
water_arrays = []
ndwi_arrays = []
extents = []
for path in paths:"*** {os.path.basename(path)} ***")
    # Open the product
    prod =, remove_tmp=True)
    # Get extents
    # Read NDWI index
    # 6Let's say we want a 60. meters resolution
    ndwi = load_ndwi(prod, res=60.)
    # Extract water
*** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
Loading bands ['NIR', 'GREEN']
Read NIR
Manage nodata for band NIR
Manage nodata for band GREEN
Loading index ['NDWI']

*** LT05_L1TP_200030_20111110_20200820_02_T1.tar ***
Loading bands ['NIR', 'GREEN']
Read NIR
Manage nodata for band NIR
Manage nodata for band GREEN
Loading index ['NDWI']

*** ***
Loading bands ['NIR', 'GREEN']
Read NIR
Manage nodata for band NIR
Manage nodata for band GREEN
Loading index ['NDWI']

Plot the result

# Plot the tiles
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import as ccrs
from shapely.errors import ShapelyDeprecationWarning
import warnings

warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning)

nrows = len(paths)
plt.figure(figsize=(6 * nrows, 6 * nrows))
cmap = ListedColormap(['khaki', 'lightblue'])
for i in range(nrows):
    # Compute cartopy projection (EOReader always is in UTM)
    extent = extents[i]
    str_epsg = str(
    zone = str_epsg[-2:]
    is_south = str_epsg[2] == 7
    proj = ccrs.UTM(zone, is_south)

    # Get extent values
    # The extents must be defined in the form (min_x, max_x, min_y, max_y)
    bounds = extent.bounds
    extent_val = [bounds.minx[0], bounds.maxx[0], bounds.miny[0], bounds.maxy[0]]

    # Plot NDWI
    axes = plt.subplot(nrows, 2, 2*i+1, projection=proj)
    axes.set_extent(extent_val, proj)
    ndwi_arrays[i][0, ::10, ::10].plot.imshow(origin='upper', extent=extent_val, transform=proj, robust=True)
    # Plot water
    axes = plt.subplot(nrows, 2, 2*i+2, projection=proj)
    axes.set_extent(extent_val, proj)
    water_arrays[i][0, ::10, ::10].plot.imshow(origin='upper', extent=extent_val, transform=proj, cmap=cmap, cbar_kwargs={'ticks': [0, 1]})