Remote Sensing Analysis in QGIS

Introduction

Remote sensing is a broad discipline involving the observation of an object or phenomenon without physically interacting with it. In a remote sensing context, this generally involves the analysis of remotely sensed images from a variety of sources such as RADAR or multispectral satellite imagery. However, many commercially available remote sensing softwares can be very expensive and not accessible to the average user outside of major private or public institutions.

The purpose of this tutorial is to showcase the remote sensing capabilities of the free GIS software QGIS. To do this, we will be going through a workflow to delineate icings from Landsat imagery. Icings, thick masses of ice that build up through successive overflows, can endanger remote communities or vital infrastructure. The workflow was originally designed using expensive proprietary software such as ENVI and ArcGIS. By reproducing the finished products using free software, this tutorial hopes to demonstrate the capabilities of QGIS in a real scientific application.

Obtaining Data and Software

To begin, the user must download the most recent version of QGIS. The current version is QGIS 3.2.3. 'Bonn', originally released on September 14th 2018. The most up to date version can be downloaded using the following [1]. Select your operating system and ensure that you download the proper version (32-Bit or 64-Bit for windows machines).

Landsat imagery is preprocessed and made publically available by the United States Geological Survey (USGS). Images can be browsed and downloaded through the EarthExplorer application. To get access to the imagery you need to create an account.

In EarthExplorer, we need to narrow down the vast repository of data that is available. Firstly, we must filter the data to only include summer images. This increases our ability to extract icings from the landscape, as there is little to no snow left on the ground.

Step 1: Navigate to the Search Months dropdown menu at the bottom of the Search Criteria tab and unselect all months except May, June, July and August as seen in the screenshot below.

Step 2: For our purposes, we should use the imagery with the broadest spectral range. Launched in 2013, data from the Landsat 8 platform should be suitable. Select the Datasets tab at the top of the page and navigate to Landsat, Landsat Collection 1 Level-1 and select Landsat 8 OLI/TIRS C1 Level 1.

Step 3: We can now narrow down the search in even greater detail using the Additional Criteria tab. We need to find an area with visible icings, commonly near the bottom of mountain slopes at high latitudes, where the hydrological and climate conditions favor icing formation. Landsat missions split the earth up into a grid of paths and rows for easy image identification. After some searching, a scene with good icing potential can be found at path 112 row 14, corresponding to an area in Siberia. We can also choose to limit cloud cover to 20% in order to minimize noise in the imagery. Enter these criteria as shown below.

Step 4: Finally, the results of our search can be seen under the Results tab. We can select to show the footprint of this image, view a preview and download any of the resulting images. For this tutorial, download the scene from May 15th 2013. Snow in the valley has mostly melted but large icings are still present, making it ideal for our purposes.

Loading Imagery and Atmospheric Correction

Make sure that you have unpacked the images twice in your Explorer. We can now begin working with the data in QGIS. Begin by loading the different landsat bands into QGIS. At the top of the user surface press Layer then Add Layer and choose Add Raster Layer.

Semi-automatic Classification Plugin

While the data obtained from EarthExplorer is very good, it still needs some correction applied in order to remove the haze and scattering that occurs when images are taken through earth's atmosphere. To do this, we will utilize the Semi-automatic Classification Plugin, or SCP within QGIS.

Navigate to the Plugins' menu at the top of the user and click Manage and Install Plugins. In the search bar, begin typing SEMI and the semi-automatic classification plugin will appear. Click the Install Plugin button on the bottom right begin using the plugin, as seen below.

Once the SCP is installed, activate it from the Manage and Install Plugins window. There should now be an SCP drop-down menu along the top of the user interface. This plugin contains a variety of different remote sensing tools beyond its classification abilities. Select SCP, Preprocessing, Landsat to open up the atmospheric correction window.

Once the window is open, enter the directory in which you've unzipped the Landsat files obtained from EarthExplorer into the Directory containing Landsat bands option. As this file should already contain the metadata file (.MTL) for the imagery, there is no need to specify its location. The metadata including date, elevation and the Earth to Sun distance should be automatically filled in.

After ensuring that all of the bands and metadata have been recognized by the plugin, check the Apply DOS1 atmospheric correction and click the Run button in the bottom right. You will be prompted to choose a directory for the output files, which should be named appropriately (eg. ATCOR).

Our imagery is now ready to work with!

Calculating Indices in Raster Calculator

The different bands of Landsat 8 correspond to different spectral wavelengths. By performing calculations on different combinations of bands, we can highlight specific features within a scene through the use of an index. The most commonly used index is most likely NDVI, or the Normalized Difference Vegetation Index, which is used to highlight areas of healthy vegetation. Using a similar concept, we will calculate various indices that highlight ice, snow and water in order to effectively extract icings from our scene.

NDSI

To begin, open the Raster Calculator by navigating to the Raster menu at the top of the page.

The first index we will calculate is the NDSI, or the Normalized Difference Snow Index. This is an index designed to tease out snow from the landscape. Our scene still has a lot of snow present on mountain peaks, so there should be a lot of snow present in the scene.

The formula for NDSI is ((Band 3 - Band 6) / (Band 3 + Band 6)), and can be entered into the raster calculator as seen below. Ensure that the expression is valid and that you have set the appropriate path for the output file. To identify your file later easily, I suggest naming it NDSI.tif.

The output of this calculation should look like the image below.

From this output, we must now extract the range of values that correspond to snow. Using the Identify button along the top of the user interface.

We can assess the values of the index and come up with an appropriate range of values that can be considered snow. In this case, values over 0.3 and less than 0.89 can be considered snow.

The output should give two classes of values ; a value of 1 where the conditions were met (aka snow is present) and 0 where there is no snow. We can now reclassify this output to isolate the class containing snow using the Translate (Convert Format) tool, found under Raster, Conversion, Translate. Be sure to select the correct Input Layer. Enter 0 in the field under Assign a specifierd nodata value to output bands. Define the Output data type as Byte. This will classify everything that did not meet the conditions of the raster calculator to NoData, or blank. Also be sure to specify an appropriate output folder.

The output file should look like this:

MDII

The second index we will calculate is the MDII, or the Maximized Difference ICE Index. This is an index designed to tease out pure ice from the landscape. The scene we chose is in early summer, where ice is still at its maximum extent but most surface snow in the valleys has melted. This should help to reduce noise in the final image.

Reopen the Raster Calculator from the top menu

The formula for MDII is ((Band 3)^2) - (Band 6)^2)), and can be entered into the raster calculator as seen below. Ensure that the expression is valid and that you have set the appropriate path for the output file.

The output of this file is in the same structure as the NDSI index in the previous section. From this output, we must now extract the range of values that correspond to pure ice. Using the Identify button along the top of the user interface, we can assess the values of the index and come up with an appropriate range of values that can be considered ice. In this case, values over 0.25 and less than 0.5 can be considered ice.

We can reclassify this output using the Translate tool again, except 1 equates to areas of ice and 0 has no ice instead of snow. click the No Data button and set the value to 0. This will classify everything that did not meet the conditions of the raster calculator to NoData, or blank. Also be sure to specify an appropriate output folder, possibly in an Indices folder alongside the NDSI output for organizational purposes.

NDWI

The final index that will be calculated is the NDWI, or the Normalized Difference Water Index. By extracting water from the image we can mask out bodies of ice that occur over water, thereby ensuring that all ice mapped occurs over land. This type of ice is of the most concern to this area of study, as ice spilling out over the land poses the largest potential threat to infrastructure in vulnerable areas. NDWI can be calculated using the expression ((Band 3 - Band 5) / (Band 3 + Band 5)).

Again, we can select a range of values from this index that best represent water. In this case, the values corresponding to water are greater than -0.3 but less than 0. Then we can reclassify this output using the same steps as before to create a raster file containing areas of water. Some error stems from the abundant mountain snow in our scene, which is unavoidable due to the spectral properties of melting snow.

Final Products

For display and reference purposes, we can render our Landsat scene to show standard RGB colours. This can be done using the Band Stack function within the semi-automatic classification plugin. You can find the tool at the top of the user surface. The band designations for RGB in Landsat 8 are 4,3,2 respectively.

Once rendered, we can turn on the three extracted layers in order to see the icings clearly. To edit the colors, you can find the layers panel at the top of your layers window.

After clicking the symbol, at the right side of your screen appears the Layer Styling window as you see in the following image. There you can select the layers and allocate them colours. Make sure you raise the brightness, otherwise you won't notice a different.

Your final picture could look like this (depends on the colours yo have seleceted).

Conclusions

This tutorial has showcased some of the remote sensing capabilities of QGIS in the context of scientific analysis. The products generated could be used to asses how icing distribution changes over time, calculations of water stored in icings or threat posed to local infrastructure. I hope this guide may empower those without the means to purchase expensive software to pursue scientific endeavors in the future.