Visualization and Analysis of Dongting Lake Floods 2024¶
Introduction¶
Dongting Lake, China's second-largest freshwater lake, burst its embankment on July 5, 2024, flooding more than 47 sq km (18 square miles). The floods have affected more than 364,000 people, inundated farmland, damaged over 2,600 sections of roads and bridges and caused multiple landslides in Pingjiang. For more information, please refer to the South China Morning Post.
Requirements¶
To follow this tutorial, you must first sign up for a Google Earth Engine account. Earth Engine is a cloud computing platform with a multi-petabyte catalog of satellite imagery and geospatial datasets. It is free for noncommercial use. To authenticate the Earth Engine Python API, see instructions here.
In this tutorial, we will use the geemap Python package to visualize and analyze the Dongting Lake floods. Geemap enables users to analyze and visualize Earth Engine datasets interactively within a Jupyter-based environment with minimal coding. To learn more about geemap, check out https://geemap.org.
Installation¶
Uncomment the following line to install geemap if needed.
# %pip install -U geemap
Import libraries¶
Import the earthengine-api and geemap.
import ee
import geemap.foliumap as geemap
geemap.ee_initialize()
Create an interactive map¶
m = geemap.Map(center=[29.343875, 112.986832], zoom=11)
m
Specify ROI and time range¶
In the tutorial, we will focus on the Dongting Lake in China, but the code can be easily modified to visualize and analyze floods in other countries or regions. Modify the geom variable to specify the place of interest and set the date range for the flood event. In order to extract the flood extent, we also need to specify the date range for the pre-flood period.
m = geemap.Map()
m.add_basemap("HYBRID")
geom = ee.Geometry.Polygon(
[
[
[112.755432, 29.642708],
[112.359924, 28.591757],
[113.414612, 28.586933],
[113.755188, 29.652256],
[112.755432, 29.642708],
]
]
)
roi = ee.FeatureCollection(geom)
jrc = ee.Image("JRC/GSW1_4/GlobalSurfaceWater").clipToCollection(roi)
vis_params = {
"bands": ["occurrence"],
"min": 0.0,
"max": 100.0,
"palette": ["ffffff", "ffbbbb", "0000ff"],
}
m.add_layer(jrc, vis_params, "Water Occurrence")
m.add_layer(geom, {}, "ROI")
m.center_object(geom)
m
pre_flood_start_date = "2024-01-01"
pre_flood_end_date = "2024-07-04"
post_flood_start_date = "2024-07-05"
post_flood_end_date = "2024-07-17"
Create Landsat composites¶
Create a Landsat composite for the pre-flood period using the Harmonized Landsat Data.
pre_flood_col = (
ee.ImageCollection("NASA/HLS/HLSL30/v002")
.filterBounds(roi)
.filterDate(pre_flood_start_date, pre_flood_end_date)
.filter(ee.Filter.lt("CLOUD_COVERAGE", 25))
)
print(
f"The number of images in the pre-flood collection: {pre_flood_col.size().getInfo()}"
)
The number of images in the pre-flood collection: 20
Create a Landsat composite for the post-flood period using the Harmonized Landsat Data.
post_flood_col = (
ee.ImageCollection("NASA/HLS/HLSL30/v002")
.filterBounds(roi)
.filterDate(post_flood_start_date, post_flood_end_date)
.filter(ee.Filter.lt("CLOUD_COVERAGE", 50))
)
print(
f"The number of images in the post-flood collection: {post_flood_col.size().getInfo()}"
)
The number of images in the post-flood collection: 5
Visualize the Landsat composite for the pre-flood and flood periods.
m = geemap.Map()
pre_flood_image = pre_flood_col.median().clipToCollection(roi)
post_flood_image = post_flood_col.median().clipToCollection(roi)
vis_params = {"bands": ["B6", "B5", "B4"], "min": 0, "max": 0.4}
m.add_layer(pre_flood_image, vis_params, "Landsat Pre-flood")
m.add_layer(post_flood_image, vis_params, "Landsat Post-flood")
m.center_object(roi, 9)
m
Compare Landsat composites side by side¶
Compare the pre-flood and flood composites side by side.
m = geemap.Map()
left_layer = geemap.ee_tile_layer(pre_flood_image, vis_params, "Landsat Pre-flood")
right_layer = geemap.ee_tile_layer(post_flood_image, vis_params, "Landsat Post-flood")
m.split_map(
left_layer,
right_layer,
left_label="Landsat Pre-flood",
right_label="Landsat Post-flood",
)
m.center_object(roi, 9)
m
Compute Normalized Difference Water Index (NDWI)¶
The Normalized Difference Water Index (NDWI) is a commonly used index for detecting water bodies. It is calculated as follows:
$$NDWI = \frac{Green - NIR}{Green + NIR}$$
where Green is the green band and NIR is the near-infrared band. The NDWI values range from -1 to 1. The NDWI values are usually thresholded to a positive number (e.g., 0.1-0.3) to identify water bodies.
Landsat 8 imagery has 11 spectral bands. The Landsat 8 NDWI is calculated using the green (B3) and NIR (B5) bands.
ndwi_pre = pre_flood_image.normalizedDifference(["B3", "B5"]).rename("NDWI")
ndwi_post = post_flood_image.normalizedDifference(["B3", "B5"]).rename("NDWI")
Compute the NDWI layers for the pre-flood and flood periods side by side.
m = geemap.Map()
ndwi_vis = {"min": -1, "max": 1, "palette": "ndwi"}
left_layer = geemap.ee_tile_layer(ndwi_pre, ndwi_vis, "NDWI Pre-flood")
right_layer = geemap.ee_tile_layer(ndwi_post, ndwi_vis, "NDWI Post-flood")
m.split_map(
left_layer, right_layer, left_label="NDWI Pre-flood", right_label="NDWI Post-flood"
)
m.center_object(roi, 9)
m
Extract Landsat water extent¶
To extract the water extent, we need to convert the NDWI images to binary images using a threshold value. The threshold value is usually set to 0 to 0.3. The smaller the threshold value, the more water bodies will be detected, which may increase the false positive rate. The larger the threshold value, the fewer water bodies will be detected, which may increase the false negative rate.
threshold = -0.05
water_pre = ndwi_pre.gt(threshold)
water_post = ndwi_post.gt(threshold)
Combine the pre-flood and surface water extent side by side.
m = geemap.Map()
m.add_layer(pre_flood_image, vis_params, "Landsat Pre-flood", True)
m.add_layer(post_flood_image, vis_params, "Landsat Post-flood", True)
left_layer = geemap.ee_tile_layer(
water_pre.selfMask(), {"palette": "blue"}, "Water Pre-flood"
)
right_layer = geemap.ee_tile_layer(
water_post.selfMask(), {"palette": "red"}, "Water Post-flood"
)
m.split_map(
left_layer,
right_layer,
left_label="Water Pre-flood",
right_label="Water Post-flood",
)
m.center_object(roi, 9)
m
Extract Landsat flood extent¶
To extract the flood extent, we need to subtract the pre-flood water extent from the flood water extent. The flood extent is the difference between the flood water extent and the pre-flood water extent. In other words, pixels identified as water in the flood period but not in the pre-flood period are considered as flooded pixels. The selfMask() method is used to mask out the no-data pixels.
flood_extent = water_post.subtract(water_pre).gt(0).selfMask()
Add the flood extent layer to the map.
m = geemap.Map()
m.add_layer(pre_flood_image, vis_params, "Landsat Pre-flood", True)
m.add_layer(post_flood_image, vis_params, "Landsat Post-flood", True)
left_layer = geemap.ee_tile_layer(
water_pre.selfMask(), {"palette": "blue"}, "Water Pre-flood"
)
right_layer = geemap.ee_tile_layer(
water_post.selfMask(), {"palette": "red"}, "Water Post-flood"
)
m.split_map(
left_layer,
right_layer,
left_label="Water Pre-flood",
right_label="Water Post-flood",
)
m.add_layer(flood_extent, {"palette": "cyan"}, "Flood Extent")
m.center_object(roi, 9)
m
Calculate Landsat flood area¶
To calculate the flood area, we can use the geemap.zonal_stats() function. The required input parameters are the flood extent layer and the country boundary layer. The scale parameter can be set to 30 to specify the spatial resolution of image to be used for calculating the zonal statistics. The stats_type parameter can be set to SUM to calculate the total area of the flood extent in square meters. Set return_fc=True to return the zonal statistics as an ee.FeatureCollection object, which can be converted to a Pandas dataframe.
area_pre_flood = geemap.zonal_stats(
water_pre.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
area_pre_df = geemap.ee_to_df(area_pre_flood)
area_pre_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 1.357662e+06 |
print(f"Pre-flood water area: {round(area_pre_df['sum'].values[0] * 900 / 1e6)} km2")
Pre-flood water area: 1222 km2
area_post_flood = geemap.zonal_stats(
water_post.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
area_post_df = geemap.ee_to_df(area_post_flood)
area_post_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 3.328183e+06 |
print(f"Post-flood water area: {round(area_post_df['sum'].values[0] * 900 / 1e6)} km2")
Post-flood water area: 2995 km2
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
flood_area_df = geemap.ee_to_df(flood_area)
flood_area_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 2.182365e+06 |
print(f"Flooded area: {round(flood_area_df['sum'].values[0] * 900 / 1e6)} km2")
Flooded area: 1964 km2
Create Sentinel-1 SAR composites¶
Besides Landsat, we can also use Sentinel-1 Synthetic Aperture Radar (SAR) data to extract flood extent. Radar can collect signals in different polarizations, by controlling the analyzed polarization in both the transmit and receive paths. Signals emitted in vertical (V) and received in horizontal (H) polarization would be indicated by a VH. Alternatively, a signal that was emitted in horizontal (H) and received in horizontal (H) would be indicated by HH, and so on. Examining the signal strength from these different polarizations carries information about the structure of the imaged surface. Rough surface scattering, such as that caused by bare soil or water, is most sensitive to VV scattering. Therefore, VV polarization is often used to detect water bodies.
Sentinel-1 operates in four exclusive acquisition modes:
- Stripmap (SM)
- Interferometric Wide swath (IW)
- Extra-Wide swath (EW)
- Wave mode (WV)
The Interferometric Wide swath (IW) mode allows combining a large swath width (250 km) with a moderate geometric resolution (5 m by 20 m). The IW mode is the default acquisition mode over land. In this tutorial, we will use Sentinel-1 IW mode data to extract flood extent.
The Sentinel-1 SAR data are available from 2014 to present. Let's filter the COPERNICUS/S1_GRD dataset by the date range and location.
pre_flood_start_date = "2024-06-01"
pre_flood_end_date = "2024-06-30"
post_flood_start_date = "2024-07-05"
post_flood_end_date = "2024-07-17"
Create the Sentinel-1 image collection for the pre-flood period.
s1_col_pre = (
ee.ImageCollection("COPERNICUS/S1_GRD")
.filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV"))
.filter(ee.Filter.eq("instrumentMode", "IW"))
.filterDate(pre_flood_start_date, pre_flood_end_date)
.filterBounds(roi)
.select("VV")
)
print(
f"The number of images in the pre-flood collection: {s1_col_pre.size().getInfo()}"
)
The number of images in the pre-flood collection: 4
Create the Sentinel-1 image collection for the post-flood period.
s1_col_post = (
ee.ImageCollection("COPERNICUS/S1_GRD")
.filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV"))
.filter(ee.Filter.eq("instrumentMode", "IW"))
# .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.filterDate(post_flood_start_date, post_flood_end_date)
.filterBounds(roi)
.select("VV")
)
print(
f"The number of images in the post-flood collection: {s1_col_post.size().getInfo()}"
)
The number of images in the post-flood collection: 1
Create Sentinel-1 SAR composites for the pre- and post-flood periods.
m = geemap.Map()
m.add_basemap("HYBRID")
sar_pre = s1_col_pre.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)
sar_post = s1_col_post.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)
m.add_layer(sar_pre, {"min": -25, "max": -5}, "SAR Pre-flood")
m.add_layer(sar_post, {"min": -25, "max": -5}, "SAR Post-flood")
m.center_object(roi, 9)
m
Apply speckle filtering¶
Speckle, appearing in synthetic aperture radar (SAR) images as granular noise, is due to the interference of waves reflected from many elementary scatterers. Speckle in SAR images complicates the image interpretation problem by reducing the effectiveness of image segmentation and classification (Lee et al., 1994). Therefore, speckle filtering is often applied to SAR images to reduce the speckle noise. In this example, we apply a morphological speckle filter to the Sentinel-1 SAR images. The morphological speckle filter is a non-linear filter that uses the median value of a pixel and its neighboring pixels to replace the pixel value. The kernel size is set to 100 meters.
s1_col_pre = s1_col_pre.map(lambda img: img.focal_median(100, "circle", "meters"))
s1_col_post = s1_col_post.map(lambda img: img.focal_median(100, "circle", "meters"))
m = geemap.Map()
m.add_basemap("HYBRID")
sar_pre = s1_col_pre.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)
sar_post = s1_col_post.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)
m.add_layer(sar_pre, {"min": -25, "max": -5}, "SAR Pre-flood")
m.add_layer(sar_post, {"min": -25, "max": -5}, "SAR Post-flood")
m.center_object(roi, 9)
m
Compare Sentinel-1 SAR composites side by side¶
Create a split-view map to compare the pre-flood and flood SAR composites side by side.
m = geemap.Map()
left_layer = geemap.ee_tile_layer(sar_pre, {"min": -25, "max": -5}, "SAR Pre-flood")
right_layer = geemap.ee_tile_layer(sar_post, {"min": -25, "max": -5}, "SAR Post-flood")
m.split_map(
left_layer,
right_layer,
left_label="Sentinel-1 Pre-flood",
right_label="Sentinel-1 Post_flood",
)
m.center_object(roi, 9)
m
Extract SAR water extent¶
Water usually appears dark in SAR images because radar waves are reflected differently by different surfaces. Water is a smooth, flat surface that does not reflect radar waves very well, so it appears dark in SAR images. Thresholding SAR imagery is one of the most widely used approaches to delineate water extent for its effectiveness and efficiency (Liang and Liu, 2020). Thresholding methods can be generally divided into two categories: global and local. Global thresholding methods use a single threshold value to segment the entire image. Local thresholding methods use a different threshold value for each pixel. In this example, we use a global thresholding method to extract the water extent. The threshold value is set to -16 dB.
threshold = -16
water_pre = sar_pre.lt(threshold)
water_post = sar_post.lt(threshold)
Create a split-view map to compare the pre-flood and flood water extent side by side.
m = geemap.Map()
m.add_layer(sar_pre, {"min": -25, "max": -5}, "SAR Pre-flood")
m.add_layer(sar_post, {"min": -25, "max": -5}, "SAR Post-flood")
left_layer = geemap.ee_tile_layer(
water_pre.selfMask(), {"palette": "blue"}, "Water Pre-flood"
)
right_layer = geemap.ee_tile_layer(
water_post.selfMask(), {"palette": "red"}, "Water Post-flood"
)
m.split_map(
left_layer,
right_layer,
left_label="Water Pre-flood",
right_label="Water Post-flood",
)
m.center_object(roi, 9)
m
Extract SAR flood extent¶
Similar to the Landsat approach, we can subtract the pre-flood water extent from the flood water extent to extract the flood extent.
flood_extent = water_post.subtract(water_pre).gt(0).selfMask()
The flood extent is the difference between the flood water extent and the pre-flood water extent. In other words, pixels identified as water in the flood period but not in the pre-flood period are considered as flooded pixels, which are shown in cyan.
m = geemap.Map()
m.add_layer(sar_pre, {"min": -25, "max": -5}, "SAR Pre-flood")
m.add_layer(sar_post, {"min": -25, "max": -5}, "SAR Post-flood")
left_layer = geemap.ee_tile_layer(
water_pre.selfMask(), {"palette": "blue"}, "Water Pre-flood"
)
right_layer = geemap.ee_tile_layer(
water_post.selfMask(), {"palette": "red"}, "Water Post-flood"
)
m.split_map(
left_layer,
right_layer,
left_label="Water Pre-flood",
right_label="Water Post-flood",
)
m.add_layer(flood_extent, {"palette": "cyan"}, "Flood Extent")
m.center_object(roi, 9)
m
Calculate SAR flood area¶
area_pre_flood = geemap.zonal_stats(
water_pre.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
area_pre_df = geemap.ee_to_df(area_pre_flood)
area_pre_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 2.333509e+06 |
print(f"Pre-flood water area: {round(area_pre_df['sum'].values[0] * 900 / 1e6)} km2")
Pre-flood water area: 2100 km2
area_post_flood = geemap.zonal_stats(
water_post.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
area_post_df = geemap.ee_to_df(area_post_flood)
area_post_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 2.454894e+06 |
print(f"Post-flood water area: {round(area_post_df['sum'].values[0] * 900 / 1e6)} km2")
Post-flood water area: 2209 km2
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), roi, scale=30, stat_type="SUM", return_fc=True
)
flood_area_df = geemap.ee_to_df(flood_area)
flood_area_df
Computing statistics ...
| sum | |
|---|---|
| 0 | 550634.537255 |
print(f"Flooded area: {round(flood_area_df['sum'].values[0] * 900 / 1e6)} km2")
Flooded area: 496 km2