In Lab Exercise 1, we learned that the metadata file distributed with each Landsat 8 scene contains information needed to convert the raw and uncalibrated digital numbers (DNs) contained in Landsat band images into either radiance or reflectance values. Certain remote sensing applications (e.g. image classification) require the use of calibrated imagery data containing reflectance values, and sometimes, radiance values. Therefore, it is important to know how to transform the raw imagery into calibrated form.

3.1. What are the differences between image DNs, radiance values, and reflectance values? (3 marks)

The process of calibrating raw satellite imagery into reflectance values may include several steps, often including a conversion of DN to radiance, a model for removing the effects of varying illumination within the scene (due to the sun angle, the earth-sun distance, terrain, and clouds and cloud shadows), and an atmospheric correction model. There are varying levels of sophistication that can be used during calibration and exact process will also vary depending on the exact satellite and sensor (Landsat 8 and the OLI sensor in our case). We are going to examine the basic method used to convert Landsat 8 Level 1 DNs into reflectance values. Based on this useful documentation (please read this over), we can see that the basic conversion to top-of-atmosphere reflectance values, including a correction for sun angles, is:

ρ = (M * DN + A) / sin(θ_SE) (Eq. 1)

where ρ is the transformed reflectance value (at a wavelength), M is the band-specific multiplicative rescaling factor, A is the band-specific additive rescaling factor, and θ_SE is the sun elevation, in radians. M, A, and θ_SE can each be found within the metadata file distributed with a Landsat scene and are specific to the geographic location of the scene and the time that it was acquired.

3.2. What is difference between the top-of-atmosphere reflectance and surface reflectance? (2 marks)

Open the metadata file associated with the lab data, LC08_L1TP_018030_20170624_20170713_01_T1_MTL.txt, using a text editor such as VS Code. Try to locate the sun elevation and each the reflectance rescaling parameters in the file.

3.3. How do the reflectance rescaling parameters compare to the radiance rescaling parameters? What do you think this difference may imply about the pre-processing that was done on the imagery before it was made available? (3 marks)

We could use Eq. 1 above to convert each of the bands within our scene data set manually. However, this would be a tediuous task and the opportunities for using the incorrect parameter by accident is fairly high. Instead, we're going to write a script to automatically read the parameters from the metadata file and then apply the DN-to-relectance transform to each of our bands. Notice, it is very important that you download the latest version of WhiteboxTools before attempting this part of the lab. There have been updates to the library since the last lab assignment that affect the following analysis. Using VS Code, create a Python script in your Lab 2 folder called reflectance.py. Copy the following code into your script file and run it using the VS Code terminal and don't forget to modify the wbt.work_dir appropriate for your data location.

from math import radians, sin
import os
from WBT.whitebox_tools import WhiteboxTools


def main():
    wbt = WhiteboxTools()  # Initialize WhiteboxTools
    wbt.work_dir = os.path.join(os.path.dirname(
        __file__), "lab2_data")  # Update this
    wbt.verbose = False  # This way the script won't output so many updates about progress

    # Open the metadata text file associated with the Landsat scence
    metadata = os.path.join(
        wbt.work_dir, "LC08_L1TP_018030_20170624_20170713_01_T1_MTL.txt")
    fh = open(metadata)

    # These two lists will contain the rescaling parameters.
    mult_term = []
    add_term = []
    for line in fh:
        # Read the file line-by-line looking for the reflectance transformation parameters
        if "REFLECTANCE_MULT_BAND_" in line:
            mult_term.append(float(line.split("=")[1].strip()))
        elif "REFLECTANCE_ADD_BAND_" in line:
            add_term.append(float(line.split("=")[1].strip()))
        elif "SUN_ELEVATION" in line:
            # We're also getting the sun elevation from the metadata. It has
            # to be converted to a float and radians.
            sun_elevation = radians(float(line.split("=")[1].strip()))
    fh.close()  # Be sure to close an open file

    # We'll transform the first five bands of the data set
    for i in range(0, 5):
        # The bands are indexed 1, 2, 3... but the mult_term list is indexed 0, 1, 2...
        band_num = i + 1
        print("Transforming band {}...".format(band_num))
        inFile = "band{}.tif".format(
            band_num)
        # This is where the DN-to-reflectance transform equation happens.
        # It creates two temporary rasters that are continually over-written.
        wbt.multiply(inFile, mult_term[i], "tmp1.tif")
        wbt.add("tmp1.tif", add_term[i], "tmp2.tif")
        # The final transformed raster is called bandX_reflect.tif
        wbt.divide("tmp2.tif", sin(sun_elevation),
                   "band{}_reflect.tif".format(band_num))

    # Provide some sort of indication that the job is done.
    print("Operation complete!")


main()

This script may take a fair bit of time before completing but it will at least provide you with periodic updates on progress (e.g. Transforming band 1...). Hopefully, the script runs without issue for you. It is however one of the longest scripts that we've seen so far and because it involves reading and parsing text data, there is the potential for system-specific bugs to occur. If you encounter errors when you run the above script, examine the script carefully to see if you can identify the source of the bug. Once you have gained a full understanding of what the script is doing (see the detailed comments), and have exhausted debugging options (e.g. adding print statements in the script to see how far it progresses, searching the Internet for solutions), then you may approach fellow students or the GTA for assistance.

Once you have successfully run the script, display each of the five transformed bands using your visualization software. You will notice that they look very similar to the original bands, however, their value range is very different. Each pixel now contains an approximate reflectance value.

3.4. What are the minimum and maximum values within each of the transformed images. Are these value ranges expected for reflectance data? (6 marks)

3.5. Describe the specific ways that you would need to modify the above script if you wanted to perform a DN-to-radiance transformation instead of the DN-to-reflectance we performed. (Hint: be sure to read over the documentation before answering.) (4 marks)