Lab6 - SAR Application for Studying Pakistan Floods 2022 Using Earth Engine and Geemap

Colab Open in SageMaker Studio Lab Open in Planetary Computer Binder

Lab6 - SAR Application for Studying Pakistan Floods 2022 Using Earth Engine and Geemap#

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

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

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

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 Sentinel-1 SAR composites#

Sentinel-1 Synthetic Aperture Radar (SAR) data can be used 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.

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")
)

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")
)

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.