Отображение наклона области и возврат процентов выше и ниже порога в R

Я пытаюсь изобразить нашу долю площади с уклоном 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 градусов - один цвет, назовите его хорошим или зеленым, а все остальное - другим цветом, назовите его плохим или красным.

Спасибо

Источник

2 ответа

Решение

Ты можешь использовать reclassify от raster пакет для достижения этого. Функция присваивает каждому значению ячейки, которое находится в пределах определенного интервала, определенное значение. Например, вы можете назначить значения ячеек в пределах интервала (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
Источник

Там нет отрицательных склонов, поэтому я предполагаю, что вы хотите те, которые меньше 5 градусов

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

z <- x <= 5

Теперь вы можете считать клетки с частотой

f <- freq(z)

Если у вас есть плоская система координат координат (то есть с единицами измерения в метрах или аналогичными), вы можете сделать

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

чтобы получить районы. Но для данных lon / lat вам нужно будет исправить ячейки разных размеров и сделать

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

И есть разные способы заговора, но самый основной

plot(z)
Источник
Другие вопросы по тегам