Flood Mapping with Sentinel-1 Data using SNAP and QGIS

I haven't done remote sensing work in a long while and have been wanting to get back into. Remote sensing is a great complement to GIS, and I've been looking for a case study to do a short tutorial. I thought the recent flooding on Kauai Island on April 14-15, 2018 would be a good case to do some flood mapping using remote sensing techniques.

Background:

The Island of Kauai is not called the Garden Isle for no reason - it is one of the rainiest place on Earth. On April 14-15, 2018 unprecedented rains caused major flooding and a series of landslides on Kauai. The National Weather Service recorded 27.52 inches (685 mm) of rainfall in the Hanalei area during the 24-hour period from April 14 to the 15. After the flooding, the County of Kauai was declared a state of emergency.

More info of the flooding event can be found here: http://floodlist.com/america/usa/record-rainfall-floods-hawaii0april-2018

Workflow/Requirements:

Here is the work flow plan and requirements for this tutorial:

  1. Download Sentinel-1 data - require registration/sign up (free)
  2. Download and install SNAP software (if you don't already have it)
  3. Pre-processing - calibrate and speckle filtering of data
  4. Determine flooded areas - threshold binning (binarization, unsupervised classification)
  5. Post-processing - geometric correction of data (apply map projection)
  6. View resulting flood areas in GIS (e.g. QGIS, Google Earth, or ArcGIS  assuming you have a GIS software already)

Step 1: Download Sentinel-1 Imagery

Assessment of flooded areas can be done using radar data, such as the Sentinel-1 SAR imagery provided by the European Space Agency (ESA). For this example, I'm using Sentinel-1 Level-1 Ground Range Detected (GRD) data which already has some basic pre-processing done on it. More info about the Sentinel missions, products, and uses can be found on the Sentinel Online Website. Access to the Sentinel data is provided through the Copernicus Open Access Hub. Although the data is free to use, you will need to register or sign up first before you can download the data.

Go to Copernicus Open Access Hub:

  • Register or sign up for a free account then log in
  • On the map, navigate to the Hawaiian Islands and zoom in to the island of Kauai
  • Use the Polygon tool (in the bottom left, box tool doesn't seem to work) and draw an area of interest on Kauai (northern portion of the island)
  • Click the Advance Search menu (in top left) and use the following search parameters
    • Sensing Period: From = 4/15/2018, To = 4/15/2018
    • Mission Sentinel-1: Product Type = GRD
    • Note: there is only 1 image for Kauai taken on the date specified so there's no need to include other parameter criteria
  • Click the Search button when done
  • There should only be 1 result based on the search criteria for the area of interest. Note: when you click the link to download the file, you may be asked again to login (think this might be a time out issue on the Copernicus site).
Fig0001.png
Fig0002.png

Notice the available imagery is Level -1 Ground Range Detected (GRD) Sentinel-1 C-band (SAR-C) data collected in Interferometric Wide (IW) mode, allowing for an image acquisition over a 250km swath with resolution of 10m (5x20m).

Step 2: Download and Install SNAP

SNAP is a application developed by the ESA to work with Sentinel data. You will need to download and install the appropriate software for your computer. I am using the Windows 64bit version of the Sentinel Toolbox.

Fig0003.png

Step 3: Explore the Image and Create Subset

After you've install the SNAP application, you can use it to open and explore the Sentinel-1 image you downloaded previously. The SNAP application can open the zipped file without need to unzip it first.

  • In the SNAP application >> Open Product >> S1A_IW_GRDH_1SDV_20180415.zip (file downloaded earlier) >> Open
  • The Product Explorer Panel on the left shows relevant information/metadata related to the Sentinel data. You can right click on any of the file/item to get a context menu to view properties related to that item.
  • Under Bands, you'll see that for each polarization (VV + VH or HH + HV), there are 2 bands (Amplitude and Intensity). Intensity band is a virtual band and is just a square of the amplitude. To view an image of a Band, double click it and an image will be created for viewing in the left panel. Double click on the Intensity_VV band to view it - this is the band being used.

Note: Although any of the polarization/band can be used to analyze flooding extent, some studies have shown that VV gives a slight advantage when delineate flooding when using Sentinel‐1 data. This article, "Multi‐temporal synthetic aperture radar flood mapping using change detection" provides a good comparison using different polarizations for delineating flooding using Sentinel data.

Fig0004.png
  • To create a subset a portion of the island, first zoom in to the extent of the area you want to subset - which in this case is the northern portion of the island where Hanelei town is located.
Fig0005.png
  • Create the subset, go to Raster Menu >> Subset
  • In the Specify Product Subset window, the lat/long coordinates for the area of interest should already be filled in based on the zoomed in map extent. Click OK to create subset
Fig0006.png

The newly created subset image should show up in the Product Explorer panel and look something like this.

Fig0007.png

Step 4: Pre-Processing - Calibration and Speckle Filtering

Pre-processing calibration and speckle filtering here is done on the subset image you created in step 3 above. The pixel values in a SAR imagery can be related to the radar backscatter of the scene taken - what calibration does is transform the pixel values from the digital values recorded by the sensor into backscatter coefficient values -- essentially calibrated values of the backscatter coefficient. Speckle filtering is just reducing the noise in the image so to obtain higher quality imagery.

  • Calibration Tool: Go to Radar Menu >> Radiometric >> Calibrate
  • In the Calibration window:
    • I/O Parameters tab: source = your subset image from step 3; target product = your output file
    • Processing Parameters tab: Polarizations = VV; check Output sigma0 band >> Run
    • Close the Calibration window when it finished running
    • The resulting output file should be shown in the Product Explorer panel
  • The Calibration tool will create a new image with the backscatter coefficient
Fig0008.png
  • Speckle Filtering Tool: go to Radar Menu >> Speckle Filtering >> Single Product Speckle Filter
  • In the Single Product Speckle Filter window:
    • I/O Parameters tab: source = your calibrated image from the above calibration step; target product = your output file
    • Processing Parameters tab: Source Bands = Sigma0_VV; Filter = Lee; Filter size = 3x3 >> Run
    • Close the window when it finishes running
    • The resulting output file should then be shown in the Product Explorer panel
Fig0009.png

Step 5: Determining Flooded Areas (areas covered by water)

Using a histogram of the filtered backscatter coefficient and applying a threshold, the water pixels can be separated from the non-water pixels. I think of this as a form of unsupervised classification - you apply a threshold value and let the software decide what is considered water and what isn't. After this first pass you could further refine results using a supervised classification method (not covered in this tutorial).

  • To view the histogram: first make sure you're viewing the filtered sigma0_VV image. In the left side panel (below Product Explorer), select the Color Manipulation tab
  • In the Color Manipulation tab: a histogram is displayed showing one or more peaks. Low values of backscatter will correspond to water and high values to non-water class. You will need to use a threshold value to get the best representation for water.
    • You may need to select the logarithmic display for the histogram
    • Use the threshold slider to see different threshold values -  I am using the value 1.28E-2 for this example
Fig0010.png
  • After determining a good threshold value, the next step is to create a new image - basically a binary image of water (1 = true or 255 = white) and non-water (0 = false or 0 = black). To do this you would use a raster math calculation.
    • Go to Raster Menu >> Band Math..
    • In the Band Maths window: type in a name and uncheck the Virtual file box (if you want to save to disk)
    • Click Edit Expression: then in the window that comes up, type in 255*(sigma0_VV< 1.28E-2). This expression will return a 1 or true for pixel values less than the threshold and 0 or false for values higher than the threshold, this is then multiplied by 255 to get the binary black and white image
    • Click OK to run

Fig0011.png
  • A newly create Water image will be added to your product. An example result is shown below. Notice that since Kauai is surrounded by ocean so it makes it harder to determine flooded areas on land. I guess a quick work around solution would be just to clip the image to the coastline so ocean is seen as "flooded" areas.
Fig0012.png

Step 6: Post-processing - Geometric Correction

The Sentinel-1 image that was downloaded is in the geometry of the sensor, meaning it does not have geographic coordinates so you will need to applying a coordinate system if you want to view it in a GIS with other data.

  • Go to Radar Menu >> Geometric >> Terrain Correction >> Range-Doppler Terrain Correction
  • In the Range-Doppler Terrain Correction window:
    • I/O Parameters tab: source = your filtered image, then input other parameters if needed or use defaults
    • Processing Parameters tab: source bands = Water; DEM = use your own or the default SRTM 3Sec (Auto Download) which requires an internet connection; change other parameters if needed or use defaults
    • Click Run and Close window when completed

Note: If you save as BEAM-DIMAP format the output image will be an *.img format (see optional step 7 below)

Fig0013.png
  • A newly create image will be added to the Product Explorer panel. Double click on the band to open it for viewing. Here is a sample image below
Fig0014.png

Optional Step 7: Export Image for Viewing in GIS or Google Earth

At the point you should be able to use the output result in an GIS - the default output format (from step 6 above) is an *.img format. If you want to export it to other formats or view in Google Earth you can do that as well.

  • Go to File Menu >> Export >> Other >> View as Google Earth KMZ   OR
  • Go to File Menu >> Export >> Choose an export format
Fig0015.png

Here is an example resulting flood image as view in QGIS 3.0. In QGIS just add in your *.img or exported image and use the layer's properties to change color and transparency of the image. Flooded areas are shown in RED in the image.

Fig0016.png

Assessment of the resulting flooded areas

The Dartmouth Flood Observatory at the University of Colorado (previously at Dartmouth College) have mapped the maximum observed extent of the flooding on Kauai - see DFO Flood Event 4601 USA (Kauai, Hawaii). This can be used to assess how well your results are. The snapshots below show the comparison between the more accurate DFO Observed Flooding and the tutorial result. I would say that tutorial result is not bad given that no refinement was applied to the methodology other than using the SNAP geoprocessing defaults. A more accurate result could have been attained by using different thresholds, and doing some post-editing or processing to further refine the flooded areas such as supervised classification to weed out false positives (wrong classification of water).  Also keep in mind that no other supplement data were used except the Sentinel-1 data. Anyway, the point of this tutorial was just to get you thinking about remote sensing application.

Fig0017.png
Fig0018.png

That is it for this tutorial. I hope it got you thinking more about applying or using remote sensing in addition to GIS. As always thanks for reading.

Adding Basemaps in QGIS 3.0

This post is a very quick guide on adding basemaps in QGIS 3.0. There are probably lots of options but these are the 3 that I prefer to quickly add many types of basemaps for use in QGIS.

1. Quick Map Services Plug In

Using the QuickMapServices Plugin is probably the easiest way to add basemaps in QGIS. This beloved plugin has just been updated to work with the new QGIS 3.0. To add in the plugin just do the following:

  • Open QGIS. Go to Plugins Menu >> Manage and Install Plugins...
Fig0001.png
  • In the Plugins Window, search for QuickMapServices then click Install button. Note if you have this plugin installed from a previous version of QGIS then it may have been updated as well.
Fig0002.png
  • After it installs, you should be able to find QuickMapServices button in the Web Toolbar.
Fig0003.png
  • Click on the button to get the drop down menu >> Settings. In the QuickMapServices Settings window, go to the More Services tab >> Click Get Contributed Pack to get more basemaps. Then go to the Visibility tab and turn on/off the basemaps you want to display in the menu
Fig0005.png

2. XYZ Tiles

The second way to add basemaps in QGIS is to use XYZ Tiles by connecting to a tile service. Right click on XYZ Tiles >> New Connections. In the XZY Connection window, fill in the info and you should be good to go.

Fig0006.png
Fig0008.png

But if you prefer to use a script to add in multiple basemap sources at once, here is a python script to do just that -- it is super handy. Much thanks and appreciation to the many out there for sharing their work. Go to https://raw.githubusercontent.com/klakar/QGIS_resources/master/collections/Geosupportsystem/python/qgis_basemaps.py  and copy and paste this in the python console in QGIS.

  • To open the Python Console, go to Plugins Menu >> Python Console.
Fig0009.png
  • Then in the Python Console window, paste in the script that you copied from the URL linked above. Run it and you should see a list of basemaps under the XYZ Tiles in the Browser panel.
Fig0012.png

That's it! Add basemaps to your heart's content.

3D DEM Visualization in QGIS 3.0

The other day I just happened to be looking through ESRI’s ArcUser magazine (Winter 2018) while waiting for a process to finish on my computer, and came across an article on visualizing DEM using multidirectional hillshade -  Create Amazing Hillshade Effects Quickly and Easily in ArcGIS Pro.  Inspired by this article, I thought I would try do the same thing in QGIS 3.0. Below is a short tutorial on how to visualize a DEM in 3D with QGIS 3.0. I'm thinking that the more I use QGIS (especially the newest update) the more I like it.

A quick note on traditional vs multidirectional hillshade:

Traditional Hillshade: hillshading illuminated from 315 deg azimuth - so only single source of light

Multidirectional Hillshade: a combination of hillshading illuminated from 225, 270, 315, and 360 deg azimuth - 4 different sources of light

1. Download DEM

First, you'll need a DEM for your area of interest. For this tutorial, I downloaded a 1 meter DEM for Diamond Head Crater from NOAA Data Access Viewer. Below is the specs I used for my download.

DEMSpec.jpg

After you download your DEM, unzipped it. You're going to then use QGIS 3.0 to view it

2. Create Multidirectional hillshade

I am using QGIS 3.0 and I assume you have it installed already. So here's what you're going to do. First, add the DEM to the QGIS, then make a copy of it. You'll end up having 2 DEM layers for viewing purposes: 1) colored DEM used for drapping over 2) the hillshade. You are not actually creating a new layer or doing anything to the actual DEM.

  • In QGIS open add your DEM to the Layers Panel
  • Make a copy of the DEM. Right Click on your DEM >> Duplicate
  • Optional - In your Layers Panel, rename your second copy of the DEM to indicate that this layer will have the Hillshade effect applied to it - like the screenshot below:
Fig0001A.png
  • Turn off the top layer cuz you're going to create the hillshade effect on the second layer first.
  • Right click on Hillshade_oahu_dem_2013 >> Properties >> Style
    • Render type: Hillshade
    • Multidirectional: check the box
    • Resampling: Use either Bilinear or Cubic and Average to get a smoother looking hillshade. The default Natural Neighbor produced weird cross hatching effect.
Fig0004.png

Try to see what happens when you checked (multidirectional hillshade) and unchecking (traditional hillshade) the Multidirectional hillshade box. Here's are two views of a traditional hillshade and a multidirectional one. It looks to me like there is less shadowing  in the multidirectional hillshade (which would make sense if illumination is coming from multiple directions) then the traditional one - this though I think makes it have less contrast (?) so things don't seem to pop out as much??

Fig0025.png
Fig0026.png

3. Drape Colored DEM Over Hillshade

After you create the multidirectional hillshade, you're going to just change the color of the top DEM layer and "draped" it on top.

  • Right click on the top DEM (e.g. oahu_dem_2013) >> Properties >> Style
    • Rendered type: single pseudocolor
    • Make sure min/max values are correct for your DEM
    • Color ramp: Use drop down menu >> Select all Color Ramps >> BrBG (this is the one I'm using). To invert color ramp, right click inside the color box >> Invert Color
    • Blending mode: Multiple (this produces better effects than Normal blending with transparency)
    • Resampling: again use Bilinear or Cubic to get smoother effect
    • Optional Transparency: change transparent if you want. I have it at 50%.
Fig0018.png

Here is a view with the colored DEM draped over the multidirectional hillshade. Play around with the setting options to get your desired effect.

Fig0020.png

For comparison, here's the same view I did in ArcGIS Pro. Both results look really good - although I think ArcGIS Pro's multidirectional effect (using default settings) has a slight edge over the QGIS one (also using default settings).

Fig0019.png

4. Create 3D Map View

You can now view/simulate 3D landscape (DEM) in QGIS 3.0 natively without a plugin.

  • In QGIS, go to View menu >> New 3D Map View
  • In the 3D Map 1 window: click the Configure button
    • Elevation: select your DEM or the hillshade
    • Make any changes to the other settings if you want
  • In the 3D Map 1 window: hold down the shift key and the left mouse button to zoom in/out and rotate
Fig0028.png
Fig0029.png

Did you also know that you can customize the QGIS user interface by stacking and moving panels and map views around - try dragging and docking the panels in your user interface. The image below shows my user interface setup with the main map view or canvas on the left and the 3D Map view on the right.

Fig0022.png

Here is one of my favorite 3D visualization where I burned an aerial (ESRI Aerial base map) into the multidirectional hillshade in QGIS.

Fig0023Burn.png

QGIS and ArcGIS Pro DEM Comparisons

My conclusion having used the same method for 3D DEM visualization in both QGIS 3.0 and ArcGIS Pro is that they're both very good and comparable. Of course you can do some comparisons and decide for yourself, but I think QGIS does a great job considering it is free. Anyway, below are some of my snapshots for comparison.

Fig0012.png
Fig0016.png

That's it for this post. I hope you enjoyed reading it. Feel free to drop me a note in the comments section or via email anytime.

A Simple Leaflet Web Map Example

By now I think everyone has heard of Leaflet. I've been playing around with it on and off over the last year. For this post, I thought I'd write about a simple Leaflet map example, using Leaflet and ESRI Leaflet plugin. I am not a coding expert by any means, just someone who likes to learn, so I try to annotate/comment in my code as much as I can. This post is not a step-by-step guide in creating a Leaflet Maps but rather comments and notes on what I've learned in the process of using/learning Leaflet. However, feel free to download the source files and use it however you like.

The Leaflet example uses several plugins and calls layers from the Hawaii Statewide GIS Program's ArcGIS REST service and the Pacific Islands Ocean Observing System (PacIOOS) Geoserver. You can view and download the source files on my GitHub account and also view the demonstration site as well.

Demo Web Map: https://stephsaephan.github.io/leaflet-map-example/

For the demo web map all of the data layers are turned off by default with exception of a basemap. I also include OSM 2.5D builidngs to the map as well.

Source files: https://github.com/StephSaephan/leaflet-map-example

Fig0011.png

3 Lessons I Learned:

plugins

Not all plugins will work well together, even if they come from the same source. Case in point, ESRI Render Plugin will render layers using the default symbology as defined by ArcGIS REST service. This plugin though doesn't currently work with the ESRI Cluster Plugin

Legends in Leaflet

Legends are harder to implement than it seems. I didn't include legends on the demo web map, but I have tried adding legends from ArcGIS and GeoServer. See images below.

I've used this ESRI Leaflet Legend plugin to call a legend for a feature service. Below is a a example code snippet for calling ArcGIS REST Legend. I didn't go any further than just trying to see if it'll work.

Fig0007.jpg
 ArcGIS REST Legend displayed on map, but not exactly working the way I want it to.

ArcGIS REST Legend displayed on map, but not exactly working the way I want it to.

Here is another code snippet for displaying a legend from a Geoserver WMS using the Get Legend Graphics request. This is calling from an in-house geoserver I set up.

Fig0009.jpg
Fig0010.jpg

Secured (https) and unsecured (Http) items

I have the demo leaflet web map hosted on GitHub - which is free to use if the repository is public, and by default will be secured. When I access the demo site, I get the warnings:

Internet Explorer:  This site contains both secure and nonsecure items.  Do you want to display the nonsecure item?

Firefox:  Parts of this page are not secure (such as images)

Chrome:  This page includes other resources which are not secure. 

This warnings occurs when items on the following web site (frames, images and other objects) are referenced in a non-secure manner (http).  For example the page code might include an image or code that pulls from a non-secure URL. Both the PacIOOS GeoServer and Hawaii Statewide GIS ArcGIS REST services unsecured - on http , so hence the warning messages. However, I noticed though that features from the PacIOOS unsecured geoserver will still display on the web map but not features on the ArcGIS REST service - snapshots are shown below.

 PacIOOS Geoserver is unsecured (hhtp)

PacIOOS Geoserver is unsecured (hhtp)

 Hawaii Statewide GIS ArcGIS REST service is unsecured (http)

Hawaii Statewide GIS ArcGIS REST service is unsecured (http)

 Leaflet web map hosted on secured site (https) but referencing items from unsecured (http) sites (PacIOOS and Hawaii State GIS).

Leaflet web map hosted on secured site (https) but referencing items from unsecured (http) sites (PacIOOS and Hawaii State GIS).

Solution is to use one of the following:

  • Link to https explicitly:

Some websites will use both http and https. If you can access the items by using https then use that instead. Testing explicit https,  I was able to access the Hawaii Statewide ArcGIS REST services (https://geodata.hawaii.gov/arcgis/rest/services) but not the PacIOOS Geoserver (https://geo.pacioos.hawaii.edu/geoserver)

  • Use relative linking if on your own domain
  • Use protocol relative linking to use images/items from other domains. This is the option I used in my Leaflet example, like so:

This works: '//geodata.hawaii.gov/arcgis/rest/services/Terrestrial/MapServer/34'

This does not work ( at least not from GitHub https, but it works locally when I do live preview locally from Brackets): //geo.pacioos.hawaii.edu/geoserver/wms?

So I kept the URL for the PacIOOS geoserver WMS as http as it displayed on the web map.

That's it for this post. Thanks for reading.