Using a Raycasting Algorithm for Point-In-Polygon test with latitude/longitude coordinates

Question

I already successfully used the DotSpatial.Contains function for testing whether each of my points (1.8 million) lie within my shapefile or not. However, the algorithm is very slow, as I am testing against a large number of points and very complex polygons. Here's an example: http://picload.org/image/igararl/pointshapesel.png

The boundary of Germany in the image is my shapefile, which is used for selecting points (simplified, but still 14.000 vertices), and the red rectangle is the area where my 1.8 million points are.

In search for a faster point-in-polygon test for spatial latitude/longitude coordinates, I came across the Raycasting Algorithm: http://alienryderflex.com/polygon/

I translated the code to VB.Net, and it runs without error, but it is not finding any intersection of points/shapefile. I am aware of the difficulties with lat/long coordinates - but in the area of Germany, lat/long coordinates match the standard Cartesian coordinate system.

Here is my (the) code. I declare the global variables for reasons of speed first:

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

Then I am adding my Shapefile vertices to the list of polyCorners (this works):

  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)

Before the actual search, I am calling precalc_values() as is suggested by the author:

    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

Finally, I am calling pointInPolygon() for each of my lat/lng points:

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

All variables seem to get filled correctly, the array contains my polygon corners and the list of points is checked accurately from the first point to the last. It is running through the full list of 1.8 Million points in under 20 seconds (compared to 1 hour and 30 minutes using the DotSpatial.Contains function). Anyone has an idea why it is not finding any intersecting points?

Answer 1

Well, I found my problem quicker than expected: I forgot to assign the number of Shapefile vertices to polyCorners. In my above code, just add polyCorners = x after

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

Maybe someone finds this code useful. I am really surprised how super-fast it is!