绘制一个区域的斜率并返回高于和低于 R 中阈值的百分比

Question

我正在尝试计算坡度为 0，+/- 5 度的区域的比例。另一种说法是高于 5 度和低于 5 度都是不好的。我正在尝试查找实际数字和图形。

为了实现这一点，我转向了 R 并使用了 Raster 包。让我们使用一个通用的国家/地区，在本例中为菲律宾

{list.of.packages <- c("sp","raster","rasterVis","maptools","rgeos")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)}

library(sp)  # classes for spatial data
library(raster)  # grids, rasters
library(rasterVis)  # raster visualisation
library(maptools)
library(rgeos)

现在让我们获取高度信息并绘制斜率。

elevation <- getData("alt", country = "PHL")
x <- terrain(elevation, opt = c("slope", "aspect"), unit = "degrees")
plot(x$slope)

由于规模不大，所以我们简单看一下巴拉望岛

e <- drawExtent(show=TRUE) #to crop out Palawan (it's the long skinny island that is roughly midway on the left and is oriented between 2 and 8 O'clock)
gewataSub <- crop(x,e)
plot(gewataSub, 1)## Now visualize the new cropped object

可视化效果更好一些。我感觉到斜坡的大小，并且在 5 度的限制下，我大部分时间都被限制在海岸上。但我还需要多一点分析。

我希望结果分为两部分： 1.“所选区域的 35%（构成）坡度超过 +/- 5 度”或“所选区域的 65% 在 +/- 5 度以内”。（使用代码获取它） 2. +/- 5 度以内的所有东西都是一种颜色的图片，称它为好或绿色，而其他所有东西都是另一种颜色，称其为坏或红色。

谢谢

Answer 1

您可以使用 raster 包中的 reclassify 来实现。该函数为位于定义间隔内的每个单元格值分配一个特定值。例如，您可以将区间 (0,5] 内的单元格值分配给值 0，将区间 (5, maxSlope] 内的单元格值分配给值 1.

library(raster)  
library(rasterVis)  

elevation <- getData("alt", country = "PHL")
x <- terrain(elevation, opt = c("slope", "aspect"), unit = "degrees")
plot(x$slope)

e <- drawExtent(show = TRUE)
gewataSub <- crop(x, e)
plot(gewataSub$slope, 1)

m <- c(0, 5, 0,  5, maxValue(gewataSub$slope), 1)
rclmat <- matrix(m, ncol = 3, byrow = TRUE)
rc <- reclassify(gewataSub$slope, rclmat)

levelplot(
  rc,
  margin = F,
  col.regions = c("wheat", "gray"),
  colorkey = list(at = c(0, 1, 2), labels = list(at = c(0.5, 1.5), labels = c("<= 5", "> 5")))
)

重新分类后您可以计算百分比：

length(rc[rc == 0]) / (length(rc[rc == 0]) + length(rc[rc == 1])) # <= 5 degrees
[1] 0.6628788
length(rc[rc == 1]) / (length(rc[rc == 0]) + length(rc[rc == 1])) # > 5 degrees
[1] 0.3371212

Answer 2

没有负斜率，所以我假设您想要那些小于 5 度的斜率

library(raster)
elevation <- getData('alt', country='CHE')
x <- terrain(elevation, opt='slope', unit='degrees')

z <- x <= 5

现在您可以使用 freq

对细胞进行计数

f <- freq(z)

如果您有平面坐标参考系统（即单位为米或类似单位），您可以这样做

f <- cbind(f, area=f[,2] * prod(res(z)))

获取区域。但是对于 lon/lat 数据，您需要针对不同大小的单元格进行校正并执行

a <- area(z)
zonal(a, z, fun=sum)

并且有不同的绘图方式，但最基本的一种

plot(z)

绘制一个区域的斜率并返回高于和低于 R 中阈值的百分比

Mapping slope of an area and returning percent above and below a threshold in R

r

raster

terrain

r-raster

gis