Test Whether A Point Is In A Polygon Or Not

**qetuol** · 03-15-2016, 01:51 AM

Originally Posted by Rick Rothstein

The Polygon argument can be either a two-dimensional VBA array or a two-dimensional range. If it is a range, assigning that range to a Variant variable converts it to a two-dimensional VBA array so that the rest of the code ends up with the same kind of array as when a real two-dimensional array is inputted. The reason the assignment has to take place within the code is because the coercion to an array does not take place via the argument itself even if that argument is a Variant. So, that line of code makes the overall function flexible enough to handle a real array (assigning such an array to the Variant variable does not change its structure) or a range (because the range gets converted to a real array).

I wanted the function to return either TRUE or FALSE directly in case the function was to be used directly in a MessageBox or to fill a TextBox or any other number of uses besides as the argument to an IF function... the idea being the user can implement the function without having to take special steps to format the output for the occasion. As for why PtInPoly is a Variant and not a Boolean... that was so it could accommodate the error messages when called as a UDF directly on a worksheet.

It was a design decision to help a user not make a mistake. If the user accidentally grabbed the wrong number of rows from a table of coordinates and inputted that too short range of values, the function would alert the user to the error because the last point did not equal the first point... if I made the function close the polygon automatically, it would blindly calculate values and the user would not be aware some of the answers were incorrect... by forcing the user to make sure the polygon closed, he/she could never make a mistake of accidentally inputting to few points.

Rick, thank you for the quick response and clarifications. Anyway, if i want to include each vertex of the polygon only once, not doubling the first and last point, how should i modify the code? is it correct if

Code:

For x = LBound(Poly) To UBound(Poly) - 1

is changed to

Code:

For x = LBound(Poly) To UBound(Poly) - 2

? thanks

**Rick Rothstein** · 03-15-2016, 02:29 AM

Originally Posted by qetuol

Rick, thank you for the quick response and clarifications. Anyway, if i want to include each vertex of the polygon only once, not doubling the first and last point, how should i modify the code? is it correct if

Code:

For x = LBound(Poly) To UBound(Poly) - 1

is changed to

Code:

For x = LBound(Poly) To UBound(Poly) - 2

? thanks

**hugohonda** · 09-22-2016, 01:01 AM

Originally Posted by Rick Rothstein

This function is not one that many will find a use for, but if you ever need its functionality, then here it is. The function can be called from other VB code or used as a UDF (user defined function) directly on a worksheet. What it does is tell you whether a point is located inside a polygon (simple or complex, convex or concave) or not. That's it... if you should ever need such a function, this is the code for it...

Code:

Public Function PtInPoly(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 Not (Poly(LBound(Poly), 1) = Poly(UBound(Poly), 1) And _
        Poly(LBound(Poly), 2) = Poly(UBound(Poly), 2)) Then
    If TypeOf Application.Caller Is Range Then
      PtInPoly = "#UnclosedPolygon!"
    Else
      Err.Raise 998, , "Polygon Does Not Close!"
    End If
    Exit Function
  ElseIf UBound(Poly, 2) - LBound(Poly, 2) <> 1 Then
    If TypeOf Application.Caller Is Range Then
      PtInPoly = "#WrongNumberOfCoordinates!"
    Else
      Err.Raise 999, , "Array Has Wrong Number Of Coordinates!"
    End If
    Exit Function
  End If
  For x = LBound(Poly) To UBound(Poly) - 1
    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
  Next
  PtInPoly = CBool(NumSidesCrossed Mod 2)
End Function

The theory behind the function is simplicity itself... start at the point being tested and project a line from that point outward in any direction (I chose straight up as that made some of the math easier) and count how many polygon sides it crosses... if the number of sides is odd, the point lies inside the polygon and if the number of sides is even, the point lies outside of the polygon. The function will return either True, False (for inside/outside respectively) or an error message (see below). The first argument is the X-Coordinate and the second argument is the Y-Coordinate of the point you want to test for being inside the polygon or not. The third argument is either a two-dimensional array of numbers (when called from another VB code procedure) or a range of numbers consisting of two columns and as many rows as needed (when called from a worksheet). The numbers describe the nodes composing the polygon.

NOTE #1: The polygon must be closed, meaning the first listed point and the last listed point must be the same. If they are not the same, the function will raise "Error #998 - Polygon Does Not Close!" if the function was called from other VB code or it will return #UnclosedPolygon! if called from the worksheet. Normally, if called from a worksheet, you would probably be using the function in a formula something like this...

=IF(PtInPoly(B9,C9,E18:F37),"In Polygon","Out Polygon")

In that case, the formula will return a #VALUE! error, not the #UnclosedPolygon! error, because the returned value to the IF function is not a Boolean; however, if you select the "PtInPoly(B9,C9,E18:F37)" part of the function in the Formula Bar and press F9, it will show you the returned value from the PtInPoly function as being #UnclosedPolygon!.

NOTE #2: The range or array specified for the third argument must be two-dimensional. If it is not, then the function will raise "Error #999 - Array Has Wrong Number Of Coordinates!" if the function was called from other VB code or it will return #WrongNumberOfCoordinates! if called from the worksheet. Error reporting when called from the worksheet will be the same as described in NOTE #1.

NOTE #3: Points extremely close to, or theoretically lying on, the polygon borders may or may not report back correctly... the vagrancies of floating point math, coupled with the limitations that the "significant digits limit" in VBA imposes, can rear its ugly head in those circumstances producing values that can calculate to either side of the polygon border being tested (remember, a mathematical line has no thickness, so it does not take too much of a difference in the significant digits to "move" a calculated point's position to one side or the other of such a line).

NOTE #4: While I think error checking is a good thing, the setup for this function is rather straightforward and, with the possible exception of the requirement for the first and last point needing to be the same, easy enough for the programmer to maintain control over during coding. If you feel confident in your ability to meet the needs of NOTE #1 and NOTE #2 without having the code "looking over your shoulder", then the function can be simplified down to this compact code...

Code:

Public Function PtInPoly(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
  For x = LBound(Poly) To UBound(Poly) - 1
    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
  Next
  PtInPoly = CBool(NumSidesCrossed Mod 2)
End Function

Hi Rick,

first thanks for sharing your code, it has been very useful to me. Similar to theCloud, i work with aircraft center of gravity and needed a function similar to yours. However I had some problems:
1- when a point was on a vertical boundarie
2- when a point was on an horizontal boundarie
3- when comparing two Double variable (b = Ycoord)...so I replace Double by Single and it works for now...I might need to use Double in the future, we'll see how it goes

maybe this will help somebody, if you see any problems with my changes let me know

thanks

Code:

Function PtInPoly(Xcoord As Double, Ycoord As Double, Polygon As Variant) As Variant
Dim X As Long, NumSidesCrossed As Long, m As Single, b As Single, Poly As Variant '*********changed m and b variable type from Double to Single
Dim btest As Single
Poly = Polygon
If Not (Poly(LBound(Poly), 1) = Poly(UBound(Poly), 1) And Poly(LBound(Poly), 2) = Poly(UBound(Poly), 2)) Then
    If TypeOf Application.Caller Is Range Then
        PtInPoly = "#UnclosedPolygon!"
    Else
        Err.Raise 998, , "Polygon Does Not Close!"
    End If
    Exit Function
ElseIf UBound(Poly, 2) - LBound(Poly, 2) <> 1 Then
    If TypeOf Application.Caller Is Range Then
        PtInPoly = "#WrongNumberOfCoordinates!"
    Else
        Err.Raise 999, , "Array Has Wrong Number Of Coordinates!"
    End If
    Exit Function
End If
For X = LBound(Poly) To UBound(Poly) - 1
Continu:
    If Poly(X, 1) >= Xcoord Xor Poly(X + 1, 1) >= Xcoord Then                '************changed > for >= to accomodate a point located on a vertical boundarie limit
        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 = 0 And b = Ycoord Then                                         '********* added this "if statement" to check if point is on horizontal boundaries limit
            If X < UBound(Poly) - 1 Then
                X = X + 1
                GoTo Continu:
            Else
                Exit For
            End If
        Else
            If m * Xcoord + b > Ycoord Then NumSidesCrossed = NumSidesCrossed + 1
        End If
    End If
Next X
PtInPoly = CBool(NumSidesCrossed Mod 2)
End Function

**Rick Rothstein** · 09-22-2016, 01:32 AM

Originally Posted by hugohonda

Hi Rick,

first thanks for sharing your code, it has been very useful to me. Similar to theCloud, i work with aircraft center of gravity and needed a function similar to yours. However I had some problems:
1- when a point was on a vertical boundarie
2- when a point was on an horizontal boundarie
3- when comparing two Double variable (b = Ycoord)...so I replace Double by Single and it works for now...I might need to use Double in the future, we'll see how it goes

For point numbers 1 and 2, I would refer you to my Note #3. You should not be overly fixated on points lying on the boundary of your polygon. Too many things can affect whether a point actually lies on the boundary. VB's floating point math is not infinite in scale, so the numbers have a limited precision (15 decimal places) and each number in a calculation theoretically degrades the accuracy of the calculated number to the point that you can never be sure if the calculated coordinate really lies on the boundary or not. For example, let's assume point is actually on the boundary... what is the significance if that is considered inside or outside the polygon on a practical level. Given that a theoretical boundary line for a polygon has zero thickness, then if the point was a millionth, trillions, zillions or whatever further away, it would be considered outside the polygon... similarly, if the point were as millionth, trillions, zillions or whatever closer in, it would be considered inside the polygon. At a practical level, either of those points (the further out one or the closer in one) could be considered the same point (what is the practical significance of two points that are a trillions of a inch, millimeter or whatever units of measurement?) So whether a real boundary point calculates inside or outside of the polygon should have absolutely no real world significance at all.

As for changing Doubles to Singles... I would not do that as it lessens the accuracy of all calculations that take place in the routine increasing the chance that a point in the vicinity of the boundary, but not truly "next to" or "on" it might calculate to the wrong side. Rather than making Doubles into Singles, leave them as Doubles and when comparing them, round each of the values down to the same number of decimal place (you choose the number of decimal places that is significant to you) and then compare those rounded down values instead.

**hugohonda** · 09-22-2016, 08:44 PM

Originally Posted by Rick Rothstein

For point numbers 1 and 2, I would refer you to my Note #3. You should not be overly fixated on points lying on the boundary of your polygon. Too many things can affect whether a point actually lies on the boundary. VB's floating point math is not infinite in scale, so the numbers have a limited precision (15 decimal places) and each number in a calculation theoretically degrades the accuracy of the calculated number to the point that you can never be sure if the calculated coordinate really lies on the boundary or not. For example, let's assume point is actually on the boundary... what is the significance if that is considered inside or outside the polygon on a practical level. Given that a theoretical boundary line for a polygon has zero thickness, then if the point was a millionth, trillions, zillions or whatever further away, it would be considered outside the polygon... similarly, if the point were as millionth, trillions, zillions or whatever closer in, it would be considered inside the polygon. At a practical level, either of those points (the further out one or the closer in one) could be considered the same point (what is the practical significance of two points that are a trillions of a inch, millimeter or whatever units of measurement?) So whether a real boundary point calculates inside or outside of the polygon should have absolutely no real world significance at all.

As for changing Doubles to Singles... I would not do that as it lessens the accuracy of all calculations that take place in the routine increasing the chance that a point in the vicinity of the boundary, but not truly "next to" or "on" it might calculate to the wrong side. Rather than making Doubles into Singles, leave them as Doubles and when comparing them, round each of the values down to the same number of decimal place (you choose the number of decimal places that is significant to you) and then compare those rounded down values instead.

Hi Rick,

first, thanks for your quick reply. I like your idea with using the round function with Doubles...I implemented it and it works great now.

As for my point numbers 1 and 2, I understand everything you mentioned and you are absolutely right. The thing is, I'm plotting different weight (y axis) and center of gravity (x axis) of an aircraft into a "Flight Envelope" and all point needs to be inside the envelope otherwise it isn't safe to fly. In order to find if a specific aircraft configuration (weight and CG) is dangerous to fly, I had to modify the code as mentioned before.

I know that the specific aircraft configuration touching the boundary, theoretically, won't be exactly on the line because like you mentioned, "a line has no thickness" but it will give me, I think, a very accurate aircraft configuration at which it is dangerous to fly.

I don't know if i'm clear but this is why i'm concerned about boundary.

I tested my code with many case scenarios and all of them give me a good result so far

thanks again!

**Rick Rothstein** · 09-22-2016, 09:20 PM

Originally Posted by hugohonda

As for my point numbers 1 and 2, I understand everything you mentioned and you are absolutely right. The thing is, I'm plotting different weight (y axis) and center of gravity (x axis) of an aircraft into a "Flight Envelope" and all point needs to be inside the envelope otherwise it isn't safe to fly. In order to find if a specific aircraft configuration (weight and CG) is dangerous to fly, I had to modify the code as mentioned before.

**Fikile Kentane** · 01-27-2017, 02:32 PM

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** · 01-27-2017, 03:56 PM

Originally Posted by 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

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** · 01-30-2017, 01:27 AM

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** · 02-01-2017, 02:43 PM

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