确定给定的经纬度是否属于多边形

Question

假设我有一个名为zone 的数据文件，其中1994 字符串2D 坐标表示多边形顶点的坐标，如下所示（每行RHS 上的第一个数字表示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），我可以比较它是否属于第一个多边形，第二个多边形，第三个多边形，....或第 1994 个多边形。 蛮力解决方案是：将x-coordinate (= 12) 与所有4 x-坐标和y-coordinate(= 20) to all the4y@987654342 @c1andc2, respectively. The conclusion would be whether there is a valid **sandwich** inequality for each given coordinatexandy`。

例如，通过使用上述求解过程，点(12,20) 将在c1 但不在c2。

我的问题：我怎样才能在 R 中实现这个目标？

我的尝试：感谢 Stéphane Laurent 的帮助，我能够生成所有矩阵，每个矩阵都有一定的大小，存储每个多边形的所有顶点的 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() 循环来确定一个点属于哪个区域！但是我一直无法弄清楚，所以任何人都可以帮我提供详细的代码吗？

实际数据集
Zone and polygons dataset

Lat-Lon pairs dataset

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)

现在，使用ptinpoly 包的pip2d 函数：

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

这意味着（参见?pip2d）该点在poly1 之外和poly2 之内。

注意pip2d 中的rbind(point)。我们使用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

第二次编辑

对于您的第二个问题，您的数据太大。我能看到的唯一解决方案是将数据拆分成更小的部分。

但首先，pip2d 函数似乎也会导致 R 会话崩溃。所以使用另一个函数：来自包SDMTools 的pnt.in.poly。

这里是这个函数的一个小修改，通过删除无用的输出使其更快：

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 分割成更小的部分，每个长度为 100 万，（除了最后一个，更小）：

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 中打开它们。

第三次编辑

我已经尝试了代码，但出现错误：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)
}