Masthead

Lab 10: Raster Analysis

Introduction

MtShasta aerial

In this lab, we will conduct a spatial analysis using remotely sensed data. We will use a recent Landsat 8/9 Level 2 Surface Reflectance Data of the Mt Shasta region (acquired in late February 2024), a Digital Elevation Model (DEM) and a watershed boundary of the region to conduct a raster based spatial analysis. The purpose of this analysis is to locate, quantify and analyze areas of snow cover within the seven watersheds that encompass Mt Shasta. The satellite data used is Landsat Collection 2 Level-2 Surface Reflectance Data. This Landsat data has been radiometrically calibrated and atmospherically corrected and the pixel values represent surface reflectance percentage.

 

 

Learning Outcomes

Set up your Workspace

First, we need to set up our workspace and transfer the data for the lab.

  1. On the Desktop create a new folder named “Lab_10”.
    1. Within the folder named “Lab_10”, create three subfolders as below:
      1. Originals
      2. Working
      3. Final
  2. Locate the GSP 216 Google Drive lab data folder and locate the file "Lab10 data" folder, download all of the files (three) and copy or directly extract all of the data into your "Originals” folder. Uncompress the files in your original folder (make sure they are extracted into separate folders). Now the files have been uncompressed and are ready to use.

Working with Landsat Surface Reflectance Data : Save and Subset Images and Calculate Ratio Index

  1. Start “ENVI ” and Select File → Open and navigate to your Lab 10 Originals folder where you extracted your data and find the folder that contains the Landsat data and open the metadata text file (..._mtl.txt) associated with the data.
  2. Open the file Shapefile "Shasta Watersheds". This data was extracted from the National Watershed Boundary Dataset (WBD) and includes the boundaries of the seven watersheds that encompass Mt. Shasta.
  3. Before we can process the data we will save and subset the surface reflectance data. This is necessary when working with Level-2 data. From the main menu select File → Save As.. (ENVI...). Select the file "Surface Reflectance" data from the list and click OK. Click the small arrow on the right hand side or click Spatial Subset. Select the Subset by Vector icon and select the watershed shapefile. This will save and subset the data. Save the file as MtShasta.dat in your Finals folder and click OK. Shasta Subset

Snow has high reflectance in the visible and low reflectance in the shortwave infrared. There are several different tools that can be used to enhance and extract spectral differences. We will test out a spectral index, band ratio and band difference to compare the reflectance of two bands.

  1. Open the Data Manager and remove all the original Landsat files and leave the MtShasta.dat and Shasta_Watersheds shapefile. Change the band combinations of the MtShasta.dat image to Red - SWIR, Green - NIR, Blue - Red. This will make the snow appear bright blue and easier to differentiate from other land cover.
  2. Use the spectral profile tool to investigate the reflectance curve of some of the snow areas. Make note of where in the spectrum the highest and lowest reflectance, these are the bands that should be used when selecting appropriate spectral enhancements. We will try several different methods and select the best results for differentiating snow cover.

Normalized Difference Snow Index (NDSI)

  1. From the Toolbox, select Band Algebra → Spectral Indices or type "Spectral Indices" in the Toolbox and open the Spectral Index tool. The Spectral Indices dialog appears. Select the MtShasta.dat is your Input Data and click OK. In the next window select the Normalized Difference Snow Index from the list. Instead of saving the output raster as a file select Virtual Raster. This will save the file in temporary memory allowing us to review the results without permanently saving the file.
  2. Check out the results of the NDSI calculation. Higher values generally indicate snow, although a positive NDSI value is indicative of any material that has lower reflectance in the SWIR than in the Green wavelengths. Compare the results of the NDSI to the false color composite image. Note that while theoretically NDSI values should range from -1 to 1 values outside this range can occur if there are negative values in the original data. This can occur in some atmospherically corrected imagery and sometimes additional processing is required.

Calculate the Band Ratio

  1. In the toolbox select Band Algebra → Band Ratio. For the MtShasta.dat file, select SRB2 (0.4819) which corresponds to the blue wavelength as the Numerator and SRB6 (1.6081) the SWIR1 as the Denominator and click Enter Pair.
  2. In the next window save the output results to Memory rather than file. This will create another virtual file that we can use to view the results. Keep the rest of the defaults and click OK.
  3. Compare the results of the band ratio to the false color composite and snow cover. Larger values indicate higher ratio between the visible and SWIR.

Calculate the Band Difference

  1. From the Toolbox, select Band Algebra → Band Math or type "Band Math" in the Toolbox and open the Band Math tool. The Band Math window will open. Enter the following expression: float(b1 - b2) and click OK. This will calculate the difference between any two data sets (the same expression was used to calculate dNBR).Band Math Expression
  2. The Variable to Bands Pairings window will open, this is where you will we tell ENVI which bands or data to use for the variables (b1 & b2). For B1 select SRB2 (0.4819) which corresponds to the blue and select SRB6 (1.6081) as B2, saving the output results to Memory again rather than to a file.
  3. Compare the results of the band difference to the false color composite and snow cover. A large positive value in the band difference is indicative of snow. This layer should give the a high degree of contrast between snow and other land cover. Compare all three of the images (NDSI, Band Ratio and Band Difference) and select the image that best differentiates snow from all other land cover types.
  4. Open the data manager and remove of the files except the spectral enhancement file you selected in the previous step (NDSI or Band Ratio or Band Difference), the subsetted Mt Shasta Landsat image and the watershed shapefile.

Sample Areas of Known Snow Cover

To better isolate the snow covered areas in the image we will use reference data. These are locations with known snow cover. This will help us decide what values represent snow and to select an appropriate threshold.

  1. In the main toolbar, open the file "Shasta_Snowcourses.shp" from the originals folder. This shapefile contains points with the locations of five Snow Course stations that were sampled for snow depth and water content by the California Department of Water Resources. All of these sample points had snow at the date of the imagery. We will use these points to sample the band difference values for areas that we know have snow.
  2. In the Layer Manager right click on the spectral enhancement file and select New Region of Interest. In the ROI Tool window select File → Import Vector. In the Data Selection window select the "Shasta_snowcourses.shp" file and OK.Import Vector
  3. In the Convert Vector to ROI window select the "All records to single ROI" and click OK. This will create one ROI with five points. Select the Band Math layer as the Visualization Layer if given an option.
  4. Now the ROI has five points that represent areas that have snow on the ground. In the ROI Tool window click the Compute Statistics button. The will compute the mean, minimum and maximum ratio values of the five snow sampling points. Write down the minimum value, this is the value we will use as the snow threshold (i.e. areas with a value greater than the minimum value should be snow and areas with values lower than the minimum value are not snow covered). Close the ROI window when you are done.

We will use the Raster Color slice to create a single class that isolates the snow in the image by selecting a range of values.

  1. In the Toolbox navigate to Classification → Raster Color Slices. Select the spectral enhancement file as the input and click OK.
  2. Delete all of the existing slices by clicking Clear Color Slicing. Then create one new Color Slice by clicking Add Color Slice. We will only be creating one class which represents snow.Color Slice
  3. Use the value from step 20 as the slice Min and keep the Slice Max as the default value (should be the maximum value found in the image). Snow Class Zoom in on your image and turn the Raster Color Slice layer on and off to preview your selection. Make adjustments by either changing the Slice Min value or by manually dragging the bar on the histogram to make adjustments to the range of the class/slice based on visual analysis. Try to make sure your selection includes all of the snow and no other features.
  4. In the Raster Color Slice window export the data by clicking the save icon and selecting Export as Class Image. Save the file as snow.dat in your Working folder.
  5. The last step is to mask the snow layer you just created to the boundaries of the watersheds. From the main menu select File → Save As (ENVI..) and select the snow.dat file as the input file. Click the Mask option and select the Shasta Watershed shapefile. Save the file as "shasta_snow.dat" in your Final folder.

Download and Process Digital Elevation Models

  1. Go to The National Map download interface: https://apps.nationalmap.gov/downloader/. A variety of data can be downloaded from the National Map interface. We will be looking for Digital Elevation Model data. You can also download Hydrography and Boundary as well as a variety of other data from the National Map.
  2. Data can be selected by Map Extent, Polygons or for specific coordinates. Shapefiles can also be uploaded for visualization purposes and then you may also use the polygon or extent tool to draw a rough outline around the shapefile of interest. Click on the Upload Shapefile button and select the zipped Shasta Watersheds shapefile (Shasta_Watersheds.zip). You should now see the Mt. Shasta watershed boundaries on the map.
  3. Click the extent button to draw a bounding box around the shapefile on the map. This will limit the data search to the Mt. Shasta watershed boundaries.
  4. Under Data on the right hand side check Elevation Products and select ⅓ arcsecond DEM. ⅓ arcsecond corresponds to approximately 10 meter spatial resolution. Keep the Data Extent as 1 x 1 degree and the File Format as GeoTIFF, IMG. Click the blue Search Products icon to search for data.
    5. Elevation Data on TNP
  5. In the Available Products window there should be at least two datasets, one directly cover Mt Shasta and one to the east, the corresponding map index names are: Weed E and Alturas W. Click the Footprint to make sure they cover the area. Download the most current GeoTIFF dataset for each location: USGS 1/3 arc-second n42w123 1 x 1 degree and USGS 1/3 arc-second n42w122 1 x 1 degree These files cover the area of the Mt. Shasta Watersheds.
  6. Find the datasets in your downloads folder and move them to your Originals folder. The DEM files will be .tif (GeoTIFF) files.
  7. We want to mosaic the two DEMs together into one seamless file for analysis. In the ENVI Toolbox select Mosaicking → Seamless Mosaic. This will open the Seamless Mosaic window. Note that your viewer may be blank after opening up Seamless Mosaic tool, this is normal and all of your files will re-appear when the process is complete. Click the "+" icon to add the scenes be mosaicked. Click the folder icon to locate the DEM files. Select your two DEM tif files from your originals folder and click OK (use the shift button to select more than one file). You will see the two DEMs in the viewer. Select Files
  8. Click on the "Export" tab. Under output filename name your file "DEM.dat" and save it in your "Working" folder. Check the "Show Preview Box" to see a preview of the mosaicked image. Once the mosaic process is complete the new DEM will appear in your viewer.
  9. When the new mosaicked DEM appears it will say [RE-PROJECTED], which indicates this file has a different spatial reference system than the other files. Therefore we will need to re-project the DEM raster. The Regrid Raster tool in ENVI is one way to accomplish this task. In the Toolbox, select Raster Management > Regrid Raster.

    The Regrid Raster tool allows a user to resampled and reproject raster files to a specified spatial grid, which includes a specific coordinate system and pixel sizes. Making sure all data have the same coordinate system and pixel size is necessary for many raster analysis techniques.

  10. In the Regrid Raster tool window select the DEM as the input raster. Under Grid Definition leave the default option (Coord system + extents + pixel size) and click the From Dataset icon from dataset. Select the shasta.dat multispectral file as the reference dataset. This will set the coordinate system and cell size based on this dataset. Click OK. The window should now update showing the coordinate system, extent and pixel size from the Landsat data.
  11. Save the file as DEM_UTM.dat in your finals folder and click OK to run the tool. When the process is finished a the re-projected and re-sampled DEM should appear the Layer Manager, this time without saying [REPROJECTED]. Check to make sure the following files are in your finals folder: Shasta.dat (subsetted Landsat image), Shasta_snow.dat and DEM_UTM.dat. You may close ENVI when you are done.

Raster Spatial Analysis in ArcGIS

We will use the Snow file, Watershed Boundary shapefile and the Digital Elevation Models (DEMs) to analyze the extent and area of the snow. The first step will be to determine the amount of snow in acres within the seven Mt. Shasta Watersheds. Then we will conduct a more in-depth analysis including looking at the snow cover by aspect (slope direction) and elevation. Before using many of these tools you will need to make sure the Spatial Analyst extension in enabled in ArcGIS.

Determining Area of Snow Layer

  1. Start ArcGIS Pro and create a new map project in your lab finals folder. We will need to make sure the Extensions are enabled so we can use all of the tools. In the Project tab click Licensing. You should see a list of Esri Extensions. Scroll down and click Configure your licensing options. Turn on the Spatial Analyst license by checking the licenses box.
    2. Licensing
  2. In the Geoprocessing Pane search for or select Toolboxes →Data Management → Raster → Raster Properties →Build Raster Attribute Table. Select the Shasta Snow layer (shasta_snow.dat) snow layer from your finals folder and click Run. This will add an attribute table to the files that include the class value (1, indicating snow) and the number of pixels in each class, this file should be added to the map after the process is complete. Note that inorder for the Build Raster Attribute Table function to work properly it must be done before adding the data to the map.
  3. Add the Shasta watersheds shapefile, the satellite imagery (MtShasta.dat) and the UTM DEM file from your map.
  4. Make sure the Spatial Analyst extension is enable (see below box). In the Geoprocessing pane search for Zonal Histogram or select Spatial Analyst Tools → Zonal → Zonal Histogram. The Zonal Histogram will calculate the number of pixels (count or frequency) in one dataset within classes or regions of another dataset.


    5. Tool No Licensed Error If you receive an error that looks like the one to the left, you need to enable the "Spatial Analyst" extension before using any of the Spatial Analyst tools. See step 37 for instructions on how to enable extensions.

  5. In the Input Raster or Feature zone select select the Shasta watershed shapefile. In the Zone Field select Name. This will summarize the number of snow pixel by the watershed name. In the input value raster select the Shasta Snow layer and save the table as wateshed_snow in your final folder. This table will summarize the number of snow pixels within each watershed.Zonal Histogram
  6. The table will automatically be added to your Table of Contents. Right click on the file and select open. The table shows the number of snow pixels in each of the seven watershed boundaries. We will use this information to calculate the area of snow in acres. Select the row with the data and copy the data and paste into Excel.
  7. Check the spatial resolution on the Snow layer by right clicking, selecting Properties and writing down the Pixel/Cell Size. Use the below guide to calculate the area of snow in each watershed in acres and the total area within all seven watersheds.

    8. Example: Calculating Area of Raster Class
    First determine the spatial resolution of the raster and the number of pixel in the class or category
    Total Area of class in square meters= spatial resolution2 x pixel count for class
    Note that 1 acre = 4046.86 m2

Create an Aspect Raster

The aspect of terrain refers to the direction it’s facing and is computed in degrees from due north. Aspect relates to the amount of sun exposure, in the northern hemisphere a north-facing slope faces away from the sun and therefore is generally cooler and wetter than a south-facing slope. The aspect of a slope can have significant influence on its micro-climate and vegetation. Given this influence on vegetation, aspect is frequently used in raster analysis. Aspect is always given as a direction in degrees, for example a North facing slope the aspect would be between 0-22.5 degrees.

  1. From Geoprocessing Toolbox, select Spatial Analysis Tools → Surface → Aspect.
  2. Set the "Input Raster" to be your Elevation layer (DEM). Save the "Output" raster to your working folder and name it "Aspect.tif”. Leave the default method as PLANAR and click Run to create the file.Aspect Tool
  3. You will now see your Aspect raster in your window. The Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors . Aspect can be thought of as the slope direction. The values of the output raster are the compass direction of the aspect, a color coded classified symbology is used to display the degrees as directions.

    4. Aspect in ArcGIS
    AspectBy default aspect files are usually generated with a color coded scheme where the degree values are grouped into directions (North, Southeast etc). You can also add the Colormap function to display the dataset using the color map renderer or a specified color scheme. Use the following parameter settings to produce consistent color mapping for your aspect datasets: Raster—Your aspect output Color Scheme Type—Color Ramp Color Ramp—The first color ramp (named Aspect) Click Create New Layer.

Create a Hillshade Raster

Hillshades are common background used in maps. A hillshade is a 3D representation of a surface using the sun's relative position to provide shading and depth to the image. A hillshade is a raster that essentially appears to be a photo from above the sea, shaded by the sun. Hillshades are not used for analysis but provide nice backgrounds for maps and data.

  1. Navigate to the Geoprocessing → Spatial Analysis Tools → Surface → Hillshade.
  2. In the Hillshade dialog box select the elevation layer (DEM)as the input raster file. Save the output raster in your Finals folder as "hillshade.tif”. Keep all of the default settings. Click OK to create the hillshade raster. You should now see the hillshade file in your viewer.
  3. Move the hillshade to the bottom of your layers for now. We will need it to create our map later.

Analyze Snow Cover by Aspect

  1. Run the Zonal Histogram tool to determine the number of snow pixels in each aspect class (North, Northeast etc). For the Input Zone, select the snow layer (shasta_snow.dat) and the aspect layer as the Input Value Raster. When selecting the aspect layer use the dropdown arrow to select the file. This will select the file with the classified symbology (North, East etc). Save the table as "aspect_snow" in your finals folder.aspect zonal histogram
  2. Use the same process as step 40 to calculate the area in acres of snow for each of the nine aspect classes.

Determine Snow Cover Elevation Range

Zonal Statistics by Table: This tool summarizes the values of a raster within the zones of another dataset (either a raster or shapefile) and reports the results to a table. The stats include the mean,minimum, maximum values etc.

  1. Go to the Geoprocessing Toolbox , select Spatial Analyst Tools → Zonal → Zonal Statistics as Table. Select the Snow layer (shasta_snow.dat) as the Input Zone (leave the file as Value) and Select the Input Value Raster as the elevation (DEM). This will provided statistics about the elevation data within the area where snow is present. Save the output table as elevation_stats in your final folder.
  2. Open the table and copy the statistics into your Excel spreadsheet. Note the minimum and maximum elevation values. This tells us the minimum and maximum elevation that snow was found at in the image.

Analyze Snow Cover by Elevation Range

We want to create an output raster that shows areas that have snow by elevation range. We will use the Raster Calculator to accomplish this.

  1. We want to divide up the elevation into three elevation classes, low elevation snow, mid and high elevation. The elevation ranges will be:
    • Low Elevation: 0-1500 (less than 1500m)
    • Mid Elevation: 1500-2999 (more than 1500, less than 3000m)
    • High Elevation: 3000-4500 (greater than 3000m)
  2. Start the Raster Calculator by going to the Geoprocessing , select Spatial Analyst Tools → Map Algebra → Raster Calculator.
  3. We want a raster that will display all of the areas that have snow and that are in the low elevation range (as determined above). Double-click on the snow file to have it added to the equation box, click on "==" and then type 1. This will select for area that have snow (values of 1). Be sure to also add parentheses ( ) around your expression. Now click the “&” operator. Double-click on the elevation layer and then click on "<" and then type 1500 (the desired elevation value), again be sure to also add parentheses ( ) around your expression. Save the output raster as “lowelev_snow.tif” in your Lab 10 “Final” folder. Your dialog window should look like the one below. Raster Calculator
  4. The raster file should have automatically been added to your Table of Contents. Remove 0 class from display by going into the Symbology and removing the "0" class. The areas show in color (1) meet our criteria, i.e. low elevation snow cover. You can also label the "1" class as Low Elevation Snow. Remove Value
  5. Repeat Raster Calculator process to create a raster layers for the remaining two elevation classes, mid elevation snow and high elevation snow. Keep the first first of the expression the same ( "shasta_snow.dat" ==1) but adjust the elevation criteria using the following guidelines:
    • Mid Elevation Range: Greater than or equal to 1500 and less than 3000. Note that each you will need to put () parentheses around each expression component and each elevation criteria needs to be a seperate criteria (e.g. ("DEM" >500) & ("DEM" >=1000)
    • High Elevation Range: Greater than or equal to 3000
  6. For each of the three analyses, right click on the layer and select Open Attribute Table. Record the number of pixels that are classified as snow (1) in each of the analyses. You will also need to verify the spatial resolution of the data. This can be done by right clicking on the file and selecting properties. In the Layer Properties window click the Source tab and note the spatial resolution (cell size).
  7. Use the above information to calculate the area of snow in acres for each of the elevation ranges. You will use this information to create a bar graph in Excel or Google Sheets showing the area of snow in each elevation class.

Create a Map Displaying Snow Cover by Elevation Range

  1. Re-arrange the layers so that the three snow elevation layers are on-top, with the hillshade on the bottom. You can remove or hide all of the other layers. Remove the default esri basemaps and hillshade as well.
  2. Select the low elevation snow layer in the Content and click on the Raster Layer tab in the main toolbar. Under “Transparency” change it to “50%”.  Change the transparency for all three of the snow layers to 50% as well. You may adjust the transparency % as you see fit.
    Shasta Snow by ElevationThis is a example figure for 2017 data, your final image will look different.
  3. Create a map layout that include a title, north arrow, scale and legend. The legend should include the three snow elevation classes labeled with the class name. You may focus your map on Mt Shasta if you like or include the entire study area.
  4. Save your map and export the map as a PDF or JPG. Back up your final folder when you are done. Make sure had all of the data needed to complete the questions. Use Excel/Google Sheets to perform the calculations and to create tables and graphs for the report. A tables and graphs should be properly formatted and easy to read.

Lab (Upload to Canvas)

You will be submitting two things for this lab, one is a map and the second is a written document answering the below questions.

Snow Map:

Questions: Complete a document that includes the following (it does not have to be in the format of a lab report):



Contact Info

Humboldt State University
1 Harpst Street Arcata, CA 95521
skh28@humboldt.edu

© Copyright 2020 HSU - All rights reserved.