update to use terra package, calculate hillshade in the script

NEONScience · Feb 17, 2024 · bd9c103 · bd9c103
1 parent d90f780
commit bd9c103
Show file tree

Hide file tree

Showing 16 changed files with 359 additions and 716 deletions.
diff --git a/tutorials/R/Geospatial-skills/intro-raster-r/01-Plot-Raster/01-Plot-Raster.R b/tutorials/R/Geospatial-skills/intro-raster-r/01-Plot-Raster/01-Plot-Raster.R
@@ -1,31 +1,29 @@
-## ----load-libraries----------------------------------------------------
+## ----load-libraries----------------------------------------------------------------------------------------------------------------------------------------------------------------
 # if they are not already loaded
-library(rgdal)
-library(raster)
+library(terra)
 
-# set working directory to ensure R can find the file we wish to import
-wd <- "~/Git/data/" # this will depend on your local environment environment
-# be sure that the downloaded file is in this directory
+# set working directory
+wd <- "~/data/" # this will depend on your local environment environment
 setwd(wd)
 
-# import raster
-DSM_HARV <- raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif"))
+# import raster into R
+dsm_harv_file <- paste0(wd, "DP3.30024.001/neon-aop-products/2022/FullSite/D01/2022_HARV_7/L3/DiscreteLidar/DSMGtif/NEON_D01_HARV_DP3_732000_4713000_DSM.tif")
+DSM_HARV <- rast(dsm_harv_file)
 
 
-
-## ----hist-raster, fig.cap="Digital surface model showing the continuous elevation of NEON's site Harvard Forest"----
+## ----hist-raster, fig.cap="Digital surface model showing the continuous elevation of NEON's site Harvard Forest"-------------------------------------------------------------------
 # Plot raster object
 plot(DSM_HARV,
      main="Digital Surface Model\nNEON Harvard Forest Field Site")
 
 
 
-## ----create-histogram-breaks, fig.cap="Histogram of digital surface model showing the distribution of the elevation of NEON's site Harvard Forest"----
+## ----create-histogram-breaks, fig.cap="Histogram of digital surface model showing the distribution of the elevation of NEON's site Harvard Forest"---------------------------------
 # Plot distribution of raster values 
 DSMhist<-hist(DSM_HARV,
      breaks=3,
      main="Histogram Digital Surface Model\n NEON Harvard Forest Field Site",
-     col="wheat3",  # changes bin color
+     col="lightblue",  # changes bin color
      xlab= "Elevation (m)")  # label the x-axis
 
 # Where are breaks and how many pixels in each category?
@@ -34,131 +32,57 @@ DSMhist$counts
 
 
 
-## ----plot-with-breaks, fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with three breaks"----
+## ----plot-with-breaks, fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with three breaks"-------------------------------------------------------
 # plot using breaks.
 plot(DSM_HARV, 
      breaks = c(300, 350, 400, 450), 
      col = terrain.colors(3),
-     main="Digital Surface Model (DSM)\n NEON Harvard Forest Field Site")
+     main="Digital Surface Model (DSM) - HARV")
 
 
 
-## ----add-plot-title, fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with UTM Westing Coordinate (m) on the x-axis and UTM Northing Coordinate (m) on the y-axis"----
+## ----add-plot-title, fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with UTM Westing Coordinate (m) on the x-axis and UTM Northing Coordinate (m) on the y-axis", fig.height=6----
 # Assign color to a object for repeat use/ ease of changing
 myCol = terrain.colors(3)
 
 # Add axis labels
 plot(DSM_HARV, 
      breaks = c(300, 350, 400, 450), 
      col = myCol,
-     main="Digital Surface Model\nNEON Harvard Forest Field Site", 
-     xlab = "UTM Westing Coordinate (m)", 
-     ylab = "UTM Northing Coordinate (m)")
+     main="Digital Surface Model - HARV", 
+     xlab = "UTM Easting (m)", 
+     ylab = "UTM Northing (m)")
 
 
-## ----turn-off-axes,fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with no axes"----
+## ----turn-off-axes,fig.cap="Digital surface model showing the elevation of NEON's site Harvard Forest with no axes"----------------------------------------------------------------
 # or we can turn off the axis altogether
 plot(DSM_HARV, 
      breaks = c(300, 350, 400, 450), 
      col = myCol,
-     main="Digital Surface Model\n NEON Harvard Forest Field Site", 
+     main="Digital Surface Model - HARV", 
      axes=FALSE)
 
 
 
-## ----challenge-code-plotting, include=TRUE, results="hide", echo=FALSE----
+## ----challenge-code-plotting, include=TRUE, results="hide", echo=FALSE-------------------------------------------------------------------------------------------------------------
 # Find min & max
-DSM_HARV@data
-
-# Pixel range & even category width
-(416.07-305.07)/6
-
-# Break every 18.5m starting at 305.07
+min(DSM_HARV)
 
+# Break every 10m starting at 310
 # Plot with 6 categories at even intervals across the pixel value range. 
 plot(DSM_HARV, 
-     #breaks = c(305, 323.5, 342, 360.5, 379, 397.5, 417),  #manual entry
-     breaks = seq(305, 417, by=18.5),  #define start & end, and interval
-     col = terrain.colors (6),
-      main="Digital Surface Model\nNEON Harvard Forest Field Site", 
-     xlab = "UTM Westing Coordinate", 
-     ylab = "UTM Northing Coordinate")
-
-
-
-## ----hillshade, fig.cap="Hillshade digital surface model showing the elevation of NEON's site Harvard Forest in grayscale"----
-# import DSM hillshade
-DSM_hill_HARV <- 
-  raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif"))
+     breaks = seq(310, 440, by=10),  #define start & end, and interval
+     col = terrain.colors(13),
+      main="Digital Surface Model - HARV", 
+     xlab = "Easting", 
+     ylab = "Northing")
 
-# plot hillshade using a grayscale color ramp that looks like shadows.
-plot(DSM_hill_HARV,
-    col=grey(1:100/100),  # create a color ramp of grey colors
-    legend=FALSE,
-    main="Hillshade - DSM\n NEON Harvard Forest Field Site",
-    axes=FALSE)
 
 
-
-## ----overlay-hillshade, fig.cap="Digital surface model overlaying the hillshade raster showing the 3D elevation of NEON's site Harvard Forest"----
-
-# plot hillshade using a grayscale color ramp that looks like shadows.
-plot(DSM_hill_HARV,
-    col=grey(1:100/100),  #create a color ramp of grey colors
-    legend=F,
-    main="DSM with Hillshade \n NEON Harvard Forest Field Site",
-    axes=FALSE)
-
-# add the DSM on top of the hillshade
-plot(DSM_HARV,
-     col=rainbow(100),
-     alpha=0.4,
-     add=T,
-     legend=F)
-
-
-## ----challenge-hillshade-layering, fig.cap=c("Digital surface model overlaying the hillshade raster showing the 3D elevation of NEON's site San Joaquin Experiment Range","Digital terrain model overlaying the hillshade raster showing the 3D ground surface of NEON's site San Joaquin Experiment Range"), echo=FALSE----
-# CREATE DSM MAPS
-# import DSM 
-DSM_SJER <- raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMcrop.tif"))
-# import DSM hillshade
-DSM_hill_SJER <- 
-  raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill.tif"))
-
-# plot hillshade using a grayscale color ramp that looks like shadows.
-plot(DSM_hill_SJER,
-    col=grey(1:100/100),  #create a color ramp of grey colors
-    legend=F,
-    main="DSM with Hillshade\n NEON SJER Field Site",
-    xlab = "UTM Westing Coordinate", 
-    ylab = "UTM Northing Coordinate")
-
-# add the DSM on top of the hillshade
-plot(DSM_SJER,
-     col=terrain.colors(100),
-     alpha=0.7,
-     add=T,
-     legend=F)
-
-# CREATE SJER DTM MAP
-# import DTM 
-DTM_SJER <- raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_DTMcrop.tif"))
-# import DTM hillshade
-DTM_hill_SJER <- 
-  raster(paste0(wd,"NEON-DS-Airborne-Remote-Sensing/SJER/DTM/SJER_DTMhill.tif"))
-
-# plot hillshade using a grayscale color ramp that looks like shadows.
-plot(DTM_hill_SJER,
-    col=grey(1:100/100),  #create a color ramp of grey colors
-    legend=F,
-    main="DTM with Hillshade\n NEON SJER Field Site",
-    axes=F)
-
-# add the DSM on top of the hillshade
-plot(DTM_SJER,
-     col=terrain.colors(100),
-     alpha=0.4,
-     add=T,
-     legend=F)
-
+## ----slope-aspect-hill-------------------------------------------------------------------------------------------------------------------------------------------------------------
+slope <- terrain(DSM_HARV, "slope", unit="radians")
+aspect <- terrain(DSM_HARV, "aspect", unit="radians")
+hill <- shade(slope, aspect, 45, 270)
+plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
+plot(DSM_HARV, col=terrain.colors(25, alpha=0.35), add=TRUE, main="HARV DSM with Hillshade")