Как найти точки в контурах в R?

Question

Я изучил этот вопрос и создал свою собственную четырехконтурную карту, основанную на нескольких тысячах пар точек долготы и широты, но я не получаю правильное количество точек внутри каждой из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 не делаетне работает

Спасибо за любую помощь и дайте мне знать, если я смогу что-нибудь уточнить!

Answer 1

I думаю, 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