Visualization and Analysis of Brazil Floods 2024¶
Introduction¶
The 2024 Rio Grande do Sul floods are severe floods caused by heavy rains and storms that have hit the Brazilian state of Rio Grande do Sul, and the adjacent Uruguayan cities of Treinta y Tres, Paysandú, Cerro Largo, and Salto. From 29 April 2024 through to May 2024, it resulted in over 140 fatalities, widespread landslides, and a dam collapse. It is considered the country's worst flooding in over 80 years. See the Wikipedia page for more information about the 2024 Brazil floods.
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 Pakistan 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 geemap
Import libraries¶
Import the earthengine-api and geemap.
import ee
import geemap.foliumap as geemap
geemap.ee_initialize()
url = "https://github.com/opengeos/datasets/releases/download/places/Rio_Grande_do_Sul_Brazil.geojson"
roi = geemap.geojson_to_ee(url)
# roi
Downloading... From: https://github.com/opengeos/datasets/releases/download/places/Rio_Grande_do_Sul_Brazil.geojson To: /tmp/18896865-507c-4806-bc74-6738b5fd4604.geojson 100%|███████████████████████████████████████████████████████| 61.2k/61.2k [00:00<00:00, 2.24MB/s]
centroid = roi.geometry().centroid(1)
lon, lat = centroid.getInfo()["coordinates"]
print(f"Longitude: {lon}, Latitude: {lat}")
Longitude: -53.2453436538024, Latitude: -29.778651228599525
Create an interactive map¶
Specify the center point [lat, lon] and zoom level of the map.
m = geemap.Map()
style = {"fillColor": "00000000", "color": "FF0000"}
m.add_layer(roi.style(**style), {}, "ROI")
geom = roi.geometry()
m.center_object(geom, 6)
m
In the tutorial, we will focus on Rio Grande do Sul in Brazil, but the code can be easily modified to visualize and analyze floods in other countries. Modify the place_name 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.
place_name = "Rio Grande do Sul"
pre_flood_start_date = "2024-01-01"
pre_flood_end_date = "2024-04-27"
post_flood_start_date = "2024-04-29"
post_flood_end_date = "2024-05-17"
Create Landsat composites¶
Create a Landsat 8 composite for the pre-flood period (August 1 to September 30, 2021) using the USGS Landsat 8 Collection 2 Tier 1 Raw Scenes.
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", 15))
)
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: 374
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: 34
Visualize the Landsat 8 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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.1 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.1
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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 1000 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 kilometers. 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=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_pre_flood)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 16801.019608 |
area_2022 = geemap.zonal_stats(
water_post.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2022)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 1884.207843 |
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(flood_area)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 634.188235 |
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 = "2023-10-01"
pre_flood_end_date = "2024-04-27"
post_flood_start_date = "2024-04-29"
post_flood_end_date = "2024-05-17"
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: 161
Create the Sentinel-1 image collection for the 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: 19
Create Sentinel-1 SAR composites for the pre-flood and flood periods.
m = geemap.Map()
m.add_basemap("SATELLITE")
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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("SATELLITE")
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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 -18 dB.
threshold = -18
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
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.add_layer(roi.style(**style), {}, place_name)
m.center_object(roi, 6)
m
Calculate SAR flood area¶
area_pre_flood = geemap.zonal_stats(
water_pre.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_pre_flood)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 10160.776471 |
area_2022 = geemap.zonal_stats(
water_post.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2022)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 13054.921569 |
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(flood_area)
Computing statistics ...
| area_1 | area_tot_g | codigo | geocodigo | mapid | mslink | nome | perimetro_ | sum | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 281748.155 | 281748.155 | 43 | 43 | 99 | 223 | RIO GRANDE DO SUL | 3524.58571 | 5745.360784 |