Я работаю над моделированием изменения землепользования / покрытия в R. У меня есть две растровые карты, одна с «классами землепользования / покрытия», а другая с «информацией о риске обезлесения». Я использую растр риска, чтобы определить пиксели леса (один из классов землепользования / покрытия), которые с большей вероятностью будут обезлесены. Пока у меня есть R-код, который работает, вот пример, который можно воспроизвести:
#create land-use/cover raster
real <- raster(nrows=25, ncols=25, vals=round(rnorm(625, 3), 0))
real[ real > 3 ] <- 3 #just so we end with 3 land-use/cover classes
real[ real < 1 ] <- 1 #just so we end with 3 land-use/cover classes
plot(real)
#function to create the deforestation risk raster created by someone else
rtnorm <- function(n, mean = 0, sd = 1, min = 0, max = 1) {
bounds <- pnorm(c(min, max), mean, sd)
u <- runif(n, bounds[1], bounds[2])
qnorm(u, mean, sd)
}
risk <- raster(nrows=25, ncols=25, vals=rtnorm(n = 625, .1, .9)) #deforestation risk raster
plot(risk)
#The actual analysis starts here:
forest <- real #read the land-use/cover raster
forest[ forest != 3 ] <- NA #all left is class 3 (let's say, forest) in the raster
plot(forest, col="green")
deforestation <- sum(forest, risk) #identify the forest pixels with highest risk
plot(deforestation)
deforestation[ deforestation <= 3.5 ] <- 0 #rule to deforest the forest pixels
deforestation[ deforestation > 0 ] <- 100 #mark the pixels that will be deforested
plot(deforestation)
simulation <- sum(real, deforestation)
simulation[ simulation > 100 ] <- 2 #I use this to mark the forest pixels to a different land-use/cover class
plot(simulation)
Я хотел бы изменить правило, которое я использую, чтобы выбрать пиксели леса, которые будут вырублены (т. Е. deforestation[ deforestation <= 3.5 ] <- 0
). Вместо того, чтобы выбирать пороговое значение, такое как 3.5
, Интересно, могу ли я установить определенное количество пикселей леса (например, 50) для вырубки леса, а затем выбрать 50 пикселей леса с наибольшим риском вырубки леса.
Я совершенно не знаю, как сделать что-то подобное в R, поэтому любые предложения будут высоко оценены. Спасибо.