Определить, принадлежит ли данный lat-lon многоугольнику

Question

Предположим, у меня есть файл данных с именем zone с 1994 струны 2D координаты, обозначающие координаты вершин многоугольников, как показано ниже (самое первое число на правой стороне каждой линии обозначает zone)

c1 <- "1", "1 21, 31 50, 45 65, 75 80"

c2 <- "2", "3 20, 5 15, 2 26, 70 -85, 40 50, 60 80"

.....

c1993 <- "1993", "3 2, 2 -5, 0 60, 7 -58, -12 23, 56 611, 85 152"

c1994 <- "1994", "30 200, 50 -15, 20 260, 700 -850, -1 2, 5 6, 8 15"

Теперь я хочу манипулировать этими строками таким образом, чтобы получить случайную пару lat-lon (скажем 12 а также 20), Я мог бы сравнить, попал ли он в первый полигон, второй полигон, 3-й полигон,.... или 1994-й полигон. Решение о грубой силе: сравните x-coordinate (= 12) всем 4x-координаты и y-coordinate(= 20) to all the4Y-coordinates inc1andc2, respectively. The conclusion would be whether there is a valid **sandwich** inequality for each given coordinateИксandy`.

Например, используя процесс решения, как указано выше, точка (12,20) будет в с1, но не в с2.

Мой вопрос: как мне достичь этой цели в R?

Моя попытка: благодаря помощи Стефана Лорана я смог сгенерировать все матрицы, каждая с определенными размерами, которые хранят lat-lon пары всех вершин каждого многоугольника со следующим кодом:

 zone <- read_delim("[directory path to zone.csv file]", delim = ",", col_names = TRUE)
for(i in 1:nrow(zone)){
  zone$geo[i] = substr(zone$geo[i],10,135)
}
zone <- zone[complete.cases(zone),]

 Numextract <- function(string){
    unlist(regmatches(string, gregexpr("[[:digit:]]+\\.*[[:digit:]]*", string)))
 }

for(i in 1:nrow(zone)){
        poly1 <- matrix(as.numeric(Numextract(zone$geo[i])),i, ncol=2, byrow=TRUE)
        poly2 <- cbind(poly1, c(i))
}

Однако, как вы можете видеть, мне нужно найти способ индексировать все матрицы, соответствующие каждой зоне, которые были сгенерированы во время for() петля. Причина в том, что после этого я могу использовать другой for() цикл, чтобы определить, к какой зоне принадлежит точка!! Но я не смог понять это, поэтому кто-нибудь может помочь мне с подробным кодом?

Фактический набор данных
Набор данных зон и полигонов

Набор данных пар Lat-Lon

Answer 1

Сначала определите полигоны как матрицы, каждая строка представляет вершину:

poly1 <- rbind(c(1,21), c(31,50), c(45,65), c(75,80))
poly2 <- rbind(c(3,20), c(5,15), c(2,26), c(70,-85))

Определите точку для тестирования:

point <- c(12,20)

Теперь используйте pip2d функция ptinpoly пакет:

> library(ptinpoly)
> pip2d(poly1, rbind(point))
[1] -1
> pip2d(poly2, rbind(point))
[1] 1

Это значит (см. ?pip2d) что точка находится за пределами poly1 и внутри poly2,

Обратите внимание rbind(point) в pip2d, Мы используем rbind потому что в более общем случае мы можем запустить тест для нескольких точек в одном и том же многоугольнике.

Если вам нужна помощь, чтобы преобразовать

c1 <- "1 21, 31 50, 45 65, 75 80"

в

poly1 <- rbind(c(1,21), c(31,50), c(45,65), c(75,80))

тогда, возможно, вам следует открыть еще один вопрос.

редактировать

Хорошо, не открывайте другой вопрос. Вы можете действовать следующим образом.

c1 <- "1 21, 31 50, 45 65, 75 80"

Numextract <- function(string){
  unlist(regmatches(string, gregexpr("[[:digit:]]+\\.*[[:digit:]]*", string)))
}

poly1 <- matrix(as.numeric(Numextract(c1)), ncol=2, byrow=TRUE)

Который дает:

> poly1
     [,1] [,2]
[1,]    1   21
[2,]   31   50
[3,]   45   65
[4,]   75   80

2-й править

Для вашей второй проблемы ваши данные слишком велики. Единственное решение, которое я вижу, это разделить данные на более мелкие части.

Но, во-первых, кажется, что pip2d Функция также вызывает сбой сеанса R. Так что используйте другую функцию: pnt.in.poly из пакета SDMTools,

Вот небольшая модификация этой функции, которая ускоряет удаление ненужных выходов:

library(SDMTools)
pnt.in.poly2 <- function(pnts, poly.pnts){
  if (poly.pnts[1, 1] == poly.pnts[nrow(poly.pnts), 1] && 
      poly.pnts[1, 2] == poly.pnts[nrow(poly.pnts), 2]){ 
    poly.pnts = poly.pnts[-1, ]
  }
  out = .Call("pip", pnts[, 1], pnts[, 2], nrow(pnts), poly.pnts[,1], poly.pnts[, 2], nrow(poly.pnts), PACKAGE = "SDMTools")
  return(out)
}

Теперь, как уже было сказано, раскол lat_lon мелкими кусочками, длиной 1 миллион каждый (кроме последнего, меньшего размера):

lat_lon_list <- vector("list", 70)
for(i in 1:69){
  lat_lon_list[[i]] = lat_lon[(1+(i-1)*1e6):(i*1e6),]
}
lat_lon_list[[70]] <- lat_lon[69000001:nrow(lat_lon),]

Теперь запустите этот код:

library(data.table)
for(i in 1:70){
  DT <- data.table(V1 = pnt.in.poly2(lat_lon_list[[i]], polys[[1]]))
  for(j in 2:length(polys)){
    DT[, (sprintf("V%d",j)):=pnt.in.poly2(lat_lon_list[[i]], polys[[j]])]
  }
  fwrite(DT, sprintf("results%02d.csv", i))
  rm(DT)
}

Если это работает, он должен сгенерировать 70 CSV-файлов, result01.csv,..., result70.csvкаждый размер 1000000x1944 (кроме последнего, поменьше), тогда их можно открыть в Excel.

3-е редактирование

Я попробовал код, и у меня есть ошибка: Error: cannot allocate vector of size 7.6 Mb,

Нам нужно более тонкое расщепление:

lat_lon_list <- vector("list", 2*69+1)
for(i in 1:(2*69)){
  lat_lon_list[[i]] = lat_lon[(1+(i-1)*1e6/2):(i*1e6/2),]
}
lat_lon_list[[2*69+1]] <- lat_lon[69000001:nrow(lat_lon),]

for(i in 1:(2*69+1)){
  DT <- data.table(V1 = pnt.in.poly2(lat_lon_list[[i]], polys[[1]]))
  for(j in 2:length(polys)){
    DT[, (sprintf("V%d",j)):=pnt.in.poly2(lat_lon_list[[i]], polys[[j]])]
  }
  fwrite(DT, sprintf("results%02d.csv", i))
  rm(DT)
}