我如何找到使用R不同土地使用分类的栅格的平均斜率？
我是 R GIS 新手。我有两个输入栅格（高程 DEM 和土地利用分类栅格），我的任务是重新投影、聚合成路线分辨率，最后计算每个土地利用分类的平均土地坡度。我无法找到土地用途的平均坡度；这就是我需要帮助的地方。
高程栅格 DEM 可以在此处下载，而土地利用分类栅格可以下载此处。
这是我所做的：
#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 操作，但我不完全确定。
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.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（1）
我认为最好在高分辨率数据上计算斜率。如果必须，则可以汇总这些值。如果您在汇总（因此平滑）高程数据上计算斜率，则可能会严重低估它。另外，只要可能，您应该避免投影栅格数据，因为这会导致数据质量丢失。在这种情况下，您必须为两个数据源之一做到这一点，但是您不应该为这两个数据源做到这一点。当然，请勿投射和将数据重新采样相同的源（然后您两次降低数据质量）。
在这种情况下，我将将土地覆盖数据重新置于坡度数据。土地覆盖数据的空间分辨率为30 m，对于高程数据约为9 m。分解土地覆盖数据不会对数据质量产生很大影响，而汇总坡度数据可能会影响数据质量。
这是我可能会做的：
这在笔记本电脑上运行良好。这是一种替代，不太密集但在概念上较低的方法，可以重新示例斜率数据：
好消息是差异很小。
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:
This runs fine on my laptop. Here is an alternative, less intensive but conceptually inferior, approach that resamples the slope data:
The good news is that the difference is very small.