{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "[![Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/giswqs/notebook-share/blob/master/docs/brazil_floods.ipynb)\n", "[![Open in SageMaker Studio Lab](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/giswqs/notebook-share/blob/master/docs/brazil_floods.ipynb)\n", "[![Open in Planetary Computer](https://img.shields.io/badge/Open-Planetary%20Computer-black?style=flat&logo=microsoft)](https://pccompute.westeurope.cloudapp.azure.com/compute/hub/user-redirect/git-pull?repo=https://github.com/giswqs/notebook-share&urlpath=lab/tree/notebook-share/docs/brazil_floods.ipynb&branch=master)\n", "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/giswqs/notebook-share/HEAD?labpath=docs%brazil_floods.ipynb)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Visualization and Analysis of Brazil Floods 2024\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "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](https://en.wikipedia.org/wiki/2024_Rio_Grande_do_Sul_floods) page for more information about the 2024 Brazil floods.\n", "\n", "\n", "## Requirements\n", "\n", "To follow this tutorial, you must first [sign up](https://earthengine.google.com/signup) for a [Google Earth Engine](https://earthengine.google.com/) account. Earth Engine is a cloud computing platform with a [multi-petabyte catalog](https://developers.google.com/earth-engine/datasets/) of satellite imagery and geospatial datasets. It is free for noncommercial use. To authenticate the Earth Engine Python API, see instructions [here](https://book.geemap.org/chapters/01_introduction.html#earth-engine-authentication).\n", "\n", "In this tutorial, we will use the [geemap](https://geemap.org) 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.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Installation\n", "\n", "Uncomment the following line to install geemap if needed." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# !pip install geemap" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Import libraries\n", "\n", "Import the earthengine-api and geemap." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import ee\n", "import geemap.foliumap as geemap" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "geemap.ee_initialize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Download administrative boundaries\n", "\n", "Download the administrative boundaries of Rio Grande doSul, Brazil from [here](https://github.com/opengeos/datasets/releases/tag/places)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "url = \"https://github.com/opengeos/datasets/releases/download/places/Rio_Grande_do_Sul_Brazil.geojson\"" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Downloading...\n", "From: https://github.com/opengeos/datasets/releases/download/places/Rio_Grande_do_Sul_Brazil.geojson\n", "To: /tmp/18896865-507c-4806-bc74-6738b5fd4604.geojson\n", "100%|███████████████████████████████████████████████████████| 61.2k/61.2k [00:00<00:00, 2.24MB/s]\n" ] } ], "source": [ "roi = geemap.geojson_to_ee(url)\n", "# roi" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Longitude: -53.2453436538024, Latitude: -29.778651228599525\n" ] } ], "source": [ "centroid = roi.geometry().centroid(1)\n", "lon, lat = centroid.getInfo()[\"coordinates\"]\n", "print(f\"Longitude: {lon}, Latitude: {lat}\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Create an interactive map\n", "\n", "Specify the center point `[lat, lon]` and zoom level of the map." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "style = {\"fillColor\": \"00000000\", \"color\": \"FF0000\"}\n", "m.add_layer(roi.style(**style), {}, \"ROI\")\n", "geom = roi.geometry()\n", "m.center_object(geom, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "place_name = \"Rio Grande do Sul\"\n", "pre_flood_start_date = \"2024-01-01\"\n", "pre_flood_end_date = \"2024-04-27\"\n", "post_flood_start_date = \"2024-04-29\"\n", "post_flood_end_date = \"2024-05-17\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Create Landsat composites\n", "\n", "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](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The number of images in the pre-flood collection: 374\n" ] } ], "source": [ "pre_flood_col = (\n", " ee.ImageCollection(\"NASA/HLS/HLSL30/v002\")\n", " .filterBounds(roi)\n", " .filterDate(pre_flood_start_date, pre_flood_end_date)\n", " .filter(ee.Filter.lt(\"CLOUD_COVERAGE\", 15))\n", ")\n", "print(\n", " f\"The number of images in the pre-flood collection: {pre_flood_col.size().getInfo()}\"\n", ")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The number of images in the post-flood collection: 34\n" ] } ], "source": [ "post_flood_col = (\n", " ee.ImageCollection(\"NASA/HLS/HLSL30/v002\")\n", " .filterBounds(roi)\n", " .filterDate(post_flood_start_date, post_flood_end_date)\n", " .filter(ee.Filter.lt(\"CLOUD_COVERAGE\", 50))\n", ")\n", "print(\n", " f\"The number of images in the post-flood collection: {post_flood_col.size().getInfo()}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the Landsat 8 composite for the pre-flood and flood periods." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "\n", "pre_flood_image = pre_flood_col.median().clipToCollection(roi)\n", "post_flood_image = post_flood_col.median().clipToCollection(roi)\n", "\n", "vis_params = {\"bands\": [\"B6\", \"B5\", \"B4\"], \"min\": 0, \"max\": 0.4}\n", "m.add_layer(pre_flood_image, vis_params, \"Landsat Pre-flood\")\n", "m.add_layer(post_flood_image, vis_params, \"Landsat Post-flood\")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Compare Landsat composites side by side\n", "\n", "Compare the pre-flood and flood composites side by side." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "left_layer = geemap.ee_tile_layer(pre_flood_image, vis_params, \"Landsat Pre-flood\")\n", "right_layer = geemap.ee_tile_layer(post_flood_image, vis_params, \"Landsat Post-flood\")\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Landsat Pre-flood\",\n", " right_label=\"Landsat Post-flood\",\n", ")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Normalized Difference Water Index (NDWI)\n", "\n", "The [Normalized Difference Water Index](https://en.wikipedia.org/wiki/Normalized_difference_water_index) (NDWI) is a commonly used index for detecting water bodies. It is calculated as follows:\n", "\n", "$$NDWI = \\frac{Green - NIR}{Green + NIR}$$\n", "\n", "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.\n", "\n", "Landsat 8 imagery has [11 spectral bands](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1#bands). The Landsat 8 NDWI is calculated using the green (`B3`) and NIR (`B5`) bands.\n", "\n", "![](https://i.imgur.com/yuZthc6.png)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ndwi_pre = pre_flood_image.normalizedDifference([\"B3\", \"B5\"]).rename(\"NDWI\")\n", "ndwi_post = post_flood_image.normalizedDifference([\"B3\", \"B5\"]).rename(\"NDWI\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Compute the NDWI layers for the pre-flood and flood periods side by side." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "ndwi_vis = {\"min\": -1, \"max\": 1, \"palette\": \"ndwi\"}\n", "left_layer = geemap.ee_tile_layer(ndwi_pre, ndwi_vis, \"NDWI Pre-flood\")\n", "right_layer = geemap.ee_tile_layer(ndwi_post, ndwi_vis, \"NDWI Post-flood\")\n", "m.split_map(\n", " left_layer, right_layer, left_label=\"NDWI Pre-flood\", right_label=\"NDWI Post-flood\"\n", ")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Extract Landsat water extent\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "threshold = 0.1\n", "water_pre = ndwi_pre.gt(threshold)\n", "water_post = ndwi_post.gt(threshold)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Combine the pre-flood and surface water extent side by side." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "\n", "m.add_layer(pre_flood_image, vis_params, \"Landsat Pre-flood\", True)\n", "m.add_layer(post_flood_image, vis_params, \"Landsat Post-flood\", True)\n", "\n", "left_layer = geemap.ee_tile_layer(\n", " water_pre.selfMask(), {\"palette\": \"blue\"}, \"Water Pre-flood\"\n", ")\n", "right_layer = geemap.ee_tile_layer(\n", " water_post.selfMask(), {\"palette\": \"red\"}, \"Water Post-flood\"\n", ")\n", "\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Water Pre-flood\",\n", " right_label=\"Water Post-flood\",\n", ")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Extract Landsat flood extent\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "flood_extent = water_post.subtract(water_pre).gt(0).selfMask()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Add the flood extent layer to the map." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "\n", "m.add_layer(pre_flood_image, vis_params, \"Landsat Pre-flood\", True)\n", "m.add_layer(post_flood_image, vis_params, \"Landsat Post-flood\", True)\n", "\n", "left_layer = geemap.ee_tile_layer(\n", " water_pre.selfMask(), {\"palette\": \"blue\"}, \"Water Pre-flood\"\n", ")\n", "right_layer = geemap.ee_tile_layer(\n", " water_post.selfMask(), {\"palette\": \"red\"}, \"Water Post-flood\"\n", ")\n", "\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Water Pre-flood\",\n", " right_label=\"Water Post-flood\",\n", ")\n", "\n", "m.add_layer(flood_extent, {\"palette\": \"cyan\"}, \"Flood Extent\")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate Landsat flood area\n", "\n", "To calculate the flood area, we can use the [`geemap.zonal_stats()`](https://geemap.org/common/#geemap.common.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." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.5857116801.019608
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 16801.019608 " ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "area_pre_flood = geemap.zonal_stats(\n", " water_pre.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(area_pre_flood)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.585711884.207843
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 1884.207843 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "area_2022 = geemap.zonal_stats(\n", " water_post.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(area_2022)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.58571634.188235
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 634.188235 " ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "flood_area = geemap.zonal_stats(\n", " flood_extent.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(flood_area)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Create Sentinel-1 SAR composites\n", "\n", "Besides Landsat, we can also use Sentinel-1 [Synthetic Aperture Radar (SAR)](https://www.earthdata.nasa.gov/learn/backgrounders/what-is-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. \n", "\n", "Sentinel-1 operates in four exclusive [acquisition modes](https://sentinels.copernicus.eu/web/sentinel/missions/sentinel-1/instrument-payload):\n", "\n", "* Stripmap (SM)\n", "* Interferometric Wide swath (IW)\n", "* Extra-Wide swath (EW)\n", "* Wave mode (WV)\n", "\n", "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.\n", "\n", "The [Sentinel-1 SAR data](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD) are available from 2014 to present. Let's filter the `COPERNICUS/S1_GRD` dataset by the date range and location." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pre_flood_start_date = \"2023-10-01\"\n", "pre_flood_end_date = \"2024-04-27\"\n", "post_flood_start_date = \"2024-04-29\"\n", "post_flood_end_date = \"2024-05-17\"" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The number of images in the pre-flood collection: 161\n" ] } ], "source": [ "s1_col_pre = (\n", " ee.ImageCollection(\"COPERNICUS/S1_GRD\")\n", " .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VV\"))\n", " .filter(ee.Filter.eq(\"instrumentMode\", \"IW\"))\n", " .filterDate(pre_flood_start_date, pre_flood_end_date)\n", " .filterBounds(roi)\n", " .select(\"VV\")\n", ")\n", "print(\n", " f\"The number of images in the pre-flood collection: {s1_col_pre.size().getInfo()}\"\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Create the Sentinel-1 image collection for the flood period." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The number of images in the post-flood collection: 19\n" ] } ], "source": [ "s1_col_post = (\n", " ee.ImageCollection(\"COPERNICUS/S1_GRD\")\n", " .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VV\"))\n", " .filter(ee.Filter.eq(\"instrumentMode\", \"IW\"))\n", " # .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))\n", " .filterDate(post_flood_start_date, post_flood_end_date)\n", " .filterBounds(roi)\n", " .select(\"VV\")\n", ")\n", "print(\n", " f\"The number of images in the post-flood collection: {s1_col_post.size().getInfo()}\"\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Create Sentinel-1 SAR composites for the pre-flood and flood periods." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "m.add_basemap(\"SATELLITE\")\n", "sar_pre = s1_col_pre.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)\n", "sar_post = s1_col_post.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)\n", "m.add_layer(sar_pre, {\"min\": -25, \"max\": -5}, \"SAR Pre-flood\")\n", "m.add_layer(sar_post, {\"min\": -25, \"max\": -5}, \"SAR Post-flood\")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Apply speckle filtering\n", "\n", "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](https://doi.org/10.1080/02757259409532206)). 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." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "s1_col_pre = s1_col_pre.map(lambda img: img.focal_median(100, \"circle\", \"meters\"))\n", "s1_col_post = s1_col_post.map(lambda img: img.focal_median(100, \"circle\", \"meters\"))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "m.add_basemap(\"SATELLITE\")\n", "sar_pre = s1_col_pre.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)\n", "sar_post = s1_col_post.reduce(ee.Reducer.percentile([20])).clipToCollection(roi)\n", "m.add_layer(sar_pre, {\"min\": -25, \"max\": -5}, \"SAR Pre-flood\")\n", "m.add_layer(sar_post, {\"min\": -25, \"max\": -5}, \"SAR Post-flood\")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Compare Sentinel-1 SAR composites side by side\n", "\n", "Create a split-view map to compare the pre-flood and flood SAR composites side by side." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "left_layer = geemap.ee_tile_layer(sar_pre, {\"min\": -25, \"max\": -5}, \"SAR Pre-flood\")\n", "right_layer = geemap.ee_tile_layer(sar_post, {\"min\": -25, \"max\": -5}, \"SAR Post-flood\")\n", "\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Sentinel-1 Pre-flood\",\n", " right_label=\"Sentinel-1 Post_flood\",\n", ")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Extract SAR water extent\n", "\n", "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](https://doi.org/10.1016/j.isprsjprs.2019.10.017)). 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 \n", "dB." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "threshold = -18\n", "water_pre = sar_pre.lt(threshold)\n", "water_post = sar_post.lt(threshold)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Create a split-view map to compare the pre-flood and flood water extent side by side." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "m.add_layer(sar_pre, {\"min\": -25, \"max\": -5}, \"SAR Pre-flood\")\n", "m.add_layer(sar_post, {\"min\": -25, \"max\": -5}, \"SAR Post-flood\")\n", "\n", "left_layer = geemap.ee_tile_layer(\n", " water_pre.selfMask(), {\"palette\": \"blue\"}, \"Water Pre-flood\"\n", ")\n", "right_layer = geemap.ee_tile_layer(\n", " water_post.selfMask(), {\"palette\": \"red\"}, \"Water Post-flood\"\n", ")\n", "\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Water Pre-flood\",\n", " right_label=\"Water Post-flood\",\n", ")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Extract SAR flood extent\n", "\n", "Similar to the Landsat approach, we can subtract the pre-flood water extent from the flood water extent to extract the flood extent." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "flood_extent = water_post.subtract(water_pre).gt(0).selfMask()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = geemap.Map()\n", "\n", "m.add_layer(sar_pre, {\"min\": -25, \"max\": -5}, \"SAR Pre-flood\")\n", "m.add_layer(sar_post, {\"min\": -25, \"max\": -5}, \"SAR Post-flood\")\n", "\n", "left_layer = geemap.ee_tile_layer(\n", " water_pre.selfMask(), {\"palette\": \"blue\"}, \"Water Pre-flood\"\n", ")\n", "right_layer = geemap.ee_tile_layer(\n", " water_post.selfMask(), {\"palette\": \"red\"}, \"Water Post-flood\"\n", ")\n", "\n", "m.split_map(\n", " left_layer,\n", " right_layer,\n", " left_label=\"Water Pre-flood\",\n", " right_label=\"Water Post-flood\",\n", ")\n", "\n", "m.add_layer(flood_extent, {\"palette\": \"cyan\"}, \"Flood Extent\")\n", "m.add_layer(roi.style(**style), {}, place_name)\n", "m.center_object(roi, 6)\n", "m" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate SAR flood area" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.5857110160.776471
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 10160.776471 " ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "area_pre_flood = geemap.zonal_stats(\n", " water_pre.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(area_pre_flood)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.5857113054.921569
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 13054.921569 " ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "area_2022 = geemap.zonal_stats(\n", " water_post.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(area_2022)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Computing statistics ...\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
area_1area_tot_gcodigogeocodigomapidmslinknomeperimetro_sum
0281748.155281748.155434399223RIO GRANDE DO SUL3524.585715745.360784
\n", "
" ], "text/plain": [ " area_1 area_tot_g codigo geocodigo mapid mslink nome \\\n", "0 281748.155 281748.155 43 43 99 223 RIO GRANDE DO SUL \n", "\n", " perimetro_ sum \n", "0 3524.58571 5745.360784 " ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "flood_area = geemap.zonal_stats(\n", " flood_extent.selfMask(), roi, scale=1000, stat_type=\"SUM\", return_fc=True\n", ")\n", "geemap.ee_to_df(flood_area)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.8" }, "vscode": { "interpreter": { "hash": "31f05ea183a4718249d13ada7f166c6bdba1d00716247af5c11c23af8d5923f1" } } }, "nbformat": 4, "nbformat_minor": 4 }