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. At the time of publication of this tutorial, the current version is QGIS 2.18.15 'Las Palmas' , originally released on August 12th 2017. The most up to date version can be downloaded using the following link.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.

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, Landsat 8 data 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

We can now begin working with the data in QGIS. Begin by loading the different landsat bands into QGIS by using the Add Raster Layer button on the menu on the left of the user interface

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.

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 tool, found under Raster, Conversion, Translate. Be sure to 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.

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, except in this case 1 equates to areas of ice and 0 has no ice. We can reclassify this output using the Translate tool again. 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.