Как найти точки в контурах в R?
Я изучил этот вопрос и создал свою собственную четырехконтурную карту, основанную на нескольких тысячах пар точек долготы и широты, но я не получаю правильное количество точек внутри каждого из 4 контуров, используя points.in.polygon метод, упомянутый в приведенном выше вопросе.
Вот код, использующий библиотеку MASS:
# use kde2d function to create kernel density estimates
x <- pedestrian.df$longitude
y <- pedestrian.df$latitude
dens <- kde2d(x, y, n=200)
# create the contours to plot - 70%, 50%, 25%, 10% of density contained in each contour
prob <- c(0.7, 0.5, 0.25, 0.1)
dx <- diff(dens$x[1:4])
dy <- diff(dens$y[1:4])
sz <- sort(dens$z)
c1 <- cumsum(sz) * dx * dy
levels <- sapply(prob, function(x) {
approx(c1, sz, xout = 1 - x)$y
})
#create the contour plot using smoothScatter which smooths the collisions into kernel densities
smoothScatter(x,y) + contour(dens, levels=levels, labels=prob, col = c("green", "yellow", "orange", "red"), lwd = 1.5, add=T)
Это правильно генерирует то, что я ожидал:
Затем я попытался использовать функцию points.in.polygon из библиотеки sp, как в ответе на связанный выше вопрос:
ls <- contourLines(dens, level=levels)
zone_1 <- point.in.polygon(df$longitude, df$latitude, ls[[4]]$x, ls[[4]]$y)
zone_2 <- point.in.polygon(df$longitude, df$latitude, ls[[3]]$x, ls[[3]]$y)
zone_3 <- point.in.polygon(df$longitude, df$latitude, ls[[2]]$x, ls[[2]]$y)
zone_4 <- point.in.polygon(df$longitude, df$latitude, ls[[1]]$x, ls[[1]]$y)
Но это приводит к неправильному количеству точек на зону или контур. Я знаю, что это неправильно, потому что у каждого контура должно быть постепенно больше точек, поскольку контур становится больше.
Я попытался посмотреть на ls (список, в котором хранится список всех координат x и y полигонов), но есть 15 уровней, а не 4, которые я интуитивно предположил бы, будет там. Среди 15 есть даже несколько уровней, имеющих одинаковое значение. Я подозреваю, что ответ на мою проблему заключается в правильном подборе этого списка списков, чтобы включить 4 уровня, которые соответствуют моим 4 контурам, но ls[[1:7]]$x, ls[[1:7]]$y не делает не работает
Спасибо за любую помощь и дайте мне знать, если я смогу уточнить что-нибудь!
1 ответ
Я думаю pedestrian Ваши собственные данные сравниваются с чем-то в pkg, и поскольку это не является частью вопроса, мы будем использовать другой:
library(MASS)
library(sp)
attach(geyser)
data.frame(
x = geyser$duration,
y = geyser$waiting
) -> xdf
dens <- kde2d(xdf$x, xdf$y, n = 100)
prob <- c(0.7, 0.5, 0.25, 0.1)
dx <- diff(dens$x[1:4])
dy <- diff(dens$y[1:4])
sz <- sort(dens$z)
c1 <- cumsum(sz) * dx * dy
levels <- sapply(prob, function(x) {
approx(c1, sz, xout = 1 - x)$y
})
smoothScatter(x,y) +
contour(dens, levels=levels, labels=prob, col = c("green", "yellow", "orange", "red"), lwd = 1.5, add=TRUE)
Причина "нескольких уровней" заключается в том, что каждый многоугольник в данном слое является отдельным, поэтому потенциально существует> 1 на уровень:
cl <- contourLines(dens, level=levels)
sort(table(sapply(cl, `[[`, "level")))
## 0.00519851181336958 0.00765971436995347 0.0107843979424224 0.0128423136194731
## 2 3 3 3
Итак, просто учтите это при расчете точек на полигон:
setNames(
lapply(cl, function(poly) sum(sp::point.in.polygon(xdf$x, xdf$y, poly$x, poly$y))),
sapply(cl, `[[`, "level")
) -> level_cts
str(level_cts)
## List of 11
## $ 0.00519851181336958: int 91
## $ 0.00519851181336958: int 174
## $ 0.00765971436995347: int 78
## $ 0.00765971436995347: int 57
## $ 0.00765971436995347: int 74
## $ 0.0107843979424224 : int 65
## $ 0.0107843979424224 : int 34
## $ 0.0107843979424224 : int 33
## $ 0.0128423136194731 : int 42
## $ 0.0128423136194731 : int 10
## $ 0.0128423136194731 : int 3
Тогда мы можем подвести их итог:
sapply(
split(level_cts, names(level_cts)),
function(level) sum(unlist(level))
) -> pt_cts
pt_cts
## 0.00519851181336958 0.00765971436995347
## 265 209
## 0.0107843979424224 0.0128423136194731
## 132 55
И получить %:
pt_cts / nrow(xdf)
## 0.00519851181336958 0.00765971436995347
## 0.8862876 0.6989967
## 0.0107843979424224 0.0128423136194731
## 0.4414716 0.1839465
ОБНОВИТЬ
Вместо того, чтобы просто вычислять проценты, мы также можем назначить уровень исходным данным:
do.call(
rbind.data.frame,
lapply(cl, function(poly) { # iterate over each polygon
# figure out which pts are in this polgyon
which_pts <- as.logical(sp::point.in.polygon(xdf$x, xdf$y, poly$x, poly$y))
tdf <- xdf[which_pts,] # assign them to a temp data frame
tdf$level <- poly$level # add the level
tdf
})
) -> new_xdf
dplyr::glimpse(new_xdf)
## Observations: 661
## Variables: 3
## $ x <dbl> 2.000000, 2.033333, 1.833333, 1.616667, 1.766667, 2.0000...
## $ y <dbl> 77, 77, 81, 89, 73, 83, 84, 85, 79, 75, 91, 87, 86, 78, ...
## $ level <dbl> 0.005198512, 0.005198512, 0.005198512, 0.005198512, 0.00...
# while you likely want the level value, this adds columns for level # & prob
new_xdf$level_num <- as.integer(factor(new_xdf$level, levels, labels=1:length(levels)))
new_xdf$prob <- as.numeric(as.character(factor(new_xdf$level, levels, labels=prob)))
dplyr::glimpse(new_xdf)
## Observations: 661
## Variables: 5
## $ x <dbl> 2.000000, 2.033333, 1.833333, 1.616667, 1.766667, 2....
## $ y <dbl> 77, 77, 81, 89, 73, 83, 84, 85, 79, 75, 91, 87, 86, ...
## $ level <dbl> 0.005198512, 0.005198512, 0.005198512, 0.005198512, ...
## $ level_num <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ prob <dbl> 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0....
dplyr::count(new_xdf, level, level_num, prob)
## # A tibble: 4 x 4
## level level_num prob n
## <dbl> <int> <dbl> <int>
## 1 0.00520 1 0.700 265
## 2 0.00766 2 0.500 209
## 3 0.0108 3 0.250 132
## 4 0.0128 4 0.100 55