验证 R 中多边形矩阵中的矢量坐标

Question

来自 R 中给定的矩阵和向量。验证一个点是否在一个区域中的短代码方法是什么？（目标：得到 True 或 False）

我目前尝试的是：

library(sf)
library(spatialEco)
# Lat long provided

Area <-  SpatialPolygons(list(Polygons(list(Polygon(do.call(rbind,list(
c(45.406164, 28.429255),
  c(44.182204, 27.395851),
  c(43.699651, 29.055894),
  c(45.259422,30.474260))))),'1')))

point=data.frame(t(c(44.590467, 28.057815)))
names(point)=c("lat", "long")
coordinates(point)=~lat+long

point.in.poly(point, Area)

我的期望是得到一个特征(True/False)来确认它是否在里面。

通过对我知道它不在多边形内的点执行相同操作，我看到 point.in.poly 输出中唯一发生变化的是 poly.ids 给出 NA.这是完成我想要检查的正确输出吗？

point2=data.frame(t(c(-44.590467, -28.057815)))
point.in.poly(point2, Area)

当然，还有更短的代码方式吗？

Answer 1

在我看来，有两种方法，具体取决于您需要多稳健的解决方案。如果您不需要高精度解决方案，您可以使用光线追踪 (ray-line intersection and point in polygon) 自己轻松编写函数。

我的意思是 "not accurate solution" 是您同意这样一个事实，即恰好位于线上或多边形节点等中的点将非常、非常、非常少地被错误分类。在 x64 计算机中，这个问题可以忽略不计，除非你确实构建了高精度三角测量算法或类似的东西（你可能不会在 R 中这样做）。

在这种情况下，这确实很简单，但如前所述，由于浮点算术精度缺陷，此类函数有时会失败。有关此问题的更多信息，请参阅计算几何书籍和 docs of CGAL。

基本函数可能如下所示：

section.ray.intersection <- function(p,d,a,b)
# this function tests if a ray r=[p,d] intersects the segment [a,b]
# p - point c(x,y)
# d - ray direction c(x,y)
# a - one of the nodes of the section [a,b] (not ordered with respect to computational geometry standards)
# b - remaining node of the section [a,b]
{
   v1 = p-a
   v2 = b-a
   v3 = c(-d[2],d[1])

   t1 = (v2[1]*v1[2] - v2[2]*v1[1])/(v2 %*% v3)
   t2 = (v1 %*% v3)/(v2 %*% v3) 

   return(t1 >=0. & (t2 >= 0. & t2 <= 1.0)) 
}

in.polygon <- function(px,py,Ax,Ay)
# basic function testing if the point is in the polygon
# px - x coordinate of the point
# py - y coordinate of the point
# Ax - a vector of x coordinates of the polygon
# Ay - a vector of y coordinates of the polygon    
{
    n <- length(Ax) # gent the number of nodes of the polygon
    if (n != length(Ay)) return(NA) # can be ommitted if Ax and Ay are always valid

    d = runif(2) # get the random vector of the direction
    is = 0 # number of intersections
    for (i in 2:n) # go trough all edges
        is <- is + section.ray.intersection(c(px,py),d,c(Ax[i-1],Ay[i-1]),c(Ax[i],Ay[i])) # and count TRUEs

    return( !(is %% 2==0)) # if the number of intersections is odd the point is inside
}

以及示例用法：

# libs
library(ggplot2)
library(purrr)
library(dplyr)

# polygon
Ax <- c(0,1,2,0.25,0)
Ay <- c(0,0,1,1.00,0)
# points
N = 1000 # number of points 
p <- data.frame(x = rnorm(N,0.5,.4), y= rnorm(N,0.5,.4)) # data frame od point coordinates

# apply the function to points - this is a "for loop" equivalent
p$is.inside <- apply(p, 1, function(X){return(in.polygon(X[1],X[2],Ax,Ay))})

# make a plot
ggplot() + geom_path(data=data.frame(x=Ax,y=Ay),aes(x,y)) + geom_point(data=p,aes(x,y,col = is.inside))

但是，如果您需要更准确的解决方案，您可以进行多次测试，然后使用一些投票方法。例如，假设在单个测试中分类不当的风险约为 p=0.001（实际上要低得多），那么在所有三个测试（使用随机方向向量）中分类不当的概率仅为 ppp = 0.000000001。该函数的修改版本可能如下所示：

in.polygon.robust <- function(px,py,Ax,Ay)
{
    n <- length(Ax)
    if (n != length(Ay)) return(NA) # can be ommitted if Ax and Ay are always valid

    d1 <- runif(2) # direction 1
    d2 <- runif(2) # direction 2
    d3 <- runif(2) # direction 3
    is1 <- 0 # number of intersections in direction 1
    is2 <- 0 # number of intersections in direction 2
    is3 <- 0 # number of intersections in direction 3
    for (i in 2:n) # go trough all edges
    {
        P <- c(px,py)
        Ap1 <- c(Ax[i-1],Ay[i-1])
        Ap2 <- c(Ax[i],Ay[i])
        # count intersections
        is1 <- is1 + section.ray.intersection(P,d1,Ap1,Ap2)
        is2 <- is2 + section.ray.intersection(P,d2,Ap1,Ap2)   
        is3 <- is3 + section.ray.intersection(P,d3,Ap1,Ap2)
    }
    r <- (is1 %% 2==0) + (is2 %% 2==0) + (is3 %% 2==0) # sum the even outcomes
    return( !(r>=2)) # return voting outcome
}

但是，如果仍然不够，您可以选择使用 CGAL via RCpp or (R bindings) 或者如果许可证不适合您，那么您必须编写自己的强大算法，这真的很难要做的事。

另外就是效率问题。如果您希望将该测试应用于大量多边形（特别是用大量边描述的），那么您必须使用一些特定的数据结构才能使一切都达到可接受的效率。可以找到更多详细信息 in this book or in this book。并全部用 C 或 C++ 编写。

验证 R 中多边形矩阵中的矢量坐标

Verify vector coordinates in polygon matrix in R

r

polygon