JHU Open Source GIS with Isaacs

Ref: Dr. Rachel Isaacs (2019). Open Source GIS. JHU Course AS.420.703.81. Email: risaacs6@jhu.edu  

_______________________________________________________________________

Summary

  • An overview of the fundamentals of mapping, photogrammetry, geographic information science, and remote sensing. 

  • Programs: QGIS (https://www.qgis.org/en/site/), Google Earth Engine (GEE), & R. 

__________________________________________________________________________________

Landsat Satellite Missions

  • Landsat Group 1 (LS 1-3): Equipped with Multispectral Scanner (MSS), which recorded data in 4 spectral bands; two visible and two near-IR (NIR). 

  • Landsat Group 2 (LS 4-7): Equipped with Thematic Mapper (TM) or Enhanced Thematic Mapper (ETM+) with finer resolution size, expanded spectral coverage, additional bands in the middle- IR (MIR aka SWIR- short wave IR) and thermal IR wavelengths. 

  • Landsat Group 3 (LS 8-9): Equipped with an operational land imager (OLI) and Thermal IR Sensor (TIRS). The OLI augments the TM with the addition of a deep blue and a cirrus band, while TIRS adds a second thermal band. 

    • Landsat 9 is scheduled for launch in Sep, 2021. 

  • NASA builds and launches the satellites while the USGS Operates the Missions. 

__________________________________________________________________________________

---QGIS How To---

__________________________________________________________________________________

Image Classification

Image Classification: The process of sorting pixels into a finite number of individual classes, or categories of data, based on their spectral response (the measured brightness of a pixel across the image bands, as reflected by the pixel’s spectral signature).

  1. Supervised: Uses image pixels representing regions of known, homogenous surface composition -- training areas -- to classify unknown pixels. 

  1. Determine a classification scheme

  2. Create training areas

    1. Digitizing: Drawing Polygons around areas in an image. 

    2. Seeding: Grows areas based on spectral similarity to seed pixel. 

  3. Generate training area signatures

  4. Evaluate and refine signatures

  5. Assign pixels to classes using a classifier (a.k.a., “decision rule”).

    1. Parametric: Image is classified based on a statistical representation of the data derived from the training area signatures; all image pixels are classified. 

      1. Minimum Distance: Classifies pixels based on the spectral distance between the candidate pixel and the mean value of each signature (class) in each image band. 

      2. Maximum Likelihood: Classifies pixels based on the probability that a pixel falls within a certain class.

    2. Non-parametric: Pixels are classified as objects in feature space; only those pixels within the feature space object are classified. 

      1. Parallelepiped: The pixels values are compared to upper and lower limits of each signature class (i.e., the min/max pixel values in each band, or the mean of each band +/- 2 standard deviations)

      2. Feature Space: Classifies pixels that fall within non-parametric signatures identified in the feature space image not used very often because it is difficult to accurately create and evaluate non-parametric signatures.

  1. Unsupervised (aka clustering): Identifies and assigns groups of pixels exhibiting a similar spectral response a meaning (i.e. a land-cover categorization). 

  1. Determine a general classification scheme

  2. Group pixels into groups of similar values based on pixel- value relationships (i.e. Iterative ISODATA technique).

    1. Iterative Self-Organizing Data Analysis (ISODATA) Technique: Uses “spectral distance” between image pixels in feature space to classify pixels into a specified number of unique spectral groups (or “clusters”).

  3. Assign spectral classes to informational classes. 

__________________________________________________________________________________

Vegetation Indices

Vegetation Indexes (VI): Vegetation reflects strongly in the NIR and weakly in the SWIR portions of the EM Spectrum.

  • Normalized Difference Vegetation Index (NDVI): Quantifies vegetation by measuring the difference between NIR (which chlorophyl vegetation reflects) and red light (which vegetation absorbs); NDVI = (NIR - Red) / (NIR + Red)

    • NDVI is based on the observation that green plants reflect infrared light (to avoid overheating), but absorb visible light (to power photosynthesis). This gives green plants a light reflection signature that can distinguish them from objects that absorb all spectra of light or different spectra of light from plants.

    • NDVI ranges from -1 (unhealthy) to 1 (healthy). 

  • Delta/Differenced NDVI (dNDVI): dNDVI = preNDVI – postNDVI

  • Normalized Burn Ratio (NBR): Highlight burned areas and estimate fire severity; NBR= (NIR-SWIR)/(NIR+SWIR)

    • Delta/Differenced NBR (dNBR): Assesses the amount of vegetation and soil changes resulting from fire; dNBR = preNBR – postNBR

  • Normalized Difference Snow Index (NDSI): Estimates the amount of a pixel covered in snow; NDSI = (Green – SWIR) / (Green + SWIR)

  • Normalized Difference Water Index (NDWI): Delineates open water features, enhancing their presence in SRS; NDWI = (G - NIR) / (G + NIR)

  • Normalized Pigment Chlorophyll Ratio Index (NPCRI): Numerical indicator used to determine crop and/or vegetation clorphyll content; NPCRI = (Red – Blue) / (Red + Blue)

  • Normalized Difference Glacier Index (NDGI): Used to detect and monitor glaciers; NDGI = (NIR – Green) / (NIR + Green)

  • Normalized Different Moisture Index (NDMI): Used in combination with NDVI and/or ADVI to calculate vegetation moisture; NDMI = (NIR – SWIR) / (NIR + SWIR)

  • Advanced Vegetation Index (AVI): Similar to NDVI, monitors crop and forest variations over time; AVI = (NIR * (1 – Red) * (NIR – Red))1/3

  • Bare Soil Index (BSI): Combines Blue, Red, NIR, and SWIR to capture soil variations; BSI = ((Red+SWIR) - (NIR+Blue)) / ((Red+SWIR) + (NIR+Blue))

  • Topographic Wetness Index (TWI): TWI= In (upslope contributing area/tan(slope)) in dry- wet. 

__________________________________________________________________________________

Pre-Processing Steps

  1. Geometric Correction: Exact positioning of the image. 

    1. Georeferencing: Alignment of imagery to geographic location

    2. Orthorectifying: Correction for the effects of relief and view direction on pixel location.

  2. Absolute Radiometric Correction (aka Absolute correction): Account for sensor, solar, atmospheric, and topographic effects. 

    1. Conversion to Radiance

    2. Solar Correction: Accounts for solar influences on pixel values. Solar corrections converts at-sensor radiance to top of atmosphere (TOA) reflectance by incorporating exo-atmospheric solar irradiance (power of the sun), earth distance, and solar elevation angle. 

    3. Atmospheric Correction: Accounts for scattering and absorption in the Earth’s atmosphere. 

    4. Topographic Correction: Accounts for illumination effects from slope, aspect, and elevation. 

  3. Relative Radiometric Correction: Brings each band of an image to the same radiometric scale as the corresponding band of a referenced image to account for sensor, solar, and atmospheric differences.

__________________________________________________________________________________

Terrain Analysis

  1. Add a DEM

  2. Add Hillshade

    1. Raster- Analysis- Hillshade

      1. Input DEM, change Azimuth if you want, add an output layer in save file. 

  3. Change DEM Coloring

    1. Go to DEM Layer- Symbology- color ramp- elevation green-brown-white. Change Transparency (60%). 

  4. Slope

    1. Go to Raster- Analysis- Slope

      1. Input DEM Layer, Save as a slope file. 

      2. Change Symbology: Color Ramp- Green-Red-Yellow (shows high slopes). Change Transparency (60%) to make more visible. 

  5. Aspect

    1. Go to Raster- Analysis- Aspect

      1. Input DEM Layer, Return 0 for flight, save as aspect file. 

      2. Change Symbology: Single Band Single Pseudocolor- Color Ramp- 8 Classes: Inferno

_______________________________________________________________________

Adding Lat/Long to a Map

  1. Layer- Add Layer- Add delimited text layer. 

  2. Browse to the file with your data. Make sure that csv is selected for File format. Additionally, make sure that X field represents the column for your longitude points and Y field for latitude. QGIS is smart enough to recognize longitude and latitude columns but double check! 

  3. Select your CRS. 

_______________________________________________________________________

Creating a Multi-Band Image (Combining Spectral Data)

  1. Open Raster- Build Virtual Raster.

  2. Add all Bands and Run.

_______________________________________________________________________

Vector Geoprocessing 

  • Clip: Clipping Data to create a new shape file. 

  1. Add Input Layer, Clip Layer, and name output file. 



  • Area Calculation

  1. Open Field Calculator in attribute table in editing mode.

  2. Create a new field- name output field name- name outfield field type ($Area= sqm).

_______________________________________________________________________

Creating a Shapefile from a CSV

  1. Go to Layer- Add Layer- Add Delimited Text Layer. 

  2. Ensure Coordinates on CSV align with X-field and Y-field and set CRS. 

  3. Add to Map. 

  4. Export and save as a shape file. 

_______________________________________________________________________

Re-Project Data to Match CRS

  1. Go to vector menu- Data management tools- Reproject Layer. 

_______________________________________________________________________

---GEE How To---

_______________________________________________________________________

Google Earth Engine

  1. Open GEE.

  2. Add VI Code you want to conduct. 

  3. Search and pin a location.

  4. Search and pin a GIS Image Set.

  5. Run VI Program.

_______________________________________________________________________

---R---

_______________________________________________________________________

Adding a Country Shapefile

  1. Clear memory: rm(list=ls())

  2. Set Working Directory: setwd("~/Desktop/GIS Open Source/Activity 4")

  3. Install requisite packages:

    1. install.packages("sf")

    2. install.packages("terra")

    3. install.packages("plotwidgets")

    4. install.packages("ggshadow")

    5. install.packages("ggspatial")

    6. install.packages("ggnewscale")

    7. install.packages("janitor")

    8. install.packages("tidyverse")

    9. install.packages("rgdal")

  4. Load all packages

    1. library(raster)

    2. library(terra)

    3. library(sf)

    4. library(tidyverse)

    5. library(plotwidgets)

    6. library(ggshadow)

    7. library(ggspatial)

    8. library(ggnewscale)

    9. library(janitor)

    10. library(rgdal)

    11. install.packages("devtools")

  5. Load Country

    1. World_Countries_Generalized <- "World_Countries__Generalized_"

    2. limits <- readOGR(dsn = ".", layer = "World_Countries_Generalized")

    3. limits <- ne_countries(scale = 50, returnclass = "sf")

_______________________________________________________________________

Unsupervised Kmeans Classification

  1. Create a stack with all of our bands

  2. Clear memory: rm(list=ls())

  3. Add Bands

    1. blue = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B1.tif")

    2. green = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B2.tif")

    3. red = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B3.tif")

    4. nir = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B4.tif")

    5. swir1 = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B5.tif")

    6. swir2 = raster("LT05_L1TP_166039_19950929_20181025_01_T1_B7.tif")

    7. land5_stack = stack(blue, green, red, nir, swir1, swir2)

    8. land5_matrix <- getValues(land5_stack)

    9. i <- which(!is.na(land5_matrix))

    10. land5_matrix <- na.omit(land5_matrix)

    11. set.seed(99)

  4. # This code converts a sample to a matrix and this random sample will be used to train and test the data:

  5. land5_matrix2 <-land5_matrix[sample(nrow(land5_matrix), 500),]

  6. rfl5 = randomForest(land5_matrix2)

  7. rf_proxl8 <- randomForest(land8_matrix2,ntree = 1000, proximity = TRUE)$proximity

  8. Install requisite packages:

    1. install.packages("sp")

    2. install.packages("raster")

    3. install.packages("rgdal")

    4. install.packages("RColorBrewer")

    5. install.packages("ggplot2")

    6. library(raster)

    7. library(rgdal)

    8. library(ggplot2)

  9. ## The "5" in this case refers to the number of clusters you want to create

  10. l8_rf <- kmeans(rf_proxl8, 5, iter.max = 100, nstart = 10)

  11. land8_rf <- randomForest(land8_matrix2,as.factor(l8_rf$cluster),ntree = 500)

  12. l8_rf_raster<- predict(land8_stack,land8_rf)

  13. plot(l8_rf_raster)

  14. kmnclusterl5 <- kmeans(land5_matrix, centers = 5, iter.max = 100, nstart = 5, algorithm="Lloyd")

  15. kmeans_rasterl8 <- raster(land8_stack)

  16. kmeans_rasterl8[i] <- kmnclusterl8$cluster

  17. plot(kmeans_rasterl8)

  18. writeRaster(x = kmeans_rasterl8,

    1. filename="~/Desktop/GIS Open Source/Activity 8/kmeansl8.tif",

    2. format = "GTiff", # save as a tif

    3. overwrite = TRUE)  # OPTIONAL - be careful. This will OVERWRITE previous files.

_______________________________________________________________________

Load and View Multispectral Data

  1. ### Clear memory: rm(list=ls())

  2. # Install packages (if necessary) and load your package library

  3. # You NEED to load a package library every time

  4. # To check that a package is active, click on the Packages tab to the right and

  5. # sp automatically loads with the raster package

    1. #install.packages("sp")

    2. #install.packages("raster")

    3. #install.packages("rgdal")

    4. install.packages("RColorBrewer")

    5. install.packages("ggplot2")

    6. library(raster)

    7. library(rgdal)

    8. library(ggplot2)

  6. ##### There are three different Raster structures that you choose from depending on how many bands you need for a scene.

    1. ##  RasterLayer = one file, one band

    2. ##  RasterBrick = one file, multiple bands. Commonly used with RGB images

    3. ##  RasterStack = multiple files, multiple bands (Landsat GeoTIFFs)

  7. ####### IMPORT MULTISPECTRAL DATA #################

  8. ### Download data from https://earthexplorer.usgs.gov/

  9. ### Once you have downloaded the data, unzip it on your computer and save it to a working folder

  10. ## Remember you should see all of your individual bands

  11. ##### Set your working directory using one of the methods learned to the folder where your bands are locate

  12. ### To note, you can set your working directory multiple times in a session

  13. ## If you want to load your bands individually, you create a raster stack

  14. ## Set each band name as an object attached to its associated raster

    1. blue = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B2.tif")

    2. green = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B3.tif")

    3. red = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B4.tif")

    4. nir = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B5.tif")

    5. swir1 = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B6.tif")

    6. swir2 = raster("clipped_LC08_L1TP_172039_20190520_20190604_01_T1_B7.tif")

  15. ## Create an image using three bands and place them in the rgb channels

    1. # example r=swir1

    2. blue = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B2.tif")

    3. nir = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B5.tif")

    4. swir1 = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B6.tif")

    5. rgb = stack(swir1, nir, blue) #chosen to emphasize vegetation

  16. # Plot the image stack you just created using a linear stretch so the image isn't so washed out

  17. # Look at the raster package for this function to see other types of stretches

    1. plotRGB(rgb, stretch='lin')

  18. ## Though we are plotting a 3-band image, you can also plot each band individually

    1. plot(blue)

  19. e <- extent(229185, 345195, 4383585, 4472265)

  20. # Crop the red band based on the exten

    1. redcrop <- crop(red, e)

  21. # Plot the cropped data to make sure everything worked properly

    1. plot(redcrop)

  22. # Crop the near infrared band based on the extent

    1. nircrop <- crop(nir, e)

  23. # Plot the cropped data to make sure everything worked properly

    1. plot(nircrop)

_______________________________________________________________________

Perform Unsupervised Classification

#### METHOD 1: KMEANS ####

  1. ## This method clusters observations into a chosen number of categories based on the mean spectral signature.

    1. ## The algorithm is looking for clusters with similar spectral signatures (e.g. water)

  2. ### First, let's create a stack with all of our bands

    1. land8_stack = stack(blue, green, red, nir, swir1, swir2)

    2. blue = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B2.tif")

    3. green = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B2.tif")

    4. red = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B3.tif")

    5. nir = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B4.tif")

    6. swir1 = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B5.tif")

    7. swir2 = raster("Clipped_LS8_VirtualLC08_L1TP_172039_20210626_20210707_01_T1_B6.tif")

    8. land8_stack = stack(blue, green, red, nir, swir1, swir2)

  3. ### Before we can perform a classification, we need to convert the raster dataset to a matrix (array)

  4. # We will also remove any pixels from the analysis containing NA values

    1. land8_matrix <- getValues(land8_stack)

    2. i <- which(!is.na(land8_matrix))

    3. land8_matrix <- na.omit(land8_matrix)

  5. # It is important to set the seed generator because `kmeans` initiates the centers in random locations

    1. set.seed(99)

  6. # We want to create 5 clusters, allow 500 iterations, start with 5 random sets using "Lloyd" method

    1. ## Centers refers to the number of clusters you want to create

  7. kmnclusterl8 <- kmeans(land8_matrix, centers = 5, iter.max = 100, nstart = 5, algorithm="Lloyd")

  8. # We need to convert the kmeans cluster analysis back to a raster layer

    1. kmeans_rasterl8 <- raster(land8_stack)

    2. kmeans_rasterl8[i] <- kmnclusterl8$cluster

    3. plot(kmeans_rasterl8)

  9. writeRaster(x = kmeans_rasterl8,

    1.  filename="~/Desktop/GIS Open Source/Activity 8/kmeansl8.tif",

    2. format = "GTiff", # save as a tif

    3. overwrite = TRUE)  # OPTIONAL - be careful. This will OVERWRITE previous files.



#### METHOD 2: RANDOM FOREST ####

  1. install.packages("randomForest")

    1. library(randomForest)

  2. ## With this method, unique spectral signatures clusters are created using the kmeans method. These unique clusters are then used to train a random forest model to assign pixels to each cluster

  3. # This code converts a sample to a matrix and this random sample will be used to train and test the data

    1. land8_matrix2 <-land8_matrix[sample(nrow(land8_matrix), 500),]

    2. rfl8 = randomForest(land8_matrix2)

    3. rf_proxl8 <- randomForest(land8_matrix2,ntree = 1000, proximity = TRUE)$proximity

  4. ## The "5" in this case refers to the number of clusters you want to create

  5. l8_rf <- kmeans(rf_proxl8, 5, iter.max = 100, nstart = 10)

  6. land8_rf <- randomForest(land8_matrix2,as.factor(l8_rf$cluster),ntree = 500)

  7. l8_rf_raster<- predict(land8_stack,land8_rf)

  8. plot(l8_rf_raster)

  9. # To export your raster to load it to QGIS or have a permanent file on your computer

    1. writeRaster(x = l8_rf_raster,

      1. filename="~/Desktop/GIS Open Source/Activity 8/rfl8.tif",

      2. format = "GTiff", # save as a tif

      3. overwrite = TRUE)  # OPTIONAL - be careful. This will OVERWRITE previous files.



#### METHOD 3: CLUSTERING LARGE APPLICATIONS (CLARA) ####

  1. install.packages("cluster")

    1. library(cluster)

  2. clus_l8 <- clara(land8_matrix, 5, samples=500, metric="manhattan", pamLike=T)

  3. clara_l8raster <- raster(land8_stack)

  4. clara_l8raster[i] <- clus_l8$clustering

  5. plot(clara_l8raster)

  6. # To export your raster to load it to QGIS or have a permanent file on your computer

    1. writeRaster(x = clara_l8raster,

      1. filename="~/Desktop/GIS Open Source/Activity 8/claral8.tif",

      2. format = "GTiff", # save as a tif

      3. overwrite = TRUE)  # OPTIONAL - be careful. This will OVERWRITE previous files.

  7. ### First, let's create a stack with all of our bands

    1. blue = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B1.tif")

    2. green = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B2.tif")

    3. red = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B3.tif")

    4. nir = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B4.tif")

    5. swir1 = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B5.tif")

    6. swir2 = raster("Clipped_LS5LT05_L1TP_172039_19880530_20170717_01_T1_B7.tif")

    7. land5_stack = stack(blue, green, red, nir, swir1, swir2)

  8. ### Before we can perform a classification, we need to convert the raster dataset to a matrix (array)

  9. # We will also remove any pixels from the analysis containing NA values

    1. land5_matrix <- getValues(land5_stack)

    2. i <- which(!is.na(land5_matrix))

    3. land5_matrix <- na.omit(land5_matrix)

  10. # It is important to set the seed generator because `kmeans` initiates the centers in random locations

    1. set.seed(99)

  11. # We want to create 5 clusters, allow 500 iterations, start with 5 random sets using "Lloyd" method

    1. kmnclusterl5 <- kmeans(land5_matrix, centers = 5, iter.max = 50, nstart = 5, algorithm="Lloyd")

  12. # We need to convert the kmeans cluster analysis back to a raster layer

    1. kmeans_rasterl5 <- raster(land5_stack)

    2. kmeans_rasterl5[i] <- kmnclusterl5$cluster

    3. plot(kmeans_rasterl5)

  13. writeRaster(x = kmeans_rasterl5,

    1. filename="~/Desktop/GIS Open Source/Activity 8/kmeansl5.tif",

    2. format = "GTiff", # save as a tif

    3. overwrite = TRUE)  # OPTIONAL - be careful. This will OVERWRITE previous files.

_______________________________________________________________________

How to Graph the Temperature in Somalia

  1. #Step 1: Gather Temperature Data from World Clim

    1. taverage <- getData("worldclim", var = "tmean", res = 2.5)

  2. #Step 2: Discern the country code for Somalia

    1. x <- getData("ISO3")

    2. x[x[, "NAME"] == "Somalia", ]

  3. #Step 3: Download the country vector data.

    1. somalia <- getData("GADM", country = "SOM", level = 1)

  4. #Step 4: Plot your country to make sure it worked correctly. 

    1. plot(somalia)

  5. #Step 5: Download global temperature data. 

    1. taverage <- getData("worldclim", var = "tmean", res = 2.5)

  6. #Step 6: Plot taverage to ensure you have it downloaded. 

    1. plot(taverage)

  7. #Step 7: Add month names to your taverage plot. 

    1. names(taverage) <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

  8. #Step 8: Replot to ensure the country names were added to the taverage data. 

    1. plot(taverage)

  9. #Step 9: Crop the taverage data to match your country. 

    1. taverage_crop <- crop(taverage, somalia)

  10. #Step 10: Plot the cropped temp data for your country. 

    1. plot(taverage_crop)

  11. #Step 11: Add a Mask to remove nearby countries and replot. 

    1. somalia_taverage <- mask(crop(taverage, somalia), somalia)

    2. plot(somalia_taverage)

  12. #Step 12: Extract the average temperature across your country: 

    1. mean_temperature <- cellStats(taverage_crop_mask, stat = "mean")

  13. #Step 13: Graph/plot the mean temperature by month: 

    1. mean_temperature

  14. #Step 14: Create a barplot of your monthly data (use /10 of data to get correct C)

    1. barplot((mean_temperature/10), main = "Mean temperature across Somalia", ylab = "temperature (C)", col = "blue", axis.lty = 1)

_______________________________________________________________________

How to View DEM (Raster) Slope, Aspect

  1. Clear memory: rm(list=ls())

  2. Set your working directory: 'Session' > 'Set working directory' > 'Choose directory'. 

    1. setwd("~/Desktop/GIS Open Source/R_Activity3")

  3. Load Packages Raster, sp, and rgdal. If not loaded, install. 

    1. #Potential: #install.packages("rgdal") and library(rgdal)

    2. install.packages("rgdal")

    3. library(rgdal)

  4. Import clipped Kilimanjaro DEM Data as kil_dem from what you saved it as from QGIS. 

    1. kil_dem <- raster("DEM_Clipped_KiliForest.tif")

  5. View a summary of the kil_dem data by running the raster object name.

    1. kil_dem

  6. To view this data as a colorized elevation plot:

    1. plot(kil_dem)

  7. Create a terrain function for the slope. 

    1. # Terrain Help: ?terrain

    2. #terrain formula: terrain(x, opt="slope", unit="radians", neighbors=8, filename="")

    3. # x = your dem data

    4. # opt = whichever topographical computation you would like to perform

    5. # units = whichever character type you'd like and in this case degrees

    6. # neighbors = the neighboring cells used to calculate the value of a cell

    7. terrain(kil_dem, opt= "slope", unit= "radians", neighbors= 8, filename= "kil_slope"

  8. Add the slope to the Environmental Database. 

    1. kil_slope <- raster("kil_slope")

  9. Plot the Slope: plot(kil_slope)

  10. Repeat steps 7-9 for aspect. 

    1. terrain(kil_dem, opt= "aspect", unit= "radians", neighbors= 8, filename= "kil_aspect")

    2. kil_aspect <- raster("kil_aspect")

    3. plot(kil_aspect)

  11. To get object summary data, type the data name and press enter. 

    1. kil_aspect

    2. kil_dem

    3. kil_slope

  12. To compute hillshade, both slope and aspect need to be in radians. Repeat steps #7-9 for hillshade. 

    1. #hillshade help: ?hillshade

    2. #hillshade code: hillShade(slope, aspect, angle=45, direction=0, filename='', normalize=FALSE, ...)

    3. # angle = elevation angle of the light source (the sun)

    4. # direction = azimuth angle of the light source (the sun)

    5. hillShade(kil_slope, kil_aspect, angle=45, direction=300, filename= "kil_hillshade")

    6. kil_hillshade <- raster("kil_hillshade")

    7. plot(kil_hillshade)

  13. Remove the legend, change the color scale, and add a title to the map. 

    1. plot(kil_hillshade, col=grey(0:100/100), legend=FALSE, main='Mount Kilimanjaro')

  14. Use writeRaster function so you can load data you create into QGIS. 

    1. #writeraster code: writeRaster(file you want to export, filename="output file name", format="GTiff", overwrite=TRUE)

    2. writeRaster(kil_hillshade, filename="kil_hillshade_export", format="GTiff", overwrite=TRUE)

    3. kil_hillshade_export

  15. Ensure your writeRaster goes to the correct folder by resetting your working folder and re-exporting. 

    1. setwd("~/Desktop/GIS Open Source/R_Activity3")

  16. Export the DEM of your choice using writerRaster. 

    1. writeRaster(kil_slope, filename="kil_slope_export", format="GTiff", overwrite=TRUE)

_______________________________________________________________________

Visualizing Landsat Data

  1. Go to Earth Explorer and download data. 

  2. Recent Data: Use Landsat Archive, Collection 1, Level- 1 and then Landsat 8 (L8 OLI/TIRS). 

  3. DL Landsat Look Images with Geographic Reference for color visible light imagery. 

  4. DL Level-1 GeoTIFF data if using other bands, including IR needed to calculate NDVI. Distributed as .tar.gz files

  5. Raster in R

  6. Layer(): Raster Layer= one file, one band. 

  7. Brick(): Raster Brick= one file, multiple bands. Common with RGB images. 

__________________________________________________________________________________

Terminology

  • Analysis Ready Data (ARD): USGS project to provide data in an alternative format. 

  • Aspect: Orientation of Slope in Cardinal Directions

  • Coordinate Reference System (CRS): Association of numerical coordinates with a position on the Earth’s surface. CRS’ are identified by a unique European Petroleum Search Group (EPSG) Code.  

    • Geographic CRS: Lat/Long in degrees, minutes, seconds and based on spherical modeling with unique coordinates anywhere on Earth. 

      • EPSG: 4326; WGS-84 LAT/LONG, the default in QGIS. 

    • Projected CRS: Easting and Northings in units of distance with projection of Earth onto a flat plane (Mercator). 

      • EPSG: 3857; Standard Web mapping

  • Digital Elevation Models (DEM): An elevation model created from Ground Surveying with DGPS, Stereo Photogrammetry, Contours, LIDAR, and/or RADAR Interferometry. Used to determine areas, delineate drainage areas, identify slopes, aspects, geologic structures, viewsheds, orthorectification, 3D simulations, change analysis, &c. 

    • Digital Terrain Model (DTM): Digital Terrain Model

    • Digital Surface Model (DSM): A DTM with additional features. 

  • Digitizing: Drawing Polygons around areas in an image. 

  • Essential Climate Variables (ECVs): Land Surface temperature, burned area extent, dynamic surface water, and snow cover areas. 

  • Filetypes

    • .shp: ESRI Shape files; Stores features geometry. 

    • .dbf: Stores a features attribute’s. 

    • .prj: Stores information about data’s projection

    • .xml: Extensible markup Language used to describe data.  

    • .sbx: Spatial Index.

    • .sbn: Spatial Index.

    • .dxf: drawing exchange format; enables data interoperability between AutoCAD and other programs. 

    • .dwg: Binary format for 2D and 3D design data. 

    • .tab: Proprietary vector format; the table structure is the main file for a MapInfo table, and the link between all other files and holds information about the type fo data file. 

    • .gpx: Used for waypoints, routes, and tracks for uploading to GPS. 

  • Image: Collection of spatially arranged measurements (bands) captured in the scene at a given time. 

  • Image Classification: The process of sorting pixels into a finite number of individual classes, or categories of data, based on their spectral response (the measured brightness of a pixel across the image bands, as reflected by the pixel’s spectral signature).

    • Supervised: Uses image pixels representing regions of known, homogenous surface composition -- training areas -- to classify unknown pixels. 

    • Unsupervised (aka clustering): Identifies and assigns groups of pixels exhibiting a similar spectral response a meaning (i.e. a land-cover categorization). 

  • Metadata: Information about the dataset; who made it? How did they make it? What does the data include? When was the data put together? Where does the data cover?

  • Orthorectification: Georeferencing of aerial photos

  • Parametric: Image classification based on a statistical representation of the data derived from the training area signatures; all image pixels are classified. 

    • Minimum Distance: Classifies pixels based on the spectral distance between the candidate pixel and the mean value of each signature (class) in each image band. 

    • Maximum Likelihood: Classifies pixels based on the probability that a pixel falls within a certain class.

  • Non-parametric: Pixels are classified as objects in feature space; only those pixels within the feature space object are classified. 

    • Parallelepiped: The pixels values are compared to upper and lower limits of each signature class (i.e., the min/max pixel values in each band, or the mean of each band +/- 2 standard deviations)

    • Feature Space: Classifies pixels that fall within non-parametric signatures identified in the feature space image not used very often because it is difficult to accurately create and evaluate non-parametric signatures.

  • Post Classification Change Detection: Classification of images from multiple dates and subsequent comparison of maps to identify changes. 

  • Raster Data: Grids where each cell has a common layer; pixels. 

  • Scene: Extent or footprint that exists on the ground. 

  • Seeding: Grows areas based on spectral similarity to seed pixel. 

  • Unsupervised Classification

    • KMEANS: Clusters observations into a chosen number of categories based on the mean spectral signature. 

    • Random Forest: Unique spectral signatures clusters are created using the kmeans method. These unique clusters are then used to train a random forest model to assign pixels to each cluster. 

    • Clustering Large Applications (CLARA).

  • Vector Data: Points, Lines, Areas, Polygons; shapes. 

__________________________________________________________________________________

Resources

__________________________________________________________________________________

Chronology

  • 27 Sep, 2021- ~2036: Landsat 9 Mission (NASA).

  • Jun, 2015: The European Space Agencies (ESA) Copernicus Program begins delivering Sentinel 2 data at multi-spectral 10m spatial resolution (NASA). 

  • 2013- Present: Landsat 8 Mission (NASA). 

  • 1999- Present: Landsat 7 Mission (NASA).

  • 1984- 2013: Landsat 5 Mission (NASA). 

  • 1982- 1993: Landsat 4 Mission (NASA).

  • 1978- 1983: Landsat 3 Mission (NASA).

  • 1975- 1982: Landsat 2 Mission (NASA).

  • 1972- 1978: Landsat 1 Mission (NASA).

  • 1972: NASA launches the first Landsat satellite (NASA).

__________________________________________________________________________________