Visualization and Analysis of Pakistan Floods 2022 Using Earth Engine and Geemap¶
Introduction¶
From 14 June to October 2022, floods in Pakistan killed 1,739 people, and caused 3.2 trillion Pakistani Rupees (14.9 billion US Dollars) of damage and 3.3 trillion Pakistani Rupees (15.2 billion US Dollars) of economic losses. The flooding was the world's deadliest flood since the 2020 South Asian floods and described as the worst in the country's history. On 25 August 2022, Pakistan declared a state of emergency because of the flooding. See the Wikipedia page for more information about the 2022 Pakistan 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
Create an interactive map¶
Specify the center point [lat, lon] and zoom level of the map.
m = geemap.Map(center=[29.3055, 68.9062], zoom=6)
m
Search datasets¶
Click on the globe icon in the top left corner of the map to open the search panel. Select the data tab and enter a keyword to search for datasets, e.g., countries. Press Enter to search. The search results will populate the dropdown list. Select the dataset you want to add to the map from the dropdown list, such as the LSIB 2017: Large Scale International Boundary Polygons, Simplified. Click on the import button to add a new code cell in Jupyter Notebook. Note that Google Colab and JupyterLab do not support creating new code cells programmatically. You will need to manually add a new code cell and copy the data sample code to the new cell.
In the tutorial, we will focus on Pakistan, but the code can be easily modified to visualize and analyze floods in other countries. Modify the country_name variable to specify the country 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.
country_name = "Pakistan"
pre_flood_start_date = "2021-08-01"
pre_flood_end_date = "2021-09-30"
flood_start_date = "2022-08-01"
flood_end_date = "2022-09-30"
Visualize datasets¶
Specify the country of interest and filter the dataset by the country name.
m = geemap.Map()
country = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017").filter(
ee.Filter.eq("country_na", country_name)
)
style = {"color": "black", "fillColor": "00000000"}
m.add_layer(country.style(**style), {}, country_name)
m.center_object(country, 5)
m
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.
m = geemap.Map()
landsat_col_2021 = (
ee.ImageCollection("LANDSAT/LC08/C02/T1")
.filterDate(pre_flood_start_date, pre_flood_end_date)
.filterBounds(country)
)
landsat_2021 = ee.Algorithms.Landsat.simpleComposite(landsat_col_2021).clipToCollection(
country
)
vis_params = {"bands": ["B6", "B5", "B4"], "max": 128}
m.add_layer(landsat_2021, vis_params, "Landsat 2021")
Create a Landsat 8 composite for the flood period (August 1 to September 30, 2022).
landsat_col_2022 = (
ee.ImageCollection("LANDSAT/LC08/C02/T1")
.filterDate(flood_start_date, flood_end_date)
.filterBounds(country)
)
landsat_2022 = ee.Algorithms.Landsat.simpleComposite(landsat_col_2022).clipToCollection(
country
)
m.add_layer(landsat_2022, vis_params, "Landsat 2022")
m.center_object(country, 5)
m
Compare Landsat composites side by side¶
Combine the pre-flood and flood composites side by side.
m = geemap.Map()
m.setCenter(68.4338, 26.4213, 7)
left_layer = geemap.ee_tile_layer(landsat_2021, vis_params, "Landsat 2021")
right_layer = geemap.ee_tile_layer(landsat_2022, vis_params, "Landsat 2022")
m.split_map(
left_layer, right_layer, left_label="Landsat 2021", right_label="Landsat 2022"
)
m.add_layer(country.style(**style), {}, country_name)
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_2021 = landsat_2021.normalizedDifference(["B3", "B5"]).rename("NDWI")
ndwi_2022 = landsat_2022.normalizedDifference(["B3", "B5"]).rename("NDWI")
Compute the NDWI layers for the pre-flood and flood periods side by side.
m = geemap.Map()
m.setCenter(68.4338, 26.4213, 7)
ndwi_vis = {"min": -1, "max": 1, "palette": "ndwi"}
left_layer = geemap.ee_tile_layer(ndwi_2021, ndwi_vis, "NDWI 2021")
right_layer = geemap.ee_tile_layer(ndwi_2022, ndwi_vis, "NDWI 2022")
m.split_map(left_layer, right_layer, left_label="NDWI 2021", right_label="NDWI 2022")
m.add_layer(country.style(**style), {}, country_name)
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_2021 = ndwi_2021.gt(threshold)
water_2022 = ndwi_2022.gt(threshold)
Combine the pre-flood and surface water extent side by side.
m = geemap.Map()
m.setCenter(68.4338, 26.4213, 7)
m.add_layer(landsat_2021, vis_params, "Landsat 2021", False)
m.add_layer(landsat_2022, vis_params, "Landsat 2022", False)
left_layer = geemap.ee_tile_layer(
water_2021.selfMask(), {"palette": "blue"}, "Water 2021"
)
right_layer = geemap.ee_tile_layer(
water_2022.selfMask(), {"palette": "red"}, "Water 2022"
)
m.split_map(left_layer, right_layer, left_label="Water 2021", right_label="Water 2022")
m.add_layer(country.style(**style), {}, country_name)
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_2022.subtract(water_2021).gt(0).selfMask()
Add the flood extent layer to the map.
m = geemap.Map()
m.setCenter(68.4338, 26.4213, 7)
m.add_layer(landsat_2021, vis_params, "Landsat 2021", False)
m.add_layer(landsat_2022, vis_params, "Landsat 2022", False)
left_layer = geemap.ee_tile_layer(
water_2021.selfMask(), {"palette": "blue"}, "Water 2021"
)
right_layer = geemap.ee_tile_layer(
water_2022.selfMask(), {"palette": "red"}, "Water 2022"
)
m.split_map(left_layer, right_layer, left_label="Water 2021", right_label="Water 2022")
m.add_layer(flood_extent, {"palette": "cyan"}, "Flood Extent")
m.add_layer(country.style(**style), {}, country_name)
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_2021 = geemap.zonal_stats(
water_2021.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2021)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 4205.678431 | S Asia |
area_2022 = geemap.zonal_stats(
water_2022.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2022)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 13145.027451 | S Asia |
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(flood_area)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 11065.72549 | S Asia |
The total area of the flood extent is 11,065 square kilometers based on Landsat 8 images.
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.
s1_col_2021 = (
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(pre_flood_start_date, pre_flood_end_date)
.filterBounds(country)
.select("VV")
)
247 Sentinel-1 IW mode images are available for the flood period.
# s1_col_2021
Create the Sentinel-1 image collection for the flood period.
s1_col_2022 = (
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(flood_start_date, flood_end_date)
.filterBounds(country)
.select("VV")
)
250 Sentinel-1 IW mode images are available for the flood period.
# s1_col_2022
Create Sentinel-1 SAR composites for the pre-flood and flood periods.
m = geemap.Map()
m.add_basemap("HYBRID")
sar_2021 = s1_col_2021.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
sar_2022 = s1_col_2022.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
m.add_layer(sar_2021, {"min": -25, "max": -5}, "SAR 2021")
m.add_layer(sar_2022, {"min": -25, "max": -5}, "SAR 2022")
m.center_object(country, 5)
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.
col_2021 = s1_col_2021.map(lambda img: img.focal_median(100, "circle", "meters"))
col_2022 = s1_col_2022.map(lambda img: img.focal_median(100, "circle", "meters"))
m = geemap.Map()
m.add_basemap("HYBRID")
sar_2021 = col_2021.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
sar_2022 = col_2022.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
m.add_layer(sar_2021, {"min": -25, "max": -5}, "SAR 2021")
m.add_layer(sar_2022, {"min": -25, "max": -5}, "SAR 2022")
m.center_object(country, 5)
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()
m.setCenter(68.4338, 26.4213, 7)
left_layer = geemap.ee_tile_layer(sar_2021, {"min": -25, "max": -5}, "SAR 2021")
right_layer = geemap.ee_tile_layer(sar_2022, {"min": -25, "max": -5}, "SAR 2022")
m.split_map(
left_layer, right_layer, left_label="Sentinel-1 2021", right_label="Sentinel-1 2022"
)
m.add_layer(country.style(**style), {}, country_name)
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_2021 = sar_2021.lt(threshold)
water_2022 = sar_2022.lt(threshold)
Create a split-view map to compare the pre-flood and flood water extent side by side.
m = geemap.Map()
m.setCenter(68.4338, 26.4213, 7)
m.add_layer(sar_2021, {"min": -25, "max": -5}, "SAR 2021")
m.add_layer(sar_2022, {"min": -25, "max": -5}, "SAR 2022")
left_layer = geemap.ee_tile_layer(
water_2021.selfMask(), {"palette": "blue"}, "Water 2021"
)
right_layer = geemap.ee_tile_layer(
water_2022.selfMask(), {"palette": "red"}, "Water 2022"
)
m.split_map(left_layer, right_layer, left_label="Water 2021", right_label="Water 2022")
m.add_layer(country.style(**style), {}, country_name)
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_2022.subtract(water_2021).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.setCenter(68.4338, 26.4213, 7)
m.add_layer(sar_2021, {"min": -25, "max": -5}, "SAR 2021")
m.add_layer(sar_2022, {"min": -25, "max": -5}, "SAR 2022")
left_layer = geemap.ee_tile_layer(
water_2021.selfMask(), {"palette": "blue"}, "Water 2021"
)
right_layer = geemap.ee_tile_layer(
water_2022.selfMask(), {"palette": "red"}, "Water 2022"
)
m.split_map(left_layer, right_layer, left_label="Water 2021", right_label="Water 2022")
m.add_layer(flood_extent, {"palette": "cyan"}, "Flood Extent")
m.add_layer(country.style(**style), {}, country_name)
m
Calculate SAR flood area¶
area_2021 = geemap.zonal_stats(
water_2021.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2021)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 68949.458824 | S Asia |
area_2022 = geemap.zonal_stats(
water_2022.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(area_2022)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 59224.121569 | S Asia |
flood_area = geemap.zonal_stats(
flood_extent.selfMask(), country, scale=1000, stat_type="SUM", return_fc=True
)
geemap.ee_to_df(flood_area)
Computing statistics ...
| abbreviati | country_co | country_na | sum | wld_rgn | |
|---|---|---|---|---|---|
| 0 | Pak. | PK | Pakistan | 12263.835294 | S Asia |
The total area of the flood extent is 12,263 square kilometers based on Sentinel-1 SAR images.