Creating Water Quality Maps From Remote Sensed Images With Python

Learn how to use waterquality python package to create beautiful water quality maps from Sentinel 2 imagery.

Image by author.

Introduction

Water security is one of the sustainable goals for 2030, set by the United Nations — UN [1]. Having water available in proper quantity and quality is essential for life on Earth. There are estimations that the world’s population can reach 10 billion people by 2050. Moreover, climate changes are also pushing more pressure in water availability, increasing water-related hazards, such as droughts.

Improving the water governance is essential, and stakeholders need to be fed with data to take good decisions. Water quality parameters are costly to be obtained, as it needs to be obtained in the field, and they are usually punctual, difficult to be generalized spatially through large areas, or along a river, for example. Many water bodies lack water quality measurements, especially in developing countries.

In this context, the use of satellite sensors is emerging as one viable solution to measure some water quality parameters from space. The copernicus mission from European Space Agency, for example provides medium resolution images with up to 10m, every 4–5 days. Spectral resolution is also improving as new satellites are being launched.

In my previous post Water Detection in High Resolution Satellite Images using the waterdetect python package I’ve shown how to use the waterdetect Python package to extract water masks from satellite imagery. In the current story I will present how to use another package (waterquality package) that makes it easier the creation of water quality maps from remote sensed images.

Installation

The waterquality package is available at the following github repository: https://github.com/cordmaur/WaterQuality

Before installing it, make sure you have an active environment with GDAL available.

To install waterquality using pip:

pip install git+https://github.com/cordmaur/WaterQuality.git@main

Alternatively, you can clone the project and install it in edit mode

git clone https://github.com/cordmaur/WaterQuality.git
cd WaterQuality
pip install -e .

Once installed, you can test the installation by typing waterquality -h and it will display the basic usage from console. Here we will call the package from a Jupyter Notebook.

From Jupyter

The first thing we will need to run the package is to prepare the configuration files for the waterdetect and for the waterquality. In these files it is possible to configure all the processing. I’ve prepared two samples that can be accessed in the test folder of the project. So, if you are in the waterdetect root folder, the following commands should work.

Then, we have to create an instance of the WaterQuality class and point to the correct folders. For this example it is necessary to have at least one image (Sentinel-2 or Landsat 8) downloaded and uncompressed. I will use a Sentinel-2 L2-A image obtained from ESA. These images are processed with the Sen2Cor algorithm for atmospheric correction, so I am using product = ‘S2_S2COR’. Additionally, I specified a smaller area to focus on my reservoir. This can be done by passing a shapefile with one polygon as parameter to the creation function. Also, as I am processing just one image, I selected single_mode=True. The command to create the instance is in the following snippet:

Code output:

Initializing waterquality package (version=1.1.0) 
Waterdetect version: 1.5.12 
Loading configuration file ..\tests\aux_data\WaterQuality.ini 
File ..\tests\aux_data\WaterQuality.ini verified. 
Loading configuration file ..\tests\aux_data\WaterDetect.ini 
File ..\tests\aux_data\WaterDetect.ini verified. 
Folder d:/images/22KEV/S2A_MSIL2A_20210706T133231_N0301_R081_T22KEV_20210706T172741.SAFE verified. 
File d:/images/shp/22KEV.shp verified. Folder d:/images/output verified.

You should see a message at the beginning indicating that the package has been initialized and that waterdetect has been found. It will also check for the informed folders.

Selecting the Inversion Function

Before launching the processing, it is necessary to create the functions that will be responsible for inverting from the spectral signature to actual water quality parameters (called inversion function).

These functions can be found in the literature, and they will depend on various factors such as water type, satellite sensor, atmospheric conditions, among others.

Once the correct function is selected, it is necessary to inform waterquality how to use this function to perform the inversion. You can inform more than one function for different parameters, for example, or to test different functions in the same area.

As an example, I will implement a function for total suspended matters (TSM) proposed by Nechad et al. (2010) [2]. This is the function proposed:

The functions shall be informed in a dictionary containing {‘name f1’: function 1, ‘name f2’: function 2, …}

NOTE: Before using an inversion function, it is important to check the validity of the function in the region to be analyzed, preferably by comparing the results with field measurements.

Tables 1 and 4 from article lists the parameters A, B and C to be used, depending on the sensor wavelength being used. To select the wavelength we can see the best coefficient of determination R², that in this case is 82.9%, around 710nm. Then, we can check on the satellite’s bands if there is any band close to this wavelength.

Figure 1: Sentinel-2 bands with 20m spatial resolution. Source: ESA.

From Figure 1, we can see that Band 5, which is 705nm is the closest to 710nm, so we will get the parameters according to this wavelength. In this case we will have: A = 493.65, B = 1.16 and C = 0.188.

Implementing the Inversion Function

Now that we have a function and its parameters, let’s implement it in Python and pass it to our waterquality object. The function to be implemented can use any band available in the satellite image but, as the software is generic for different satellites and atmospheric processing routines, the name must match the name used internally by the waterdetect package. So, before writing the inversion function, check the bands available for the satellite, that are informed in the InputOutput.py module from waterdetect.

{'bands_names': {'Blue': 'B02_10m', 'Green': 'B03_10m', 'Red': 'B04_10m', 'Mir': 'B11_20m', 'Mir2': 'B12_20m', 'RedEdg1': 'B05_20m', 'RedEdg2': 'B06_20m', 'RedEdg3': 'B07_20m', 'Nir': 'B08_10m', 'Nir2': 'B8A_20m'}, 'suffix': '.jp2', 'string': '', 'metadata': '*MTD_TL.xml', 'subdir': 'GRANULE/*/IMG_DATA', 'recursive': True}

From the result, we can see that the Band5 is called RedEdg1 inside waterdetect. So, to make use of this band, it enough to use this name in the function arguments. To simplify the calculation of the quality parameter, all the bands will be rescaled to the same shape, as defined in the reference_band in the config file WaterDetect.ini. Our function will look like this:

Note that we don’t have to bother with no data pixels or water vs land pixels, as the package will take care of these for us. Now, we have to tell the package to run the processing using this function as inversion function. This can be done by passing a dictionary to the run_water_quality method.

Starting WaterDetection version: 1.5.11 
Opening image in folder d:\images\22KEV\S2A_MSIL2A_20210706T133231_N0301_R081_T22KEV_20210706T172741.SAFE 
Retrieving bands for product S2_S2COR 35 were found: 
Loading band: Blue 
Loading band: Green 
Loading band: Red 
Loading band: Mir 
Loading band: Mir2 
Loading band: RedEdg1
...

Displaying the results

Once the processing has finished, the results will be available and we can plot the results on the jupyter notebook. If PDF Reports is set to True in the WaterDetect.ini, a report will be written in the output folder with all the information about the processing (Figure 2). The Report will show the original image (RGB), the water mask, the sunglint probability and the water quality parameter with a color scale.

Figure 2: Water quality report example. The graph (bottom left) represents the clustering used to identify the water pixels in the scene. Image by author.

From the notebook, the results can be accessed through the quality_parameters member. The function plot_param can also be used to create a quick view of the result:

wq.plot_param('SPM', cmap='jet', vmin=5, vmax=20)

Figure 3: Code output. Image by author.

Conclusion

The number of studies in the literature towards water quality assessment by the remote sensed images has been increasing in the last years. However, we can note that it is still lacking operational tools to make it available to a broader audience. The softwares are usually developed within specific researches and stay experimental, locked in the universities or research centers. What I hope that from this initiative, and others like the waterdetect package is that more and more tools for operational remote sensing are available for people to use, to study and to share.

Thanks and see you in the next story.