在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)来确认它是否在里面。

通过做同样的点，我知道它是不是在里面的poligon，我看到，唯一的变化在进入输出的 point.in.poly 是 poly.ids 给予 NA. 这样的输出是否能达到我想检查的目的？

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

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

Answer 1

在我看来，有两种方法，取决于你需要多么强大的灵魂。如果你不需要一个高精度的解决方案，你可以很容易地自己用射线跟踪写函数 (射线交点 和 点).

我所说的 "不精确的解决方案 "是指你同意这样一个事实，即精确地位于线上或多边形节点等处的点会被非常、非常、非常少地不当分类。在x64计算机中，这个问题是可以忽略的，除非你确实建立了一个高精度的三角测量算法或类似的东西（你可能不会在R中这样做）。

在这种情况下，要做的事情真的很简单，但正如前面提到的，这种函数有时会因为浮点算术精度缺陷而失败。关于这个问题的更多信息可以在计算几何学的书籍中找到，也可以在 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 通过RCpp或(R 绑定)或者如果授权不适合你，那么你就必须自己写一个强大的算法，这真的是一件很难的事情。

另外就是效率问题。如果你想把这个测试应用于大量的多边形（尤其是描述了大量的边），你就必须使用一些特定的数据结构，以使每件事都能被接受的效率。更多关于这方面的细节可以找到 本书中 或 本书中. 而且都是用C或C++写的。