r - 如何使用 R 找到具有不同土地利用分类的栅格的平均坡度？

Question

我是 R 中的 GIS 新手。我有两个输入栅格(高程 DEM 和土地利用分类栅格)，我的任务是重新投影、聚合为路线分辨率，最后计算每 block 土地的平均土地坡度使用分类。我无法找到土地用途的平均坡度；这就是我需要帮助的地方。

高程栅格DEM可下载here ，而土地利用分类栅格可以下载 here .

这是我所做的:

#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)

调用freq(mask)命令显示有17个土地利用值(不包括NA值)。输出如下所示:

      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

如何找到每个土地使用值的平均坡度？我认为解决方案涉及组合栅格并运行 cellStats 操作，但我不完全确定。

Answer 1

我认为最好计算高分辨率数据的斜率。如果必须的话，您可以汇总这些值。如果您根据聚合(因此平滑的)高程数据计算坡度，您可能会严重低估它。此外，只要有可能，您就应该避免投影栅格数据，因为这会导致数据质量损失。在这种情况下，您必须对两个数据源之一执行此操作，但不应同时对两个数据源执行此操作。当然，不要对同一源的数据进行投影和重新采样(然后数据质量就会恶化两次)。

在这种情况下，我会将土地覆盖数据重新采样为坡度数据。土地覆盖数据的空间分辨率为30 m，高程数据的空间分辨率为~9 m。分解土地覆盖数据不会对数据质量产生太大影响，而聚合坡度数据则会影响数据质量。

这是我可能会做的:

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

这在我的笔记本电脑上运行良好。这是一种替代方法，强度较低但概念上较差，对斜率数据进行重新采样:

# 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

好消息是差异非常小。