Функция для определения того, пересекаются ли ряд точек и ряд многоугольников

Question

У меня есть список точек и список многоугольников, и я хочу найти многоугольник, в котором находится каждая точка. Хотя использовать команду st_intersects в пакете sf довольно просто, я считаю, что это не хорошоздесь, поскольку точек и многоугольников так много, что результирующая матрица слишком велика для моего компьютера.

Я пытаюсь обойти это, создав функцию, которая может просматривать все многоугольники, чтобы найтиправильный для каждой точки.Я сделал один, но он не работал хорошо, и я понятия не имел, в чем проблема с ним.Так что мне нужна помощь!Спасибо!

Ниже описан порядок создания списка точек

pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 
2.5))) %>% st_sf() %>% mutate(namepoint=c(1,2,3))

Simple feature collection with 3 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 2.5 ymax: 2.5
epsg (SRID):    NA
proj4string:    NA
   namepoint        geometry
1         1 POINT (0.5 0.5)
2         2 POINT (1.5 1.5)
3         3 POINT (2.5 2.5)

Ниже описано, как создать список многоугольников

pol_down = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0))))
pol_up=pol_down+c(2,1) 
pol_add <- list(pol_down,pol_up) %>% st_sfc() %>% st_sf %>% 
mutate(namepgn=c('A','B'))

Simple feature collection with 2 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 3
epsg (SRID):    NA
proj4string:    NA
   namepgn                       geometry
1       A POLYGON ((0 0, 2 0, 2 2, 0 ...
2       B POLYGON ((2 1, 4 1, 4 3, 2 ...

Функция, которую я сделалis

 fun_pgn <- function(x,y){
                location_m<- st_intersects(x,
                                           y,
                                           sparse = F)
                name <- pull(y,1) %>% .[location_m]
                return(name)
             }

Используйте команду st_intersects, чтобы получить их отношения

st_intersects(pts,pol,sparse = F)

      [,1]  [,2]
[1,]  TRUE FALSE
[2,]  TRUE FALSE
[3,] FALSE  TRUE

Затем я хочу добавить переменную, например pgnname в данные ptsчтобы записать многоугольник для каждой точки,

        namepoint    pgnname    geometry
1               1       "A"     POINT (0.5 0.5)
2               2       "A"     POINT (1.5 1.5)
3               3       "B"     POINT (2.5 2.5)

Однако моя функция не работает.После запуска следующих кодов

pts %>% mutate(pgnname=fun_pgn(geometry,pol))

следует

Simple feature collection with 3 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 2.5 ymax: 2.5
epsg (SRID):    NA
proj4string:    NA
  namepoint pgnname        geometry
1         1       A POINT (0.5 0.5)
2         2       B POINT (1.5 1.5)
3         3    <NA> POINT (2.5 2.5)   

Имена полигонов неверны!Но ключевая часть моей функции работает хорошо.Например

> pull(pol,1) %>% .[st_intersects(pts[3,],pol,sparse = F)]
[1] "B"

Так что я в замешательстве.Что происходит, когда я запускаю свои функции?

Answer 1

Вы говорите, что не можете использовать st_intersects, потому что "результирующая матрица слишком велика для моего компьютера".Самое простое решение для этого - не генерировать плотную матрицу, которую вы используете с помощью sparse = F.Поведение по умолчанию st_intersects без этого аргумента заключается в создании разреженного списка, который вы пытаетесь сделать со своей функцией.

x = st_intersects(pts, pol)
pts %>% mutate(pgnname = pol$namepgn[unlist(x)])
# Simple feature collection with 3 features and 2 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 0.5 ymin: 0.5 xmax: 2.5 ymax: 2.5
# epsg (SRID):    NA
# proj4string:    NA
#   namepoint pgnname        geometry
# 1         1       A POINT (0.5 0.5)
# 2         2       A POINT (1.5 1.5)
# 3         3       B POINT (2.5 2.5)

Обратите внимание, что в случае вероятности, что точка можетпопадают в несколько полигонов, и вы хотите назначить только один полигон для этой точки, вместо этого вы можете использовать следующее:

pgn = sapply(x, `[`, 1)
pts %>% mutate(pgnname = pol$namepgn[pgn])