Measuring the water level of a Theewaterskloof Dam in South Africa¶

Natebook showcases an example of Earth observation processing chain that determines water levels of any water body (dam, reservoir, lake, ...) from satellite imagery. The entire processing chain is performed using the eo-learn package. The user simply needs to provide a polygon with water body's nominal water extent.

In [1]:

%reload_ext autoreload
%autoreload 2
%matplotlib inline

Imports¶

eo-learn imports¶

In [2]:

from eolearn.core import EOTask, EOPatch, LinearWorkflow, Dependency, FeatureType

# We'll use Sentinel-2 imagery (Level-1C) provided through Sentinel Hub
# If you don't know what `Level 1C` means, don't worry. It doesn't matter.
from eolearn.io import S2L1CWCSInput 
from eolearn.core import LoadFromDisk, SaveToDisk

# cloud detection
from eolearn.mask import AddCloudMaskTask, get_s2_pixel_cloud_detector
from eolearn.mask import AddValidDataMaskTask

# filtering of scenes
from eolearn.features import SimpleFilterTask

# burning the vectorised polygon to raster
from eolearn.geometry import VectorToRaster

/home/matej/.local/share/virtualenvs/eo-learn-Gx1GfgPz/lib/python3.7/site-packages/numba/types/containers.py:3: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working
  from collections import Iterable

Other imports¶

In [3]:

# The golden standard: numpy and matplotlib
import numpy as np

# import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

# For manipulating geo-spatial vector dataset (polygons of nominal water extent)
import geopandas as gpd

# Image manipulations
# Our water detector is going to be based on a simple threshold 
# of Normalised Difference Water Index (NDWI) grayscale image
from skimage.filters import threshold_otsu

# Loading polygon of nominal water extent
import shapely.wkt
from shapely.geometry import Polygon

# sentinelhub-py package
from sentinelhub import BBox, CRS

Water level extraction EOWorkflow¶

Our basic logic of the example workflow is:

Download all available Sentinel-2 sattelite imagery of Theewaterskloof Dam from beginning of 2016 and today
- the following layers will be downloaded:
  - TRUE_COLOR for nicer visualisations
  - NDWI for water detection
Clouds are very often obscuring the view of the ground. In order to correctly determine the water level of the dam all images with clouds need to be filtered out.
Apply adaptive thresholding to NDWI grayscale images
Extract water level from a comparison of measured water extent with the nominal one

Each step in the above overview of the workflow is accomplished by adding an EOTask to the EOWorkflow

Load the Polygon of nominal water extent and define a BBOX¶

The BBOX defines an area of interest and will be used to create an EOPatch.

In [4]:

# The polygon of the dam is written in wkt format and WGS84 coordinate reference system
with open('theewaterskloof_dam_nominal.wkt', 'r') as f:
    dam_wkt = f.read()

dam_nominal = shapely.wkt.loads(dam_wkt)

# inflate the BBOX 
inflate_bbox = 0.1
minx, miny, maxx, maxy = dam_nominal.bounds

delx = maxx - minx
dely = maxy - miny
minx = minx - delx * inflate_bbox
maxx = maxx + delx * inflate_bbox
miny = miny - dely * inflate_bbox
maxy = maxy + dely * inflate_bbox
    
dam_bbox = BBox([minx, miny, maxx, maxy], crs=CRS.WGS84)

dam_bbox.geometry - dam_nominal

Out[4]:

Step 1: Intialize (and implement workflow specific) EOTasks¶

Create an EOPatch and add all EO features (satellite imagery data)¶

In [5]:

input_task = S2L1CWCSInput('TRUE-COLOR-S2-L1C', resx='20m', resy='20m', maxcc=0.5, instance_id=None)

add_ndwi = S2L1CWCSInput('NDWI')

Burn in the nominal water extent¶

The VectorToRaster task expects the vectorised dataset in geopandas dataframe.

In [6]:

dam_gdf = gpd.GeoDataFrame(crs=CRS.WGS84.pyproj_crs(), geometry=[dam_nominal])

In [7]:

dam_gdf.plot();

In [8]:

dam_gdf.plot();

In [9]:

add_nominal_water = VectorToRaster(dam_gdf, (FeatureType.MASK_TIMELESS, 'NOMINAL_WATER'), values=1, 
                                   raster_shape=(FeatureType.MASK, 'IS_DATA'), raster_dtype=np.uint8)

Run s2cloudless cloud detector and filter out scenes with cloud coverage >20%¶

To speed up the process the cloud detection is executed at lower resolution (160m). The resulting cloud probability map and binary mask are stored as CLP and CLM features in EOPatch.

In [10]:

cloud_classifier = get_s2_pixel_cloud_detector(average_over=2, dilation_size=1, all_bands=False)

cloud_detection = AddCloudMaskTask(cloud_classifier, 'BANDS-S2CLOUDLESS', cm_size_y='160m', cm_size_x='160m', 
                                   cmask_feature='CLM', cprobs_feature='CLP', instance_id=None)

Define a VALID_DATA layer: pixel has to contain data and should be classified as clear sky by the cloud detector (CLM equals 0)

In [11]:

def calculate_valid_data_mask(eopatch):
    return np.logical_and(eopatch.mask['IS_DATA'].astype(np.bool), 
                          np.logical_not(eopatch.mask['CLM'].astype(np.bool)))

add_valid_mask = AddValidDataMaskTask(predicate=calculate_valid_data_mask)

Calculate fraction of valid pixels per frame and store it as SCALAR feature

In [12]:

def calculate_coverage(array):
    return 1.0 - np.count_nonzero(array) / np.size(array)

class AddValidDataCoverage(EOTask):
    
    def execute(self, eopatch):
        
        valid_data = eopatch.get_feature(FeatureType.MASK, 'VALID_DATA')
        time, height, width, channels = valid_data.shape
        
        coverage = np.apply_along_axis(calculate_coverage, 1, valid_data.reshape((time, height * width * channels)))
        
        eopatch.add_feature(FeatureType.SCALAR, 'COVERAGE', coverage[:, np.newaxis])
        return eopatch
    
add_coverage = AddValidDataCoverage()

Filter out too cloudy scenes

In [13]:

class ValidDataCoveragePredicate:
    
    def __init__(self, threshold):
        self.threshold = threshold
        
    def __call__(self, array):
        return calculate_coverage(array) < self.threshold
    
remove_cloudy_scenes = SimpleFilterTask((FeatureType.MASK, 'VALID_DATA'), ValidDataCoveragePredicate(0.2))

Apply Water Detection¶

In [14]:

class WaterDetector(EOTask):
    
    @staticmethod
    def detect_water(ndwi):
        """
        Very simple water detector based on Otsu thresholding method of NDWI.
        """
        otsu_thr = 1.0
        if len(np.unique(ndwi)) > 1:
            otsu_thr = threshold_otsu(ndwi)

        return ndwi > otsu_thr

    def execute(self, eopatch):
        water_masks = np.asarray([self.detect_water(ndwi[...,0]) for ndwi in eopatch.data['NDWI']])
        
        # we're only interested in the water within the dam borders
        water_masks = water_masks[...,np.newaxis] * eopatch.mask_timeless['NOMINAL_WATER']
        
        water_levels = np.asarray([np.count_nonzero(mask)/np.count_nonzero(eopatch.mask_timeless['NOMINAL_WATER']) 
                                   for mask in water_masks])
        
        eopatch.add_feature(FeatureType.MASK, 'WATER_MASK', water_masks)
        eopatch.add_feature(FeatureType.SCALAR, 'WATER_LEVEL', water_levels[...,np.newaxis])
        
        return eopatch
    
water_detection = WaterDetector()

Step 2: Define the EOWorkflow¶

In [15]:

workflow = LinearWorkflow(input_task, add_ndwi, cloud_detection, add_nominal_water, add_valid_mask, add_coverage,
                          remove_cloudy_scenes, water_detection)

Step 3: Run the workflow¶

Process all Sentinel-2 acquisitions from beginning of 2016 and until end of August 2018.

In [16]:

time_interval = ['2016-01-01','2018-08-31']

In [17]:

result = workflow.execute({
    input_task: {
        'bbox': dam_bbox,
        'time_interval': time_interval
    },
})

In [18]:

patch = list(result.values())[-1]

Print content of eopatch at the end of the workflow execution

In [19]:

patch

Out[19]:

EOPatch(
  data: {
    CLP: numpy.ndarray(shape=(64, 618, 928, 1), dtype=float32)
    NDWI: numpy.ndarray(shape=(64, 618, 928, 1), dtype=float32)
    TRUE-COLOR-S2-L1C: numpy.ndarray(shape=(64, 618, 928, 3), dtype=float32)
  }
  mask: {
    CLM: numpy.ndarray(shape=(64, 618, 928, 1), dtype=bool)
    IS_DATA: numpy.ndarray(shape=(64, 618, 928, 1), dtype=bool)
    VALID_DATA: numpy.ndarray(shape=(64, 618, 928, 1), dtype=bool)
    WATER_MASK: numpy.ndarray(shape=(64, 618, 928, 1), dtype=uint8)
  }
  scalar: {
    COVERAGE: numpy.ndarray(shape=(64, 1), dtype=float64)
    WATER_LEVEL: numpy.ndarray(shape=(64, 1), dtype=float64)
  }
  label: {}
  vector: {}
  data_timeless: {}
  mask_timeless: {
    NOMINAL_WATER: numpy.ndarray(shape=(618, 928, 1), dtype=uint8)
  }
  scalar_timeless: {}
  label_timeless: {}
  vector_timeless: {}
  meta_info: {
    maxcc: 0.5
    service_type: 'wcs'
    size_x: '20m'
    size_y: '20m'
    time_difference: datetime.timedelta(days=-1, seconds=86399)
    time_interval: ['2016-01-01', '2018-08-31']
  }
  bbox: BBox(((19.10818927, -34.08851246), (19.30962163, -33.977424140000004)), crs=CRS('4326'))
  timestamp: [datetime.datetime(2016, 4, 6, 8, 38, 18), ..., datetime.datetime(2018, 8, 24, 8, 40, 30)], length=64
)

Plot results¶

In [20]:

from skimage.filters import sobel
from skimage.morphology import disk
from skimage.morphology import erosion, dilation, opening, closing, white_tophat

In [21]:

def plot_rgb_w_water(eopatch, idx):
    ratio = np.abs(eopatch.bbox.max_x - eopatch.bbox.min_x) / np.abs(eopatch.bbox.max_y - eopatch.bbox.min_y)
    fig, ax = plt.subplots(figsize=(ratio * 10, 10))
    
    ax.imshow(eopatch.data['TRUE-COLOR-S2-L1C'][idx])
    
    observed = closing(eopatch.mask['WATER_MASK'][idx,...,0], disk(1))
    nominal = sobel(eopatch.mask_timeless['NOMINAL_WATER'][...,0])
    observed = sobel(observed)
    nominal = np.ma.masked_where(nominal == False, nominal)
    observed = np.ma.masked_where(observed == False, observed)
    
    ax.imshow(nominal, cmap=plt.cm.Reds)
    ax.imshow(observed, cmap=plt.cm.Blues)
    ax.axis('off')

In [22]:

plot_rgb_w_water(patch, 0)

In [23]:

plot_rgb_w_water(patch, -1)

In [24]:

def plot_water_levels(eopatch, max_coverage=1.0):
    fig, ax = plt.subplots(figsize=(20, 7))

    dates = np.asarray(eopatch.timestamp)
    ax.plot(dates[eopatch.scalar['COVERAGE'][...,0] < max_coverage],
            eopatch.scalar['WATER_LEVEL'][eopatch.scalar['COVERAGE'][...,0] < max_coverage],
            'bo-', alpha=0.7)
    ax.plot(dates[eopatch.scalar['COVERAGE'][...,0] < max_coverage],
            eopatch.scalar['COVERAGE'][eopatch.scalar['COVERAGE'][...,0] < max_coverage],
            '--', color='gray', alpha=0.7)
    ax.set_ylim(0.0, 1.1)
    ax.set_xlabel('Date')
    ax.set_ylabel('Water level')
    ax.set_title('Theewaterskloof Dam Water Levels')
    ax.grid(axis='y')
    return ax

In [25]:

ax = plot_water_levels(patch, 1.0)

One observation at the end of 2017 sticks out. Let's visualise it.

In [26]:

idx = np.where(patch.scalar['WATER_LEVEL'] > 0.8)[0][0]

In [27]:

plot_rgb_w_water(patch, idx)

As can be seen the image from this date is covered with clouds, which interfer with water detection algorithms. Let's filter the cloudy scene out and plot only water levels for observations with cloud coverage below 2%.

In [28]:

ax = plot_water_levels(patch, 0.02)