如何使用 R 找到具有不同土地利用分类的栅格的平均坡度？
How can I find the average slope of a raster with different land use classifications using R?
我是 R 中的 GIS 新手。我有两个输入栅格（一个高程 DEM 和一个土地利用分类栅格），我的任务是重新投影，聚合成更粗的分辨率，最后计算平均土地坡度每个土地利用分类。我无法找到土地用途的平均坡度；这就是我需要帮助的。
可以下载高程栅格 DEM here, whereas the land use classification raster can be downloaded 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 操作，但我不太确定。
我认为最好在高分辨率数据上计算斜率。如果必须，您可以聚合这些值。如果您根据聚合（因此平滑）的高程数据计算坡度，您可能会严重低估它。此外，应尽可能避免投影栅格数据，因为这会导致数据质量下降。在这种情况下，您必须对两个数据源之一执行此操作，但不应该对两个数据源都执行此操作。当然不要投影 和 对相同来源的数据重新采样（然后你会使数据质量恶化两次）。
在这种情况下，我会将土地覆盖数据重新采样为坡度数据。土地覆盖数据的空间分辨率为 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
好消息是差别很小。
我是 R 中的 GIS 新手。我有两个输入栅格（一个高程 DEM 和一个土地利用分类栅格），我的任务是重新投影，聚合成更粗的分辨率，最后计算平均土地坡度每个土地利用分类。我无法找到土地用途的平均坡度；这就是我需要帮助的。
可以下载高程栅格 DEM here, whereas the land use classification raster can be downloaded 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 操作，但我不太确定。
我认为最好在高分辨率数据上计算斜率。如果必须，您可以聚合这些值。如果您根据聚合（因此平滑）的高程数据计算坡度，您可能会严重低估它。此外，应尽可能避免投影栅格数据，因为这会导致数据质量下降。在这种情况下，您必须对两个数据源之一执行此操作，但不应该对两个数据源都执行此操作。当然不要投影 和 对相同来源的数据重新采样（然后你会使数据质量恶化两次）。
在这种情况下，我会将土地覆盖数据重新采样为坡度数据。土地覆盖数据的空间分辨率为 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
好消息是差别很小。