Using R, how can I find the average slope of a raster with land use area information?

Question

I am new to GIS in R. I have two input rasters (an elevation DEM and a land use classified raster), and I am tasked with reprojecting, aggregating into courser resolution, and finally calculating the average land slope for each land use classification. I am having trouble finding the average slope for the land uses; this is what I need help with.

For ease assisting with this question, I have provided both of the rasters on a Google Drive for users to download. The elevation raster DEM can be downloaded here, whereas the land use classification raster can be downloaded here.

Here is what I have done so far:

#Load libraries
library(raster)
library(sf)

#Load data
cdl19 <- raster("CDL_2019_clip_20220329135051_1368678123.tif")
elev <- raster("elevation.tif")

#View information about the rasters
plot(cdl19)
hist(cdl19)

plot(elev)
hist(elev)

#Check coordinate reference systems
crs(cdl19)
crs(elev)

#Reproject each raster so that they will be in the same grid system
cor_ref = '+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

cdl19_proj <- projectRaster(cdl19, crs = cor_ref, method = 'ngb') #Use 'ngb' for categorical data

elev_proj <- projectRaster(elev, crs = cor_ref, method = 'bilinear') #use 'bilinear' for continuous data

#The references should now be the same
#####################################################

#Resample to eliminate resolution problems
elev_resamp <- resample(elev_proj, cdl19_proj, method = 'bilinear')

#Aggregate to larger, 1km x 1km cell size
#Note: This is not exactly 1km x 1km, but projecting to a finer coordinate system would overwhelm my computer's memory capabilities.

cdl19_1k <- aggregate(cdl19_proj, fact = 33.333333333, fun = modal, na.rm=TRUE)

elev_1k <- aggregate(elev_resamp, fact = 33.333333333, fun = modal, na.rm=TRUE)


#Calculate slope of DEM
slope <- terrain(elev_1k, opt = 'slope', units = 'degrees')

#Crop and mask categorical raster to the DEM
cdl19_1k_crop <- crop(cdl19_1k, extent(slope))
mask <- mask(cdl19_1k_crop, mask = slope)

Here is where I need assistance. Calling the freq(mask) shows that there are 17 land use values (not including NA values). The output can be seen below:

      value count
 [1,]     1   137
 [2,]     4     9
 [3,]     5   230
 [4,]    24    16
 [5,]    28     1
 [6,]    36     7
 [7,]    37     1
 [8,]    61    13
 [9,]   111   136
[10,]   121     7
[11,]   122    52
[12,]   123     7
[13,]   124     2
[14,]   141   351
[15,]   143     1
[16,]   176  2021
[17,]   195    22
[18,]    NA  1342

How can I find the average slope of each land use value? I would presume that the solution involves combining the rasters and running a cellStats operation, but I'm not entirely sure.

Answer 1

I would point out that slope, like aspect, is an angle. Taking the mean of an angle is not as simple as mean(slope). To get as mean angle generally, you take the atan2 of the mean cos and sin in radians.

Here is a function to calculate the angle mean (in development version of spatialEco). The default input is degree but, you can specify radians if your slope or aspect are output in radians.

mean_angle <- function(a, angle=c("degree", "radians")) {
  angle=angle[1]
  deg2rad <- function(x) { x * pi/180} 
  rad2deg <- function(x) { x * 180/pi }
  deg2vec <- function(x, ang = c("degree", "radians")) { 
    if(ang == "degree") {
      a <- c(sin(deg2rad(x)), cos(deg2rad(x)))
    } else if(ang == "radians") {
      a <- c(sin(x), cos(x))
    }
    return(a)
  }
  vec2deg <- function(x) {
    res <- rad2deg(atan2(x[1], x[2]))
      if (res < 0) { res <- 360 + res }
    return(res)
  }
  mean_vec <- function(x) {
    y <- lapply(x, deg2vec, ang=angle)
    Reduce(`+`, y)/length(y)
  }
  return( vec2deg(mean_vec(a)) ) 
}

When we compare the mean to the angle mean you can see that in one case the mean is the same but with a different distribution the angle mean is different.

# mean angle = 95 and mean = 95
mean_angle(c(180, 10))
  mean(c(180, 10))

# mean angle = 93.22 and mean = 100
mean_angle(c(90, 180, 70, 60))
  mean(c(90, 180, 70, 60))

For your application you would simply replace the function call in:

zonal(x=slope, z=polys, fun=mean_angle)

Answer 2

You are really close (or maybe already you did it), you just have to calculate slope out of your elevation data and then use the zonal statistics function:

that is summarized values of a Raster* object for each "zone" defined by a RasterLayer.

# obtain slope
x_slope = terrain(elev_resamp, v = "slope")

# calculate zonal statistics
zstats = zonal(x = x_slope, z = cdl19_proj, fun = "mean")

zstats
#      zone        mean
# [1,]    0 0.058269462
# [2,]    1 0.025891800
# [3,]    2 0.036698524
# [4,]    4 0.028648946
# [5,]    5 0.026401503