Detecting water with radar

Sign up to the DEA Sandbox to run this notebook interactively from a browser
Compatibility: Notebook currently compatible with the DEA Sandbox environment only
Products used: s1_gamma0_geotif_scene

Background

Over 40% of the world’s population lives within 100 km of the coastline. However, coastal environments are constantly changing, with erosion and coastal change presenting a major challenge to valuable coastal infrastructure and important ecological habitats. Up-to-date data on the position of the coastline is essential for coastal managers to be able to identify and minimise the impacts of coastal change and erosion.

While coastlines can be mapped using optical data (demonstrated in the Coastal Erosion notebook), these images can be strongly affected by the weather, especially through the presence of clouds, which obscure the land and water below. This can be a particular problem in cloudy regions (e.g. southern Australia) or areas where wet season clouds prevent optical satellites from taking clear images for many months of the year.

Sentinel-1 use case

Radar observations are largely unaffected by cloud cover, so can take reliable measurements of areas in any weather. Radar data is readily available from the ESA/EC Copernicus program’s Sentinel-1 satellites. The two satellites provide all-weather observations, with a revisit time of 6 days. By developing a process to classify the observed pixels as either water or land, it is possible to identify the shoreline from radar data.

Description

In this example, we use data from the Sentinel-1 satellites to build a classifier that can determine whether a pixel is water or land in radar data. Specifically, this notebook uses an analysis-ready radar product known as backscatter, which describes the strength of the signal recieved by the satellite. The worked example takes users through the code required to:

Load Sentinel-1 backscatter data for an area of interest
Visualise the returned data.
Perform pre-processing steps on the Sentinel-1 bands.
Design a classifier to distinguish land and water.
Apply the classifier to the area of interest and interpret the results.
Investigate how to identify coastal change or the effect of tides.

Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

After finishing the analysis, return to the “Analysis parameters” cell, modify some values (e.g. choose a different location or time period to analyse) and re-run the analysis. There are additional instructions on modifying the notebook at the end.

Load packages

Load key Python packages and supporting functions for the analysis.

[1]:

%matplotlib inline

import datacube
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from scipy.ndimage import uniform_filter

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.plotting import display_map

Connect to the datacube

Activate the datacube database, which provides functionality for loading and displaying stored Earth observation data.

[2]:

dc = datacube.Datacube(app="Radar_water_detection")

Analysis parameters

The following cell sets the parameters, which define the area of interest and the length of time to conduct the analysis over. The parameters are

latitude: The latitude range to analyse (e.g. (-11.288, -11.086)). For reasonable loading times, make sure the range spans less than ~0.1 degrees.
longitude: The longitude range to analyse (e.g. (130.324, 130.453)). For reasonable loading times, make sure the range spans less than ~0.1 degrees.

If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results. The example covers Melville Island, which sits off the coast of the Northern Territory, Australia. The study area also contains an additional small island, which will be useful for assessing how well radar data distinguishes between land and water.

To run the notebook for a different area, make sure Sentinel-1 data is available for the chosen area using the DEA Sandbox Explorer.

[3]:

# Define the area of interest
latitude = (-11.288, -11.086)
longitude = (130.324, 130.453)

View the selected location

The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point.

[4]:

display_map(x=longitude, y=latitude)

[4]:

Make this Notebook Trusted to load map: File -> Trust Notebook

Preprocessing the data through filtering

Speckle Filtering using Lee Filter

You may have noticed that the water in the individual VV and VH images isn’t a consistent colour. The distortion you’re seeing is a type of noise known as speckle, which gives the images a grainy appearence. If we want to be able to easily decide whether any particular pixel is water or land, we need to reduce the chance of misinterpreting a water pixel as a land pixel due to the noise.

Speckle can be removed through filtering. If interested, you can find a technical introduction to speckle filtering here. For now, it is enough to know that we can filter the data using the Python function defined in the next cell:

[12]:

# Adapted from https://stackoverflow.com/questions/39785970/speckle-lee-filter-in-python
def lee_filter(img, size):
    """
    Applies the Lee filter to reduce speckle noise in an image.

    Parameters:
    img (ndarray): Input image to be filtered.
    size (int): Size of the uniform filter window.

    Returns:
    ndarray: The filtered image.
    """
    img_mean = uniform_filter(img, size)
    img_sqr_mean = uniform_filter(img**2, size)
    img_variance = img_sqr_mean - img_mean**2

    overall_variance = np.var(img)

    img_weights = img_variance / (img_variance + overall_variance)
    img_output = img_mean + img_weights * (img - img_mean)
    return img_output

Now that we’ve defined the filter, we can run it on the VV and VH data. You might have noticed that the function takes a size argument. This will change how blurred the image becomes after smoothing. We’ve picked a default value for this analysis, but you can experiement with this if you’re interested.

[13]:

# Define a function to apply the Lee filter to a DataArray
def apply_lee_filter(data_array, size=7):
    """
    Applies the Lee filter to the provided DataArray.

    Parameters:
    data_array (xarray.DataArray): The data array to be filtered.
    size (int): Size of the uniform filter window. Default is 7.

    Returns:
    xarray.DataArray: The filtered data array.
    """
    data_array_filled = data_array.fillna(0)
    filtered_data = xr.apply_ufunc(
        lee_filter, data_array_filled,
        kwargs={"size": size},
        input_core_dims=[["latitude", "longitude"]],
        output_core_dims=[["latitude", "longitude"]],
        dask_gufunc_kwargs={"allow_rechunk": True},
        vectorize=True,
        dask="parallelized",
        output_dtypes=[data_array.dtype]
    )
    return filtered_data

[14]:

# Apply the Lee filter to both VV and VH data
ds_s1["filtered_vv"] = apply_lee_filter(ds_s1.vv, size = 7)
ds_s1["filtered_vh"] = apply_lee_filter(ds_s1.vh, size = 7)

Visualise filtered data

We can now visualise the filtered bands in the same way as the original bands. Note that the filtered values must also be converted to decibels before being displayed.

Visualise filtered VH band

[15]:

# Scale to plot data in decibels
ds_s1["filtered_vh_dB"] = dB_scale(ds_s1.filtered_vh)

# Plot all filtered VH observations for the year
ds_s1.filtered_vh_dB.plot(cmap="Greys_r", robust=True, col="time", col_wrap=5)
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_34_0.png

[16]:

# Plot the average of all filtered VH observations
mean_filtered_vh_dB = ds_s1.filtered_vh_dB.mean(dim="time")

fig = plt.figure(figsize=(7, 9))
mean_filtered_vh_dB.plot(cmap="Greys_r", robust=True)
plt.title("Average filtered VH")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_35_0.png

Visualise filtered VV band

[17]:

# Scale to plot data in decibels
ds_s1["filtered_vv_dB"] = dB_scale(ds_s1.filtered_vv)

# Plot all filtered VV observations for the year
ds_s1.filtered_vv_dB.plot(cmap="Greys_r", robust=True, col="time", col_wrap=5)
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_37_0.png

[18]:

# Plot the average of all filtered VV observations
mean_filtered_vv_dB = ds_s1.filtered_vv_dB.mean(dim="time")

fig = plt.figure(figsize=(7, 9))
mean_filtered_vv_dB.plot(cmap="Greys_r", robust=True)
plt.title("Average filtered VV")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_38_0.png

Now that you’ve finished filtering the data, compare the plots before and after and you should be able to notice the impact of the filtering. If you’re having trouble spotting it, it’s more noticable in the VH band.

Plotting VH and VV histograms

Another way to observe the impact of filtering is to view histograms of the pixel values before and after filtering. Try running the next two cells to view the histograms for VH and VV.

[19]:

fig = plt.figure(figsize=(15, 3))
ds_s1.filtered_vh_dB.plot.hist(bins=1000, label="VH filtered")
ds_s1.vh_dB.plot.hist(bins=1000, label="VH", alpha=0.5)
plt.legend()
plt.xlabel("VH (dB)")
plt.title("Comparison of filtered VH bands to original")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_42_0.png

[20]:

fig = plt.figure(figsize=(15, 3))
ds_s1.filtered_vv_dB.plot.hist(bins=1000, label="VV filtered")
ds_s1.vv_dB.plot.hist(bins=1000, label="VV", alpha=0.5)
plt.legend()
plt.xlabel("VV (dB)")
plt.title("Comparison of filtered VV bands to original")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_43_0.png

You may have noticed that both the original and filtered bands show two peaks in the histogram, which we can classify as a bimodal distribution. Looking back at the band images, it’s clear that the water pixels generally have lower VH and VV values than the land pixels. This lets us conclude that the lower distribution corresponds to water pixels and the higher distribution corresponds to land pixels. Importantly, the act of filtering has made it clear that the two distributions can be separated, which is especially obvious in the VH histogram. This allows us to confidently say that pixel values below a certain threshold are water, and pixel values above it are land. This will form the basis for our classifier in the next section.

Designing a threshold-based water classifier

Given that the distinction between the land and water pixel value distributions is strongest in the VH band, we’ll base our classifier on this distribution. To separate them, we can choose a threshold: pixels with values below the threshold are water, and pixels with values above the threshold are not water (land).

There are a number of ways to determine the threshold; one is to estimate it by looking at the VH histogram. From this, we might guess that \(\text{threshold} = -20.0\) is a reasonable value. Run the cell below to set the threshold.

[21]:

threshold = -20.0

The classifier separates data into two classes: data above the threshold and data below the threshold. In doing this, we assume that values of both segments correspond to the same water and not water distinctions we make visually. This can be represented with a step function:

\[\begin{split}\text{water}(VH) = \left\{ \begin{array}{lr} \text{True} & : VH < \text{threshold}\\ \text{False} & : VH \geq \text{threshold} \end{array} \right.\\\end{split}\]

Visualise threshold

To check if our chosen threshold reasonably divides the two distributions, we can add the threshold to the histogram plots we made earlier. Run the next two cells to view two different visualisations of this.

[22]:

fig = plt.figure(figsize=(15, 3))
plt.axvline(x=threshold, label=f"Threshold at {threshold}", color="red")
ds_s1.filtered_vh_dB.plot.hist(bins=1000, label="VH filtered")
ds_s1.vh_dB.plot.hist(bins=1000, label="VH", alpha=0.5)
plt.legend()
plt.xlabel("VH (dB)")
plt.title("Histogram Comparison of filtered VH bands to original")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_50_0.png

[23]:

fig, ax = plt.subplots(figsize=(15, 3))
ds_s1.filtered_vh_dB.plot.hist(bins=1000, label="VH filtered")
ax.axvspan(xmin=-40.0, xmax=threshold, alpha=0.25, color="green", label="Water")
ax.axvspan(xmin=threshold,
           xmax=-0.5,
           alpha=0.25,
           color="red",
           label="Not Water")
plt.legend()
plt.xlabel("VH (dB)")
plt.title("Effect of the classifier")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_51_0.png

If you’re curious about how changing the threshold impacts the classifier, try changing the threshold value and running the previous two cells again.

Detecting change between two images

The standard deviation we calculated before gives us an idea of how variable a pixel has been over the entire period of time that we looked at. It might also be interesting to look at which pixels have changed between any two particular times in our sample.

In the next cell, we choose the images to compare. Printing the dataset should show you that there are 9 time-steps, so the first has an index value of 0, and the last has an index value of 8. You can change these to be any numbers in between, as long as the start is earlier than the end.

[29]:

start_time_index = 0
end_time_index = ds_s1.water.sizes["time"] - 1

Next, we can define the change as the difference in the classified pixel value at each point. Land becoming water will have a value of -1 and water becoming land will have a value of 1.

[30]:

change = np.subtract(ds_s1.water.isel(time=start_time_index),
                     ds_s1.water.isel(time=end_time_index),
                     dtype=np.float32)

# Set all '0' entries to NaN, which prevents them from displaying in the plot.
change = change.where(change != 0)
ds_s1["change"] = change

Now that we’ve added change to the dataset, you should be able to plot it below to look at which pixels changed. You can also plot the original mean VH composite to see how well the change matches our understanding of the shoreline location.

[31]:

plt.figure(figsize=(15, 12))
ds_s1.filtered_vh_dB.mean(dim="time").plot(cmap="Blues")
ds_s1.change.plot(cmap="RdBu", levels=2)
plt.title(f"Change in pixel value between time={start_time_index} and time={end_time_index}")
plt.show()

../../../_images/notebooks_Real_world_examples_Radar_water_detection_71_0.png

Coastal change or tides?

Tides can greatly affect the appearance of the coastline, particularly along northern Australia where the tidal range is large (up to 12 m). Without additional data, it is difficult to determine whether the change above is due to the coastline having eroded over time, or because the two radar images were captured at different tides (e.g. low vs. high tide). The radar water classifier in this notebook could potentially be combined with tidal modelling from the Coastal Erosion notebook to look into this question in more detail.

Drawing conclusions

Here are some questions to think about:

What are the benefits and drawbacks of the possible classification options we explored?
How could you extend the analysis to extract a shape for the coastline?
How reliable is our classifier?
Is there anything you can think of that would improve it?

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on GitHub.

Last modified: August 2024

Compatible datacube version:

[32]:

print(datacube.__version__)

1.8.18

Detecting water with radar

Background

Sentinel-1 use case

Description

Getting started

Load packages

Connect to the datacube

Analysis parameters

View the selected location

Load and view Sentinel-1 data

Visualise loaded data

Visualise VH band

Visualise VV band

Preprocessing the data through filtering

Speckle Filtering using Lee Filter

Visualise filtered data

Visualise filtered VH band

Visualise filtered VV band

Plotting VH and VV histograms

Designing a threshold-based water classifier

Visualise threshold

Build and apply the classifier

Assessment with mean

Interpreting the mean classification

Assessment with standard deviation

Interpreting the standard deviation of the classification

Detecting change between two images

Coastal change or tides?

Drawing conclusions

Next steps

1. Select another Australian area using DEA Data

2. Investigate Global Areas (excl. Antarctica) using Microsoft Planetary Computer Data

Loading data using `odc-stac`

Additional information

Tags

Detecting water with radar

Background

Sentinel-1 use case

Description

Getting started

Load packages

Connect to the datacube

Analysis parameters

View the selected location

Load and view Sentinel-1 data

Visualise loaded data

Visualise VH band

Visualise VV band

Preprocessing the data through filtering

Speckle Filtering using Lee Filter

Visualise filtered data

Visualise filtered VH band

Visualise filtered VV band

Plotting VH and VV histograms

Designing a threshold-based water classifier

Visualise threshold

Build and apply the classifier

Assessment with mean

Interpreting the mean classification

Assessment with standard deviation

Interpreting the standard deviation of the classification

Detecting change between two images

Coastal change or tides?

Drawing conclusions

Next steps

1. Select another Australian area using DEA Data

2. Investigate Global Areas (excl. Antarctica) using Microsoft Planetary Computer Data

Loading data using odc-stac

Additional information

Tags

Loading data using `odc-stac`