Lab3b - 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.
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:
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.
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
# import geemap.foliumap as geemap
c:\Users\C00553090\AppData\Local\miniconda3\envs\hypercoast\lib\site-packages\pandas\core\computation\expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).
from pandas.core.computation.check import NUMEXPR_INSTALLED
# This function is used to initialize the Google Earth Engine API,
# which is necessary before you start working with any Earth Engine data.
geemap.ee_initialize()
Download administrative boundaries#
Download the administrative boundaries of Rio Grande doSul, Brazil from here.
# Specifying a URL and the name of the file being accessed.
# The file Rio_Grande_do_Sul_Brazil.geojson is likely a GEOJSON file containing geographical data specific to Rio Grande do Sul, a state in Brazil.
# GEOJSON is a format for encoding various geographical data structures.
url = "https://github.com/opengeos/datasets/releases/download/places/Rio_Grande_do_Sul_Brazil.geojson"
# geemap refers to the geemap library
# geojson_to_ee(url) is used to convert a GEOJSON file into a Google Earth Engine object.
# The url variable, the GEOJSON file from the URL is converted into a format (an Earth Engine object) that can be used within the Google Earth Engine framework.
#roi = geemap.geojson_to_ee(url)
roi = ee.FeatureCollection(ee.Geometry.BBox(-91.9264, 28.6081, -88.3441, 30.6893))
# roi
# geometry(): This method extracts the geometry from the roi object.
# In the context of Google Earth Engine, geometry represents the shapes (points, lines, polygons) that define the region.
# centroid(1): This method calculates the geometric center (centroid) of the region.
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.
# This line initializes a new interactive map using the geemap library.
# Map() is a class provided by geemap that creates an interactive map widget, which can be displayed within a Jupyter notebook or a Python environment supporting interactive widgets.
m = geemap.Map()
style = {"fillColor": "00000000", "color": "FF0000"}
# m.add_layer(...): This method adds the styled roi as a new layer to the map m.
# The empty dictionary {} could be used for additional layer-specific options if needed, and "ROI" is the name given to this layer, which will appear in the map's layer control.
m.add_layer(roi.style(**style), {}, "ROI")
# extract geometry of the roi object.
geom = roi.geometry()
m.center_object(geom, 6)
m
m.user_roi_coords()
[-91.9264, 28.6081, -88.3441, 30.6893]
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-09-30"
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: 344
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:
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
# .gt(threshold): This is a method applied to the NDWI data. gt stands for "greater than".
# This method checks each value in ndwi_pre to see if it is greater than the threshold (0.1 in this case).
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": "yellow"}, "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": "yellow"}, "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)
area_2022 = geemap.zonal_stats(
water_post.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2022)
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), roi, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(flood_area)