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

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 the4</code>y<code>-coordinates inc1andc2, 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)
}

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

Determine if a given lat-lon belong to a polygon

r

point-in-polygon

编辑

第二次编辑

第三次编辑