Detect water with several constellations

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

/home/docs/checkouts/readthedocs.org/user_builds/eoreader/envs/stable/lib/python3.9/site-packages/dask/dataframe/__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

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-05-05 11:29:42,977 - [INFO] - *** LC08_L1TP_200030_20201220_20210310_02_T1.tar ***
2025-05-05 11:29:42,981 - [WARNING] - There is no existing products in EOReader corresponding to /home/prods/LANDSATS_COL2/LC08_L1TP_200030_20201220_20210310_02_T1.tar.
2025-05-05 11:29:42,982 - [INFO] - Your given path may not be a satellite image. If it is, maybe the product isn't handled by EOReader. If you are sure this product is handled, it is either corrupted or you may need to go deeper in the filetree to find the correct path to give.
2025-05-05 11:29:42,982 - [DEBUG] - Please look at what folder you should give to EOReader by accessing the documentation: https://eoreader.readthedocs.io/latest/main_features.html#recognized-paths
2025-05-05 11:29:42,983 - [INFO] - None

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 16
     13 logger.info(prod)
     15 # Get extents
---> 16 extents.append(prod.extent())
     18 # Read NDWI index
     19 # Let's say we want a 60. meters pixel_size
     20 ndwi = load_ndwi(prod, pixel_size=60.)

AttributeError: 'NoneType' object has no attribute 'extent'

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)
../_images/8a4c377c1d40de2f4e61d6657cce1bd0175ca20dad52de8f88e19b1d8618973e.png
Contents