Detect water with several constellations#
Let’s use EOReader for water detection, using several constellations over the same area.
Imports#
# Imports
import os
import numpy as np
import xarray as xr
from eoreader.reader import Reader
from eoreader.bands import NDWI
Core functions#
A function loading the Normalized Water Index (
NDWI)A function extracting water through a threshold
def load_ndwi(prod, pixel_size):
"""
Load NDWI index (and rename the array)
"""
# Read NDWI index
ndwi = prod.load(NDWI, pixel_size=pixel_size)[NDWI]
ndwi_name = f"{ndwi.attrs['constellation']} 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: http://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html
water = water.where(~np.isnan(ndwi))
# Plot a subsampled version
water_name = f"{ndwi.attrs['constellation']} WATER"
water = water.rename(water_name)
water.attrs["long_name"] = "Water detection"
return water.rename(water_name)
Create logger#
# Create logger
import logging
from sertit import logs
logger = logging.getLogger("eoreader")
logs.init_logger(logger)
Linking some data#
Let’s take 3 products covering approximately the same area (over DAX in France):
One Landsat-8 OLI-TIRS collection 2
One Landsat-5 TM collection 2
One Sentinel-2 L1C
prod_folder = os.path.join("/home", "prods")
paths = [
# Landsat-8 OLI-TIRS collection 2
os.path.join(prod_folder, "LANDSATS_COL2", "LC08_L1TP_200030_20201220_20210310_02_T1.tar"),
# Landsat-5 TM collection 2
os.path.join(prod_folder, "LANDSATS_COL2", "LT05_L1TP_200030_20111110_20200820_02_T1.tar"),
# Sentinel-2 L2A
os.path.join(prod_folder, "S2", "PB 02.07+", "S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE"),
]
Process these products#
# Create the reader
reader = Reader()
# Loop on all the products
water_arrays = []
ndwi_arrays = []
extents = []
for path in paths:
logger.info(f"*** {os.path.basename(path)} ***")
# Open the product
prod = reader.open(path, remove_tmp=True)
logger.info(prod)
# Get extents
extents.append(prod.extent())
# Read NDWI index
# Let's say we want a 60. meters pixel_size
ndwi = load_ndwi(prod, pixel_size=60.)
ndwi_arrays.append(ndwi)
# Extract water
water_arrays.append(extract_water(ndwi))
logger.info("\n")
2025-12-23 12:04:17,953 - [INFO] - *** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
2025-12-23 12:04:18,511 - [INFO] - eoreader.LandsatProduct 'LC08_L1TP_200030_20201220_20210310_02_T1'
Attributes:
condensed_name: 20201220T104856_L8_200030_OLI_TIRS
path: /home/prods/LANDSATS_COL2/LC08_L1TP_200030_20201220_20210310_02_T1.tar
constellation: Landsat-8
sensor type: Optical
product type: L1
default pixel size: 30.0
default resolution: 30.0
acquisition datetime: 2020-12-20T10:48:56
band mapping:
COASTAL_AEROSOL: 1
BLUE: 2
GREEN: 3
RED: 4
NIR: 5
NARROW_NIR: 5
SWIR_CIRRUS: 9
SWIR_1: 6
SWIR_2: 7
THERMAL_IR_1: 10
THERMAL_IR_2: 11
PANCHROMATIC: 8
needs extraction: False
cloud cover: 16.36
tile name: 200030
2025-12-23 12:04:18,708 - [DEBUG] - Loading bands ['GREEN', 'NIR']
2025-12-23 12:04:18,710 - [DEBUG] - Read GREEN
2025-12-23 12:04:18,776 - [DEBUG] - Manage nodata for band GREEN
2025-12-23 12:04:18,826 - [DEBUG] - Converting GREEN to reflectance (if needed)
2025-12-23 12:04:18,828 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
/opt/conda/lib/python3.11/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
dest = _reproject(
2025-12-23 12:04:26,301 - [DEBUG] - Read NIR
2025-12-23 12:04:26,374 - [DEBUG] - Manage nodata for band NIR
2025-12-23 12:04:26,421 - [DEBUG] - Converting NIR to reflectance (if needed)
2025-12-23 12:04:26,423 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
2025-12-23 12:04:33,880 - [DEBUG] - Loading indices ['NDWI']
2025-12-23 12:04:43,970 - [INFO] -
2025-12-23 12:04:43,970 - [INFO] - *** LT05_L1TP_200030_20111110_20200820_02_T1.tar ***
2025-12-23 12:04:44,005 - [INFO] - eoreader.LandsatProduct 'LT05_L1TP_200030_20111110_20200820_02_T1'
Attributes:
condensed_name: 20111110T103612_L5_200030_TM
path: /home/prods/LANDSATS_COL2/LT05_L1TP_200030_20111110_20200820_02_T1.tar
constellation: Landsat-5
sensor type: Optical
product type: L1
default pixel size: 30.0
default resolution: 30.0
acquisition datetime: 2011-11-10T10:36:12
band mapping:
BLUE: 1
GREEN: 2
RED: 3
NIR: 4
NARROW_NIR: 4
SWIR_1: 5
SWIR_2: 7
THERMAL_IR_1: 6
THERMAL_IR_2: 6
needs extraction: False
cloud cover: 19.0
tile name: 200030
2025-12-23 12:04:44,080 - [DEBUG] - Loading bands ['GREEN', 'NIR']
2025-12-23 12:04:44,082 - [DEBUG] - Read GREEN
2025-12-23 12:04:44,131 - [DEBUG] - Manage nodata for band GREEN
2025-12-23 12:04:44,180 - [DEBUG] - Converting GREEN to reflectance (if needed)
2025-12-23 12:04:44,182 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
2025-12-23 12:04:50,575 - [DEBUG] - Read NIR
2025-12-23 12:04:50,665 - [DEBUG] - Manage nodata for band NIR
2025-12-23 12:04:50,726 - [DEBUG] - Converting NIR to reflectance (if needed)
2025-12-23 12:04:50,737 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
2025-12-23 12:04:57,618 - [DEBUG] - Loading indices ['NDWI']
2025-12-23 12:05:07,384 - [INFO] -
2025-12-23 12:05:07,384 - [INFO] - *** S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE ***
2025-12-23 12:05:07,460 - [INFO] - eoreader.S2Product 'S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155'
Attributes:
condensed_name: 20191215T110441_S2_T30TXP_L1C_114155
path: /home/prods/S2/PB 02.07+/S2A_MSIL1C_20191215T110441_N0208_R094_T30TXP_20191215T114155.SAFE
constellation: Sentinel-2
sensor type: Optical
product type: MSIL1C
default pixel size: 10.0
default resolution: 10.0
acquisition datetime: 2019-12-15T11:04:41
band mapping:
COASTAL_AEROSOL: 01
BLUE: 02
GREEN: 03
RED: 04
VEGETATION_RED_EDGE_1: 05
VEGETATION_RED_EDGE_2: 06
VEGETATION_RED_EDGE_3: 07
NIR: 08
NARROW_NIR: 8A
WATER_VAPOUR: 09
SWIR_CIRRUS: 10
SWIR_1: 11
SWIR_2: 12
needs extraction: False
cloud cover: 0.6668
tile name: T30TXP
2025-12-23 12:05:07,537 - [DEBUG] - Loading bands ['GREEN', 'NIR']
2025-12-23 12:05:07,559 - [DEBUG] - Read GREEN
2025-12-23 12:05:08,220 - [DEBUG] - Manage nodata for band GREEN
2025-12-23 12:05:08,482 - [DEBUG] - Converting GREEN to reflectance (if needed)
2025-12-23 12:05:08,523 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
2025-12-23 12:05:19,652 - [DEBUG] - Read NIR
2025-12-23 12:05:19,745 - [DEBUG] - Manage nodata for band NIR
2025-12-23 12:05:19,979 - [DEBUG] - Converting NIR to reflectance (if needed)
2025-12-23 12:05:19,981 - [DEBUG] - Clip the reflectance array to 0 as minimum value (in some cases, reflectance can have higher value than 1)
2025-12-23 12:05:31,340 - [DEBUG] - Loading indices ['NDWI']
2025-12-23 12:05:40,737 - [INFO] -
Plot the result#
# Plot the tiles
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import cartopy.crs as ccrs
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(extent.crs.to_epsg())
proj = ccrs.UTM(
zone=str_epsg[-2:],
southern_hemisphere=str_epsg[2] == 7
)
# 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, cmap="GnBu", robust=True)
axes.coastlines(linewidth=1)
axes.set_title(ndwi_arrays[i].name)
# 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]})
axes.coastlines(linewidth=1)
axes.set_title(water_arrays[i].name)
