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.
The elevation raster DEM can be downloaded here, whereas the land use classification raster can be downloaded here.
Here is what I have done:
#Load libraries
library(raster)
#Load data
cdl19 <- raster("CDL_2019_clip_20220329135051_1368678123.tif")
elev <- raster("elevation.tif")
#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
#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)
Calling the freq(mask) command 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 think that the solution involves combining the rasters and running a cellStats operation, but I'm not entirely sure.
I think it would be best to compute slope on the high resolution data. If you must, you can then aggregate these values. If you compute slope on aggregated (hence smoothed) elevation data you may grossly underestimate it. Also, whenever possible you should avoid projecting raster data as that leads to data quality loss. In this case you have to do that for one of the two data sources, but you should not do it for both. Certainly do not project and resample the data same source (then you deteriorate the data quality twice).
In this case, I would resample the land cover data to the slope data. The land cover data has a spatial resolution of 30 m, and it is ~9 m for the elevation data. Disaggregating the land cover data does not affect the data quality much, whereas aggregating the slope data could.
Here is what I might do:
library(terra)
cdl19 <- rast("CDL_2019_clip_20220329135051_1368678123.tif")
names(cdl19) <- "cdl"
elev <- rast("elevation.tif")
# to speed things up, remove elevation data not overlapping with cdl
e <- crop(elev, c(-97, -96.35, 39.08, 39.6))
s <- terrain(e, "slope")
cd <- project(cdl19, s, "near")
z <- zonal(cd, s, mean)
head(z)
# cdl slope
#1 0 4.065172
#2 1 1.734034
#3 2 2.064955
#4 4 1.872033
#5 5 1.717542
#6 6 1.155761
This runs fine on my laptop. Here is an alternative, less intensive but conceptually inferior, approach that resamples the slope data:
# first aggregate a little to get just below the output resolution
sa <- aggregate(s, 2, mean)
scd <- project(sa, cdl19)
zz <- zonal(scd, cdl19, mean)
head(zz)
# cdl slope
#1 0 4.072662
#2 1 1.758178
#3 2 2.085544
#4 4 1.885815
#5 5 1.739788
#6 6 1.169922
The good news is that the difference is very small.