ASF Logo
North American Space Agency (NASA) Logo

Alaska Satellite Facility - Distributed Active Archive Center

How to Map Regional Inundation with Spaceborne L-band SAR using QGIS – NEEDS TESTING

A frame of an animation of inundation in the Amazon rain forest. Credit: Chapman, Guritz 2016; RTC: ASF 2015; Includes Material © JAXA/METI 2007.

Source: Bruce Chapman, JPL, and Rick Guritz, ASF. Adapted to QGIS and GIMP by Eric Anderson, ESSC, English translation Jagan Nautiyal

◆ Advanced (OSX): Use QGIS and GIMP to classify and map regional inundation.

Background

The all-weather ability of synthetic aperture radar (SAR) to penetrate cloud cover and low-light conditions to acquire imagery of the Earth’s surface is well known. Given the availability of high-resolution terrain models at 30 meters from the Shuttle Radar Topography Mission (SRTM) and an interest to make SAR data easier to use, the Alaska Satellite Facility (ASF) offers radiometric terrain corrected (RTC) L-band image products from ALOS PALSAR data.

At longer wavelengths such as L-band, SAR can penetrate surface vegetation including the Amazon rain forest, and flooding events can be mapped very accurately at high resolution. With the availability of RTC data, it is now relatively easy to do accurate flood mapping using PALSAR data. In this recipe, we describe methods used to (a) create an inundation image animation from 18 PALSAR RTC products and (b) produce an inundation map that quantifies the number of flooding events from the given set of data.

Prerequisites

Materials List for Generating an Animation:

  • Windows, Mac OS X
  • ALOS-1 PALSAR RTC products
  • QGIS (version 2.18.12 used in this recipe) with Orfeo Toolbox (Windows users see How to Access Orfeo Toolbox (Windows))
  • GIMP (version 2.8 used in this recipe)

Materials List for Generating a Colorized Map:

Sample Granule List

Click on the product names to download each granule individually.

Granule Name Product Name
ALPSRP048877090 AP_04887_FBS_F7090_RT1.zip
ALPSRP075717090 AP_07571_FBD_F7090_RT1.zip
ALPSRP082427090 AP_08242_FBD_F7090_RT1.zip
ALPSRP102557090 AP_10255_FBS_F7090_RT1.zip
ALPSRP115977090 AP_11597_FBS_F7090_RT1.zip
ALPSRP122687090 AP_12268_FBD_F7090_RT1.zip
ALPSRP129397090 AP_12939_FBD_F7090_RT1.zip
ALPSRP136107090 AP_13610_FBD_F7090_RT1.zip
ALPSRP156237090 AP_15623_FBS_F7090_RT1.zip
ALPSRP183077090 AP_18307_FBD_F7090_RT1.zip
ALPSRP189787090 AP_18978_FBD_F7090_RT1.zip
ALPSRP209917090 AP_20991_FBS_F7090_RT1.zip
ALPSRP223337090 AP_22333_FBS_F7090_RT1.zip
ALPSRP230047090 AP_23004_FBD_F7090_RT1.zip
ALPSRP243467090 AP_24346_FBD_F7090_RT1.zip
ALPSRP263597090 AP_26359_FBS_F7090_RT1.zip
ALPSRP270307090 AP_27030_FBS_F7090_RT1.zip
ALPSRP236757090 AP_23675_FBD_F7090_RT1.zip

Part 1: Generating an Inundation Animation

Download and Organize the Data

    1. Download the granules from the Prerequisites section.
    2. Extract data from the downloaded .zip files and copy the HH image from each granule into a new folder (e.g., HH).

Convert Images to Decibels in QGIS

    1. Create a new directory to place your output images (e.g., db_HH)
    2. Open QGIS and use the Open Raster button to load the HH images from the folder you created in the Download and Organize the Data step.
    3. Ensure that your HH images are open in the Layers Panel (there should be 18 of them)
GDAL Raster Calculator with batch processing
    1. Ensure that your Processing Toolbox is open and ready (Processing menu > Toolbox)
    2. Within the Processing Toolbox, search for “raster calculator”
    3. Right-click on Raster calculator (GDAL) and execute as batch process.
    4. In the Parameters tab, click on “…” next to the first box within Input Layer A
    5. Select from open layers and click Select All
    6. Click OK to add all image layers to Input Layer A
    7. Drag the navigation bar over until you see “Calculation in gdalnumeric syntax…”
    8. In the first box, enter “10*log10(A)”
    9. Double-click on the title above where you just entered your formula
    10. To the right of the Calculation header, find “Output raster type“
    11. Under Output raster type, change the selection from “Byte” to “Float32”
    12. Double-click on Output raster type to apply to all
    13. To the right of Output raster type, find Calculated
    14. Navigate to the directory/folder you created and enter “db_” and click Save
    15. When prompted, select Fill with parameter values and click OK
    16. To the right of Calculated find “Load in QGIS”
    17. Change the first “Yes” to “No” and double-click on Load in QGIS
    18. Click Run to begin batch processing
    19. Close Batch Processing dialog window once complete
    20. Close all images in QGIS using the Layers Panel close button.

Apply Despeckling to dB Converted Images

    1. Create a new directory to place your output images (e.g., sf_db_HH)
    2. In QGIS, use the Open Raster button to load the decibel-converted images you created in the Convert Images to Decibels in QGIS step.
    3. Within the Processing Toolbox, search for “Despeckle”
    4. Right-click on Despeckle (lee) (Orfeo) and execute as batch process.
Despeckle (lee) tool with batch processing
    1. Under “Input Image,” click the “…” and select from open Layers
    2. Click “Select All” and click OK to add all images
    3. Scroll right until you find “Output Image” and select “…“
    4. Navigate to the directory/folder you created and enter “sf_” and click Save
    5. When prompted, select “Fill with parameter values” and click OK
    6. To the right of “Fill with parameters values” find “Load in QGIS”
    7. Select “No” and double-click on “Load in QGIS” to apply to all
    8. Click Run to begin batch processing
    9. Close Batch Processing dialog window once complete
    10. Close all images in QGIS using the Layers Panel close button

Exclude Erroneous Pixels at the Edge of Your Images

The Despeckle tool leaves a strip of extra pixels with very negative values to the edge of the images, which must be removed before the following analysis. When you open the despeckle images, you will find that the images are not correctly rendering due to these negative values.

    1. Create a new directory to place your output images (e.g., xsf_db_HH)
    2. In QGIS, use the Open Raster button to load the images created in the Apply Despeckling to dB Converted Images step.
    3. Within the Processing Toolbox, search for “raster calculator”
    4. Right-click on Raster calculator (GDAL) and execute as batch process
    1. In the Parameters tab, click on “…” next to the first box within Input Layer A
    2. Select from open layers and click Select All
    3. Click OK to add all image layers to Input Layer A
    4. Drag the window bar over until you see “Calculation in gdalnumeric syntax…”
    5. In the first box, enter “(A>-9999)*A”
    6. Double-click on Calculation in gdalnumeric syntax… to apply to all
    7. To the right of Calculation, find “Set output nodata value”
Raster Calculator with correct setup
    1. Enter “0” in the first box and double-click on Set output nodata value
    2. To the right of the Calculation, find “Output raster type”
    3. Under Output raster type, change the selection from “Byte” to “Float32.”
    4. Double-click on Output raster type to apply to all
    5. To the right of Output raster type, find “Calculated” 
    6. Click on “…” next to the first box within Calculated
    7. Navigate to the directory/folder you created and enter “x” and click Save
    8. When prompted, select Fill with parameter values and click OK
    9. To the right of Calculated find “load in QGIS”
    10. Change the first “Yes” to “No” and double-click on Load in QGIS
    11. Click Run to begin batch processing
    12. Close the Batch Processing dialog window once complete
    13. Close all images in QGIS using the Layers Panel close button.

Classify the Images using a Threshold

This step assumes that your images have a bimodal distribution of pixel values in the image. If you are using different data, the class values will vary and it may take several iterations to determine the appropriate thresholds. For the process of determining the appropriate thresholds, it is good practice to consider one image representing the maximum flooding and another image representing the minimum flooding. Knowledge of the geographical coverage of seasonal flooding helps confirm whether the results are correct.

    1. In QGIS, use the Open Raster button to load the images created in the Exclude Erroneous Pixels at the Edge of Your Images step.
    2. Right-click on the first image in the Layers Panel and choose Properties
    3. In the Style menu, change the Render type to Singleband pseudocolor
    4. Click on the Histogram tab from the Layer Properties window you have open.
    5. View the histogram to interpret the distribution of pixel values in the image

Note: You should see a clearly bi-modal histogram with water pixels appearing significantly darker than the main image data.

    1. To separate water from the rest of the image, pick a threshold at (or near) the minimum between the two modes of the distribution (e.g., -13)
    2. Return to the Style tab of the Layer Properties window and change the Interpolation to Discrete
    3. Next to Classify, click on the green cross button three times to add your class value.
Histogram showing bimodal distribution of pixel values.
Classification of the image with determined thresholds
    1. Change the value of the classes to the following: “-13; “-2”; “inf”
    2. Double-click on the Color to change the class colors
    3. Set “-13” to blue to represent the main channels of water. Refine the breakpoint values until water is mostly contained within river channels and a very small number of blue dots show on land outside of river channels.
    4. Set “-2” to black to represent land
    5. Set “inf” to yellow to represent brighter areas because of double-bounce, which indicates flooded vegetation. To increase the prominence of double-bounce, you must adjust “-2” to a lower value, which adjusts the breakpoint of your second threshold value.
    6. Adjust the following values until you are satisfied. The class values will be determined by the pixel values in the image.
    7. Optional — You may adjust Resampling in the Style menu to smooth rendered pixels
    8. Apply the same Style properties to all images by right-clicking on the image with the desired classification and select Styles > Copy Style
    9. Click the up button in the Layers Panel to collapse all layers and select all other layers (using Shift)
    10. Right-click on any of the selected images and click on Paste Style

Export Classified Images as PNG Image Files

    1. Create a new directory to save your exported images (e.g., Flood_PNG)
    2. In QGIS, ensure that all of your images are still open in the Layers Panel
    3. After applying the correct thresholds and colors to all images, turn off visibility for all layers except the first one (xsf_db_AP_27030_FBS_F7090_RT1_HH) by clicking on the checkbox next to the raster layer’s name
    4. Use the magnification button in the Map Navigation toolbar to center your image to prepare for PNG exporting.
All of the images with the same symbology
    1. Navigate to Project > Save as Image

Tip: These steps can be cumbersome if you manually type in each file name. To speed up the process, right-click on each image and select “Rename.” Once highlighted, you can copy the name of the file before saving your image and simply paste the filename in the box, then add a “.png” file extension.

    1. Save as xsf_db_AP_27030_FBS_F7090_RT1_HH.png and click OK
    2. Turn off the first image and select the next image
    3. Repeat steps 1-5 for all images until you have 18 PNG image files

Generate and Display Animated GIF with GIMP

    1. Select or create a directory to save your animated GIF file
    2. Open GIMP 2 and select File > Open as Layers…
    3. Navigate to the directory where your image files are located and select all
    4. Once all images are loaded, select File > Export As
    5. Select the directory you chose to save your animation and save the animation as Flood_Animation.gif and an Export GIF dialog window will appear.
GIMP Export Image as GIF dialog
    1. Check the As animation box and set Delay between frames to 500
    2. Check Use delay entered above for all frames and click Export
    3. Navigate to the Flood_Animation.gif file and open the animation in a browser by using right-click Open as…

Sample Inundation Animation

Part 2: Generating a Colorized Inundation Map

Classify the Images Using a Threshold in QGIS

Classify each image to separate water (1) from land (0).

    1. Create a new directory to save your output files (e.g., Reclass _HH)
    2. In QGIS, use the Open Raster button to load the images created in Part 1 — Exclude Erroneous Pixels at the Edge of Your Images (e.g., xsf_db_AP_10255_FBS_F7090_RT1_HH)
    3. Within the Processing Toolbox, search for “reclassify”
    4. Right-click on Reclassify values (simple) (SAGA) and run as batch process.
    5. Click on “…” next to the first box within Grid
    6. Select from open layers and click Select All
    7. Click OK to add all image layers to Grid
    8. To the right of Grid, find “Replace Condition”
    9. Change the condition to (1) Low value < grid value < high value
    10. Double-click on Replace Condition to apply to all
    11. To the right of Replace Condition find lookup Table
Reclassify values (simple) with batch processing
Table to reclassify the values
    1. Click on “…” to the right of the first box within Lookup Table
    2. In the Fixed table window, assign the following values seen in the image to the right.
    3. Repeat for all images (batch process won’t allow applying to all)
    4. To the right of Lookup Table, find Changed Grid
    5. Click on “…” next to the first box within Changed Grid
    6. Navigate to the directory you created and enter “r_” and click Save
    7. When prompted, fill with parameter values and click OK
    8. To the right of Changed Grid, find “Load in QGIS”
    9. Change the first option to “No” and double-click “Load in QGIS” to apply to all
    10. Click Run to begin batch processing
    11. Close Batch Processing dialog window once complete
    12. Close all images in QGIS using the Layers Panel close button

Create the Final Product

    1. Choose a directory to save your colorized inundation map
    2. In QGIS, use the Open Raster button to load your reclassified images created in Part 2 — Classify the Images Using a Threshold in QGIS
    3. Open the Raster Calculator by navigating to Raster > Raster Calculator
    4. Double-click on the first image and click on the “+” operator
    5. Continue adding each image and symbol until all images are in the equation
  1.  
Raster Calculator with correct setup
    1. Set the Output Layer to a name and directory of your choosing (e.g., count.tif)
    2. Click OK to run the Raster Calculator

Note: The combined image will display as a range of gray values. To change this, you may classify the pixel values of the image for analysis. It is good to assign values 0 and 1 to black because these values are usually a result of speckling that remained after filtering. It helps to choose colors that are distinguishable and have logic in terms of the combined image (higher values indicate more frequent flooding).

    1. Right-click on the image and click on Properties
    2. In the Style tab, change render type to singleband pseudocolor
    3. Under Color options, select a color ramp from the list or choose random colors
    4. Change the classification method to Equal Interval
    5. Change the number of classes to “19”
    6. Optional: adjust resampling type to smooth the rendered image
    7. When you are satisfied, click OK to see the results

Sample Colorized Inundation Maps

Map of inundation in Amazon rain forest. The color key provides the color for the number of times an area was flooded, and 18 represents constant river. Numbers less than 18 are the number of times an area floods. Credit: Anderson 2017; RTC: ASF 2015; Includes Material © JAXA/METI 2007

How to Access Orfeo Toolbox (Windows)

To use the despeckle tool, you must have access to the Orfeo Toolbox, which is not included with the Windows standalone version of QGIS 2.18.12. To add the toolbox, you must download it using the OSGeo4W installer.

Note: For these steps to work properly, you must have QGIS pre-installed.

    1. Navigate to QGIS.org
    2. Select Download Now
    3. Under Advanced Users, download OSGeo4W setup
    4. Open OSGeo4W setup
    5. From the dialog window, select Advanced Install and click Next
    6. Select Install from Internet and click Next
    7. Select Root Install Directory (e.g., C:\OSGeo4W) and click Next 
    8. Select Local Package Directory (default works fine) and click Next
    9. Select Direct Connection and click Next
    10. Click on https://download.osgeo.org to select and click Next
    11. In Select Packages menu, type in “otb” in the search box
    12. Click on Libs to expand
    13. Under Libs, click on each Skip once and click Next
    14. Leave all defaults and click Next
    15. Click on the checkbox and click Next
    16. Your download will begin
    17. Click Finish once installation completes
    18. Open QGIS 2.18.12
    19. Select Processing from the main menu bar
OSGeo4W Setup
Processing Options: Providers
    1. Click on Options…
    2. Click on Providers to expand
    3. Click on GRASS commands to expand
    4. Double-click in blank space to the right of Msys folder
    5. Delete the existing path for “Msys folder”
    6. Click on Orfeo Toolbox to expand
    7. Ensure Activate is checked
    8. Double-click in blank space to the right of OTB applications folder
    9. Enter the path of your OTB applications folder in the box (e.g., “C:\OSGeo4W64\bin”)
    10. Double-click in blank space to the right of OTB command line tools folder
    11. Enter the path of the OSGeo4W64 bin folder in the box (e.g., “C:\OSGeo4W64\bin”)
    12. Exit QGIS and restart the application
    13. In the Toolbox menu on the right, the Orfeo Toolbox will now appear
    14. Return to Prerequisites