使用 Raycasting 算法进行具有纬度/经度坐标的 Point-In-Polygon 测试

Posted 2023-02-23

【中文标题】使用 Raycasting 算法进行具有纬度/经度坐标的 Point-In-Polygon 测试

【英文标题】：Using a Raycasting Algorithm for Point-In-Polygon test with latitude/longitude coordinates

【发布时间】：2015-03-26 16:28:52

【问题描述】：

我已经成功使用 DotSpatial.Contains 函数来测试我的每个点（180 万）是否位于我的 shapefile 中。但是，该算法非常慢，因为我正在测试大量点和非常复杂的多边形。这是一个例子： http://picload.org/image/igararl/pointshapesel.png

图片中德国的边界是我的shapefile，用来选点（简化了，不过还是14000个顶点），红色矩形是我的180万个点所在的区域。

在为空间纬度/经度坐标寻找更快的多边形点测试时，我遇到了光线投射算法： http://alienryderflex.com/polygon/

我将代码翻译成 VB.Net，它运行没有错误，但它没有找到任何点/shapefile 的交集。我知道经纬度坐标的困难 - 但在德国地区，经纬度坐标与标准笛卡尔坐标系相匹配。

这是我的（该）代码。出于速度的原因，我首先声明了全局变量：

Public polyCorners As Integer
Public polyX() As Double
Public polyY() As Double
Public xP, yP As Double
Public constant() As Double
Public multiple() As Double

然后我将我的 Shapefile 顶点添加到 polyCorners 列表中（这有效）：

  Dim ShapefilePoly As Shapefile = Shapefile.OpenFile(TextBox4.Text)
    Dim x As Long = 1
    For Each MyShapeRange As ShapeRange In ShapefilePoly.ShapeIndices
        For Each MyPartRange As PartRange In MyShapeRange.Parts
            For Each MyVertex As Vertex In MyPartRange
                If MyVertex.X > 0 AndAlso MyVertex.Y > 0 Then
                    pointsShape.Add(New PointLatLng(MyVertex.Y, MyVertex.X))
                    ReDim Preserve polyY(x)
                    ReDim Preserve polyX(x)
                    polyY(x) = MyVertex.Y
                    polyX(x) = MyVertex.X
                    x = x + 1
                End If
            Next
        Next
    Next
    ReDim constant(x)
    ReDim multiple(x)

在实际搜索之前，我按照作者的建议调用 precalc_values()：

    Private Sub precalc_values()

    Dim i As Integer, j As Integer = polyCorners - 1

    For i = 0 To polyCorners - 1
        If polyY(j) = polyY(i) Then
            constant(i) = polyX(i)
            multiple(i) = 0
        Else
            constant(i) = polyX(i) - (polyY(i) * polyX(j)) / (polyY(j) - polyY(i)) + (polyY(i) * polyX(i)) / (polyY(j) - polyY(i))
            multiple(i) = (polyX(j) - polyX(i)) / (polyY(j) - polyY(i))
        End If
        j = i
    Next
End Sub

最后，我为每个 lat/lng 点调用 pointInPolygon()：

Function LiesWithin(latP As Double, lngP As Double) As Boolean
    LiesWithin = False
    xP = lngP
    yP = latP
    If pointInPolygon() = True Then LiesWithin = True
End Function

Private Function pointInPolygon() As Boolean

    Dim i As Integer, j As Integer = polyCorners - 1
    Dim oddNodes As Boolean = False

    For i = 0 To polyCorners - 1
        If (polyY(i) < yP AndAlso polyY(j) >= yP OrElse polyY(j) < yP AndAlso polyY(i) >= yP) Then
            oddNodes = oddNodes Xor (yP * multiple(i) + constant(i) < xP)
        End If
        j = i
    Next

    Return oddNodes
End Function

所有变量似乎都正确填充，数组包含我的多边形角，点列表从第一个点到最后一个点都被准确检查。它在 20 秒内运行了 180 万个点的完整列表（相比之下，使用 DotSpatial.Contains 函数需要 1 小时 30 分钟）。任何人都知道为什么它没有找到任何相交点？

【参考方案1】：

嗯，我比预期更快地发现了我的问题： 我忘记将 Shapefile 顶点的数量分配给 polyCorners。 在我上面的代码中，只需在

之后添加 polyCorners = x

   ReDim constant(x)
   ReDim multiple(x)
   polyCorners = x

也许有人觉得这段代码很有用。我真的很惊讶它的速度有多快！

【讨论】：

我一定会试一试的。