Test Whether A Point Is In A Polygon Or Not

[Fikile Kentane]

Hi Rick

Thanks so much for the "PtInPoly" function. Awesome.

I'm doing field data collection in a meter audit project. My field workers complete a Meter Audit FORM and attach lat/long GIS coordinates for the location of the meters. I cover many areas and most of the time I have to name the areas using my own names. This function helps me locate if the lat/long are inside designated areas.

How I use the "PtInPoly" function is the following: I call it is inside my function in an excel spreadsheet and pass it "X-coord:lat", "Y-coord:long"d a "2-dimensional array of polygon points". My polygon is a closed polygon.

There are about 5000 waypoints points my field database. So, for each waypoint I call the "PtInPoly" and test on which area the waypoint is located. I do this by running a loop for all areas and come out of the loop after finding the matching area.

Wanted to ask if I'm allowed to modify you function and pass it four variables instead of three. Fourth one will be used for looping through all the areas?

When I'm done. I will post the modified function here for you and everyone to see.

Best
Fikile

[Rick Rothstein]

Hi Rick

Thanks so much for the "PtInPoly" function. Awesome.

I'm doing field data collection in a meter audit project. My field workers complete a Meter Audit FORM and attach lat/long GIS coordinates for the location of the meters. I cover many areas and most of the time I have to name the areas using my own names. This function helps me locate if the lat/long are inside designated areas.

How I use the "PtInPoly" function is the following: I call it is inside my function in an excel spreadsheet and pass it "X-coord:lat", "Y-coord:long"d a "2-dimensional array of polygon points". My polygon is a closed polygon.

There are about 5000 waypoints points my field database. So, for each waypoint I call the "PtInPoly" and test on which area the waypoint is located. I do this by running a loop for all areas and come out of the loop after finding the matching area.

Wanted to ask if I'm allowed to modify you function and pass it four variables instead of three. Fourth one will be used for looping through all the areas?

When I'm done. I will post the modified function here for you and everyone to see.

Best
Fikile

Note... I removed your quote of my entire original post from your post... no need to burden this forum's servers with duplicate data that does not serve to any real reference purpose in your message.

I am glad that you found my code useful to you. As to your question... sure, feel free to modify it in anyway that is useful to you. None of the code I have posted in this sub-forum is copyrighted, so you can use or modify any of it as needed.

[Fikile Kentane]

Hi Rick

I have to thank you again for the "PtInPoly" function. I need some help please. First let me explain how I use the "PtInPoly" function.

I have a number of areas where my filed workers collect electricity meter audit data using Tablets. Each data set collected from the filed where a meter is located has a lat/long (GPS coordinates). There are a number of different areas with meters are located. I use the "PtInPoly" function to get the name of the area after checking if lat/long is in the area polygon. So for each lat/long I run a loop through all areas (18 times. There are 18 areas) and come out with the area name where the lat/long is inside an area.

Below is an example of one of the areas involved. The area name is called "eBodi" and I have 4 points defining he polygon of the Ebodi erea.

---------------------------------------------------------------------------------------
'ebodi
areaName(totNumAreasElm - 6) = "ebodi"
numOfPolygonPoints(totNumAreasElm - 6) = 4

'Define the (x:y) poligon points. The shape must be closed, so there must be one array element more than shape corners and
'the first and last coordinates must be the same

'Initialize polygon point1
pPoint(totNumAreasElm - 6, 0) = -31.2064891207045:
pPoint(totNumAreasElm - 6, 1) = 28.2321971696144

'Initialize polygon point2
pPoint(totNumAreasElm - 6, 2) = -31.2045994633557:
pPoint(totNumAreasElm - 6, 3) = 28.2303846938358

'Initialize polygon point3
pPoint(totNumAreasElm - 6, 4) = -31.2080756213509:
pPoint(totNumAreasElm - 6, 5) = 28.2279014276724

'Initialize polygon point4
pPoint(totNumAreasElm - 6, 6) = -31.209143785413:
pPoint(totNumAreasElm - 6, 7) = 28.2299482763019
------------------------------------------------------------------------------------------

Before I pass these points to "PtInPoly" I first pass them into a 2 dimensional array (pP) as I understand is the requirement for the third argument of the "PtInPoly" function. I have managed t do this successfully and checked this when i was debugging my code.

The problem is, the PtInPoly function does not work when I pass it pP in my VB code. But it works when I use and array in the worksheet. I've tested it using the same polygon points.



Below is my code that I use to first assign my polygon points (pPoint) to Rick's two dimensional polygon (pP)
I pass pP to the PtInPoly function together with the lat/long GPS coordinates
-------------------------------------------------------------------------------------------------------------------------
Dim areaCounter As Integer 'areaCounter loops through all area from 0 to 17
Dim areaLocated As String '

For areaCounter = 0 To totNumAreasElm Step 1 ' Loop through all the areas

MsgBox ("testing area : " + areaName(areaCounter))

'Array of polygon points (pP). Each point is lat/long (x/y). This is the polygon that passed to Rick function
ReDim pP(0 To (numOfPolygonPoints(areaCounter) - 1), 0 To 1) As Double

'Initialize the polygon Rick polygon
Dim ppCounter As Integer
Dim pointNumber As Integer
pointNumberIncrement = 0

For ppCounter = 0 To (numOfPolygonPoints(areaCounter) - 2) Step 1 'use the for loop to assign the my polygon points (pPoint) to Rick's two dimensional polygon (pP)
pP(ppCounter, 0) = pPoint(areaCounter, 0 + pointNumberIncrement)
pP(ppCounter, 1) = pPoint(areaCounter, 1 + pointNumberIncrement)
If (ppCounter = 0) Then
pP(ppCounter + (numOfPolygonPoints(areaCounter) - 1), 0) = pPoint(areaCounter, 0 + pointNumberIncrement)
pP(ppCounter + (numOfPolygonPoints(areaCounter) - 1), 1) = pPoint(areaCounter, 1 + pointNumberIncrement)
End If
pointNumberIncrement = 2
Next ppCounter

If (PtInPoly(x, y, pP)) Then 'if the lat/long is in the area under test, PtInPoly returns a TRUE and I use this to get out of the loop

areaLocated = areaName(areaCounter)
loopCounter = totNumAreasElm

End If

areaCounter = areaCounter + 1

Next areaCounter
-------------------------------------------------------------------------------------------------------------------
Please help me understand the problem that makes "PtInPoly" not to function

[Fikile Kentane]

Hi Rick

I have solved problem above.

Had to modify my code to use the polygon array exactly the way you specified.

Its now working

Thanks a million

[Fikile Kentane]

Hi Rick

I have another question.

I'm busy with a VBA code that determines if a lat/long point lies in a specified area and I have succesfuly used your "PtInPoly" function. It works well tnx again.

My areas are defined in kml files from Google earth. How do I uploaded an example of one one of my kml area files "Bhekela - v1.kml"?

My question is, in VBA, how do I automate the process of extracting the coordinates that form polygon points of the area. What I want to do is to read the "Bhekela - v1.kml" file and extract the coordinates inside the "<coordinates>" tags (see below)
"<coordinates>
28.25959909739643,-31.20619763209652,0 28.26071545752565,-31.19702465776928,0 28.25085444969957,-31.19927874343989,0 28.24923203195886,-31.20652374573175,0 28.25959909739643,-31.20619763209652,0
</coordinates>“.

I want to read the data inside the tags "<coordinates>" and present it is the format below. pPoint is my static array for polygon points representing an area. The data below represents polygon points of a Bhekela area.

pPoint(17, 1) = -31.2061976320965: pPoint(17, 2) = 28.2595990973964
pPoint(17, 3) = -31.1970246577692: pPoint(17, 4) = 28.2607154575256
pPoint(17, 5) = -31.1992787434398: pPoint(17, 6) = 28.2508544496995
pPoint(17, 7) = -31.2065237457317: pPoint(17, 8) = 28.2492320319588
pPoint(17, 9) = -31.2061976320965: pPoint(17, 10) = 28.2595990973964

Can anyone help show me how to do this? I know how to program in VB but not an expert.

Tnx in anticipation

[gonurvia]

Hi Rick, I am trying to use your function (PtInPoly) inside a function i've made (APPROACH). It's inside an If;

Function APPROACH(P1 As Range, x1 As Double, y1 As Double) As Double

Dim A1(3, 3) As Double
Dim A2(5, 2) As Double

For i = 1 To 5 Step 2
A1(1, 1) = P1(i, 1)
A1(1, 2) = P1(i, 2)
A1(1, 3) = P1(i, 3)
A1(2, 1) = P1(i + 1, 1)
A1(2, 2) = P1(i + 1, 2)
A1(2, 3) = P1(i + 1, 3)
A1(3, 1) = P1(i + 2, 1)
A1(3, 2) = P1(i + 2, 2)
A1(3, 3) = P1(i + 2, 3)

A2(1, 1) = P1(i, 1)
A2(1, 2) = P1(i, 2)
A2(2, 1) = P1(i + 1, 1)
A2(2, 2) = P1(i + 1, 2)
A2(3, 1) = P1(i + 3, 1)
A2(3, 2) = P1(i + 3, 2)
A2(4, 1) = P1(i + 2, 1)
A2(4, 2) = P1(i + 2, 2)
A2(5, 1) = P1(i, 1)
A2(5, 2) = P1(i, 2)

If PtInPoly(x1, y1, A2) = 0 Then
    APPROACH = AltXY(A1, x1, y1)
    Exit For
End If
Next i

End Function

I want the function, in case the point (x1,y1) is INSIDE the polygon, to apply function AltXY (If PtInPoly(x1, y1, A2) = 0 Then...) and if not, check for other areas but for some reason its giving me the wrong result. AltXY is not badly written, as if I wanted it to give me APPROACH = 10, it still keeps giving me zeros but I don't understand why. Any help is appreciated

[Rick Rothstein]

Hi Rick, I am trying to use your function (PtInPoly) inside a function i've made (APPROACH). It's inside an If;

Function APPROACH(P1 As Range, x1 As Double, y1 As Double) As Double

Dim A1(3, 3) As Double
Dim A2(5, 2) As Double

For i = 1 To 5 Step 2
A1(1, 1) = P1(i, 1)
A1(1, 2) = P1(i, 2)
A1(1, 3) = P1(i, 3)
A1(2, 1) = P1(i + 1, 1)
A1(2, 2) = P1(i + 1, 2)
A1(2, 3) = P1(i + 1, 3)
A1(3, 1) = P1(i + 2, 1)
A1(3, 2) = P1(i + 2, 2)
A1(3, 3) = P1(i + 2, 3)

A2(1, 1) = P1(i, 1)
A2(1, 2) = P1(i, 2)
A2(2, 1) = P1(i + 1, 1)
A2(2, 2) = P1(i + 1, 2)
A2(3, 1) = P1(i + 3, 1)
A2(3, 2) = P1(i + 3, 2)
A2(4, 1) = P1(i + 2, 1)
A2(4, 2) = P1(i + 2, 2)
A2(5, 1) = P1(i, 1)
A2(5, 2) = P1(i, 2)

If PtInPoly(x1, y1, A2) = 0 Then
    APPROACH = AltXY(A1, x1, y1)
    Exit For
End If
Next i

End Function

I want the function, in case the point (x1,y1) is INSIDE the polygon, to apply function AltXY (If PtInPoly(x1, y1, A2) = 0 Then...) and if not, check for other areas but for some reason its giving me the wrong result. AltXY is not badly written, as if I wanted it to give me APPROACH = 10, it still keeps giving me zeros but I don't understand why. Any help is appreciated

You did not supply an x1,y1 point that fails to work so that I could test it out; however, I am suspicious of your A1 array declaration (probably your A2 array's declaration as well). You have the A1 array declared as...

Dim A1(3, 3) As Double

but then you start filling it at element number 1. Unless you are using Option Base 1, your declaration set the lower bound for the array at 0, not 1 and that might be screwing up the calculations in my code (not sure, but I think it would). Try changing your declaration to this and see if that solves the problem...

Dim A1(1 To 3, 1 To 3) As Double

This declaration set the lower bound to 1, not 0. While I don't know if it matters to your AltXY function or not, but for consistency, I would change the declaration for the A2 array in a similar way.

[Phil_n_IN]

Hi Rick,
I know this is a very old post but it's exactly what I'm looking and works great in excel! Thanks!
But now I need to figure out how to modify it so it can works in MS Access (I'm a newbie to Access).

I have two tables in Access 2016, one of points (x/y coordinates) and another of polygons (consisting of an Poly_ID and a series of x/y nodes defining each poly perimeter). I need to discover what points are in what polygon. The actual database has over 250K points and 4K polygons.

I'd like to use this function in Access but can't get my head around how to refer to a range in access. I've tried numerous methods but all are clearly not heading toward a viable solution.
Can you point me in the right direction?

Attached is a pic and txt delimited data with a few poly and points tables that i'm using for the test. (The system won't let me attach a accdb file)
2017-07-14 12_47_21-Access - TEST Point-in-Polygon _ Database- E__MLS_RETS_DATA_TEST Point-in-Po.jpg

[Phil_n_IN]

Hi Rick,
I've figure out how to make it work in Access for a single polygon with multiple points, but this is only half the solution. I'm now stuck on how to get it to work with multiple polygons.
Using the previously attached files I created a query that limits the results to one polygon. Using your code as a base, it lists which points are within the polygon and identifies the polygons ID.

I'd appreciate any pointers you may have on how I could proceed with the next step of making it work with multiple polygons?

Thanks so much in advance

Option Explicit

Public Function PtInPoly(Xcoord As Double, Ycoord As Double) As Variant
Dim X As Long, inPoly As String, NumSidesCrossed As Long, m As Double, b As Double, Poly As Variant
Dim Xx As Long, Yy As Long, Xupper As Long, Yupper As Long, transposeArray As Variant
Dim dbs As DAO.Database
Dim Polyrst As DAO.Recordset

Set dbs = CurrentDb
Set Polyrst = dbs.OpenRecordset("SELECT x_nodes, y_nodes ,Poly_ID FROM Poly_ID_2only", dbOpenSnapshot)
    With Polyrst
        .MoveLast
        .MoveFirst
        Poly = .GetRows(.RecordCount)
    End With
'GetRows() is weird in that it returns rows & columns horizontally,
' the code below "transposes" the data to read down instead of across
    Xupper = UBound(Poly, 2)
    Yupper = UBound(Poly, 1)
    ReDim transposeArray(Xupper, Yupper)
    For Xx = 0 To Xupper
        For Yy = 0 To Yupper
             transposeArray(Xx, Yy) = Poly(Yy, Xx)
        Next Yy
    Next Xx
    Poly = transposeArray
'-----------------------------------------------------------
Debug.Print UBound(Poly) + 1 & " records retrieved."
    For X = LBound(Poly) To UBound(Poly) - 1
        If Poly(X, 0) > Xcoord Xor Poly(X + 1, 0) > Xcoord Then
            m = (Poly(X + 1, 1) - Poly(X, 1)) / (Poly(X + 1, 0) - Poly(X, 0))
            b = (Poly(X, 1) * Poly(X + 1, 0) - Poly(X, 0) * Poly(X + 1, 1)) / (Poly(X + 1, 0) - Poly(X, 0))
            If m * Xcoord + b > Ycoord Then NumSidesCrossed = NumSidesCrossed + 1
        End If
    Next
Debug.Print NumSidesCrossed + 1; "Lines Crossed"

If CBool(NumSidesCrossed Mod 2) = True Then
    inPoly = Poly(0, 2)
Else
    inPoly = "not in polygon"
End If

PtInPoly = inPoly
End Function

TEST-Results Point-in-Polygo.jpg

results shown in pic

[vstepaniuk]

This is a working code, which doesn't require the first and last points to be the same. The function is called IsInside() here.

Public Function IsInside(Xcoord As Double, Ycoord As Double, Polygon As Variant) As Variant
  Dim x As Long, NumSidesCrossed As Long, m As Double, b As Double, Poly As Variant
  Poly = Polygon
  If (UBound(Poly, 2) - LBound(Poly, 2)) <> 1 Then
    If TypeOf Application.Caller Is Range Then
      IsInside = "#WrongNumberOfCoordinates!"
    Else
      Err.Raise 999, , "Array Has Wrong Number Of Coordinates!"
    End If
    Exit Function
  End If
  For x = LBound(Poly) To UBound(Poly)
    If x < UBound(Poly) Then
      If Poly(x, 1) > Xcoord Xor Poly(x + 1, 1) > Xcoord Then
        m = (Poly(x + 1, 2) - Poly(x, 2)) / (Poly(x + 1, 1) - Poly(x, 1))
        b = (Poly(x, 2) * Poly(x + 1, 1) - Poly(x, 1) * Poly(x + 1, 2)) / (Poly(x + 1, 1) - Poly(x, 1))
        If m * Xcoord + b > Ycoord Then NumSidesCrossed = NumSidesCrossed + 1
      End If
    Else
      If Poly(UBound(Poly), 1) > Xcoord Xor Poly(LBound(Poly), 1) > Xcoord Then
        m = (Poly(LBound(Poly), 2) - Poly(UBound(Poly), 2)) / (Poly(LBound(Poly), 1) - Poly(UBound(Poly), 1))
        b = (Poly(UBound(Poly), 2) * Poly(LBound(Poly), 1) - Poly(UBound(Poly), 1) * Poly(LBound(Poly), 2)) / (Poly(LBound(Poly), 1) - Poly(UBound(Poly), 1))
        If m * Xcoord + b > Ycoord Then NumSidesCrossed = NumSidesCrossed + 1
      End If
    End If
  Next
  IsInside = CBool(NumSidesCrossed Mod 2)
End Function