Android Question Counting geopoints in map area

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Using OSM_Droid, version 4.1.
I have a list of geopoints and want to count the number of those geopoints that are in a certain map area.
I can do this no problem if the area is rectangular and parallel to a latitude line, so for example the 4 corners are lat 52, lon 1 -- lat 52, lon 2 -- lat 51, lon 2, lat 51, lon 1.
In that case I can simply do (in SQL):
select count(*) from addresses where latitude < ? And latitude > ? And longitude > ? And longitude < ?

Can I do this though is the area is no rectangular and/or not parallel to a latitude line, say the area is for example a circle? I know there is SQLite Spatialite and I take it that can handle this, but I need encryption, so using SQLCipher.

RBS
 

MarkusR

Well-Known Member
Licensed User
Longtime User
say the area is for example a circle?
maybe something like this, means one x,y is the position at center of your circle area, the other x,y is your position, d = delta / difference
B4X:
dx=x2-x1
dy=y2-y1
length=sqrt(dx*dx + dy*dy)
if length<= radius
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
maybe something like this, means one x,y is the position at center of your circle area, the other x,y is your position, d = delta / difference
B4X:
dx=x2-x1
dy=y2-y1
length=sqrt(dx*dx + dy*dy)
if length<= radius

Yes, circle was the wrong example as it is simple mathematics. How about for example a complex polygon though?
I found some code actually that may do this and will see if I can convert this to B4A.

RBS
 
Upvote 0

MarkusR

Well-Known Member
Licensed User
Longtime User
keyword from ray tracing
intersection point polygon
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
No, didn't see that one, but I have a different function, which I will test later.

RBS

This one seems to work well with good performance, but tested no other ones yet.

B4X:
Sub GeoPointInPolygon(dPointLat As Double, dPointLon As Double, arrPolygonLatLon() As OSMDroid_GeoPoint ) As Boolean

    'http://www.mrexcel.com/forum/excel-questions/713005-does-point-fall-within-polygon-visual-basic-applications-function.html

    Dim i As Int
    Dim j As Int
    Dim PolySides As Int
    Dim OddNodes As Boolean

    PolySides = arrPolygonLatLon.Length - 1
    j = PolySides

 For i = 0 To PolySides
  If (((arrPolygonLatLon(i).Longitude < dPointLon And arrPolygonLatLon(j).Longitude >= dPointLon) _
             Or (arrPolygonLatLon(j).Longitude < dPointLon And arrPolygonLatLon(i).Longitude >= dPointLon)) _
             And (arrPolygonLatLon(i).Latitude <= dPointLat Or arrPolygonLatLon(j).Latitude <= dPointLat)) Then
   OddNodes = OddNodes <> (arrPolygonLatLon(i).Latitude + _
                                    (dPointLon - arrPolygonLatLon(i).Longitude) / _
                                    (arrPolygonLatLon(j).Longitude - arrPolygonLatLon(i).Longitude) * _
                                    (arrPolygonLatLon(j).Latitude - arrPolygonLatLon(i).Latitude) < dPointLat)
  End If
  j = i
 Next

    Return OddNodes

End Sub


RBS
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Found cases where my posted function didn't work (false negatives) and using your function now (after minor edit)
and that works indeed fine. Thanks for posting that.

RBS

Seen the new function fail as well now, I think on areas with more that about 8 sides.
Will check some further and post some data.
Looks a difficult one to solve.

RBS
 
Upvote 0

klaus

Expert
Licensed User
Longtime User
I think on areas with more that about 8 sides.
Sounds strange, the test project in the link in my post #6 uses 14 sides.
The red crosses are specific points for checking.

upload_2018-6-5_16-10-35.png


Can you post some data making the routine to fail?
What changes did you make in the routine?
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Sounds strange, the test project in the link in my post #6 uses 14 sides.
The red crosses are specific points for checking.

View attachment 68628

Can you post some data making the routine to fail?
What changes did you make in the routine?

I noticed your check points are all on the edge of the area. Did you check inside the area?

This is some data that failed:

52.605686 -2.16525
52.605849 -2.164698
52.606191 -2.164875
52.606833 -2.163662
52.607768 -2.162176
52.607618 -2.161898
52.607315 -2.161989
52.607113 -2.162337
52.6068 -2.162235
52.6065 -2.162482
52.6061 -2.16311
52.605829 -2.163711
52.605546 -2.164424

52.6071468 -2.1647307

Last row is the point inside the area to check.
Latitude in first column, Longitude in second.

There are no relevant changes, but this is my code:

B4X:
Sub CheckPointInPolygone(dPointLat As Double, dPointLon As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean
 Dim x1 As Double
 Dim y1 As Double
 Dim x2 As Double
 Dim y2 As Double
 Dim D As Double
 Dim i As Int
 Dim ni As Int
 Dim iPoints As Int
 
 ni = 0
 x1 = arrPolygonLatLon(0).Latitude
 y1 = arrPolygonLatLon(0).Longitude
 iPoints = arrPolygonLatLon.Length
 For i = 1 To iPoints
  If i < iPoints Then
   x2 = arrPolygonLatLon(i).Latitude
   y2 = arrPolygonLatLon(i).Longitude
  Else
   x2 = arrPolygonLatLon(0).Latitude ' checks the last line
   y2 = arrPolygonLatLon(0).Longitude
  End If
       
  If dPointLon >= Min(y1, y2) Then
   If dPointLon <= Max(y1, y2) Then
    If dPointLat <= Max(x1, x2) Then
     If (dPointLat = x1 And dPointLon = y1) Or (dPointLat = x1 And dPointLat = x2) Then ' checks vertices and vertical lines
      Return True
     End If
     If y1 <> y2 Then
      D = (dPointLon - y1) * (x2 - x1) / (y2 - y1) + x1
      If x1 = x2 Or dPointLat <= D Then
       ni = ni + 1
      End If
     End If
    End If
   End If
  End If
  x1 = x2
  y1 = y2
 Next
 
 If ni Mod 2 = 0 Then
  Return False
 Else
  Return True
 End If
End Sub


RBS
 
Upvote 0

klaus

Expert
Licensed User
Longtime User
Did you check inside the area?
Of course YES!
If you test the demo project you can move the finger and see if you are inside or outside the polygon.
I noticed your check points are all on the edge of the area.
The critical points are always on the edges!

I have never tested with lat/lon points.
The algorithm works only with a closed polygon.
I will look at it tomorrow with your data.
 
Last edited:
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Of course YES!
If you test the demo project you can move the finger and see if you are inside or outside the polygon.

The critical points are always on the edges!

I have never tested with lat/lon points.
The algorithm works only with a closed polygon.
I will look at it tomorrow with your data.

If it fails it always seems to fail on all given points, so iow it's the polygon that causes the problems.
I had a look at another algo:
http://geomalgorithms.com/a03-_inclusion.html
Using a different approach.
Converted to VBA and the B4A, but also fails on occasions:

B4X:
Sub WindingPointInPolygon(dLatitude As Double, dLongitude As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean

    Dim i As Int
    Dim wn As Int
    Dim iUBound As Int
 
 iUBound = arrPolygonLatLon.Length - 1

 For i = 0 To iUBound - 1 '- 1 as we do i + 1
  If arrPolygonLatLon(i).Longitude <= dLongitude Then
   If arrPolygonLatLon(i + 1).Longitude > dLongitude Then
    If IsLeft(arrPolygonLatLon(i), arrPolygonLatLon(i + 1), dLatitude, dLongitude) Then
     wn = wn + 1
    End If
   End If
  Else
   If arrPolygonLatLon(i + 1).Longitude <= dLongitude Then
    If IsLeft(arrPolygonLatLon(i), arrPolygonLatLon(i + 1), dLatitude, dLongitude) = False Then
     wn = wn - 1
    End If
   End If
  End If
 Next
    
    Return wn <> 0

End Sub

Sub IsLeft(P0 As OSMDroid_GeoPoint, P1 As OSMDroid_GeoPoint, dLatitude As Double, dLongitude As Double) As Boolean
 
 Dim dReturn As Double

    dReturn = (P1.Latitude - P0.Latitude) * (dLongitude - P0.Longitude) - _
              (dLatitude - P0.Latitude) * (P1.Longitude - P0.Longitude)
     
 Return dReturn > 0

End Sub

This code needs the path to be closed and I think also the last line crossing the first.


RBS
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Of course YES!
If you test the demo project you can move the finger and see if you are inside or outside the polygon.

The critical points are always on the edges!

I have never tested with lat/lon points.
The algorithm works only with a closed polygon.
I will look at it tomorrow with your data.

Found a bug somewhere else and your function works OK, so don't look at this.
I have tested 4 different algo's now and will do a bit of testing and see which one is fastest
and/or most accurate and will report that back.
Apologies for causing confusion.

RBS
 
Upvote 0

RB Smissaert

Well-Known Member
Licensed User
Longtime User
Found a bug somewhere else and your function works OK, so don't look at this.
I have tested 4 different algo's now and will do a bit of testing and see which one is fastest
and/or most accurate and will report that back.
Apologies for causing confusion.

RBS

Tested the 4 algo's I have used now and all performing well and fast with similar times.
Not tested unusual polygons, eg with holes in it. All giving the same result (mappoints found).
Tested with 12135 map points and 96 polygon points, on a Samsung S9, debugging off.
Code posted below. I think 3 of those may be similar, but WindingPointInPolygon is based on a different
principle. Times in average milliseconds run on the above sample 30 times.

GeoPointInPolygon 256
WindingPointInPolygon 194
PtInPoly 230
CheckPointInPolygone 243

The author of WindingPointInPolygon (see URL) claims it is more accurate and as it is the fastest in my test
I am using that one for now.

B4X:
Sub GeoPointInPolygon(dPointLat As Double, dPointLon As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean

 'http://www.mrexcel.com/forum/excel-questions/713005-does-point-fall-within-polygon-visual-basic-applications-function.html
 '--------------------------------------------------------------------------------------------------------------------------

 Dim i As Int
 Dim j As Int
 Dim PolySides As Int
 Dim OddNodes As Boolean

 PolySides = arrPolygonLatLon.Length - 1
 j = PolySides

 For i = 0 To PolySides
  If (((arrPolygonLatLon(i).Longitude < dPointLon And arrPolygonLatLon(j).Longitude >= dPointLon) _
             Or (arrPolygonLatLon(j).Longitude < dPointLon And arrPolygonLatLon(i).Longitude >= dPointLon)) _
             And (arrPolygonLatLon(i).Latitude <= dPointLat Or arrPolygonLatLon(j).Latitude <= dPointLat)) Then
   OddNodes = OddNodes <> (arrPolygonLatLon(i).Latitude + _
                                    (dPointLon - arrPolygonLatLon(i).Longitude) / _
                                    (arrPolygonLatLon(j).Longitude - arrPolygonLatLon(i).Longitude) * _
                                    (arrPolygonLatLon(j).Latitude - arrPolygonLatLon(i).Latitude) < dPointLat)
  End If
  j = i
 Next

 Return OddNodes

End Sub

'------------------------------------------------------------------------

Sub WindingPointInPolygon(dLatitude As Double, dLongitude As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean
 
 'http://geomalgorithms.com/a03-_inclusion.html
 '2012 Dan Sunday
 '--------------------------------------------

    Dim i As Int
    Dim wn As Int
    Dim iUBound As Int
 
 iUBound = arrPolygonLatLon.Length - 1

 For i = 0 To iUBound - 1 '- 1 as we do i + 1
  If arrPolygonLatLon(i).Longitude <= dLongitude Then
   If arrPolygonLatLon(i + 1).Longitude > dLongitude Then
    If IsLeft(arrPolygonLatLon(i), arrPolygonLatLon(i + 1), dLatitude, dLongitude) Then
     wn = wn + 1
    End If
   End If
  Else
   If arrPolygonLatLon(i + 1).Longitude <= dLongitude Then
    If IsLeft(arrPolygonLatLon(i), arrPolygonLatLon(i + 1), dLatitude, dLongitude) = False Then
     wn = wn - 1
    End If
   End If
  End If
 Next
    
    Return wn <> 0

End Sub

Sub IsLeft(P0 As OSMDroid_GeoPoint, P1 As OSMDroid_GeoPoint, dLatitude As Double, dLongitude As Double) As Boolean
 
 Return (P1.Latitude - P0.Latitude) * (dLongitude - P0.Longitude) - _
           (dLatitude - P0.Latitude) * (P1.Longitude - P0.Longitude) > 0
     
End Sub

'----------------------------------------------------------------------------------------------

Sub PtInPoly(dLatitude As Double, dLongitude As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean

    Dim i As Int
    Dim NumSidesCrossed As Int
    Dim m As Double
    Dim b As Double
    Dim Boolean1 As Boolean
    Dim Boolean2 As Boolean
 Dim iUBound As Int
 
 'Rick RothStein, http://www.excelfox.com/forum/showthread.php/1579-Test-Whether-A-Point-Is-In-A-Polygon-Or-Not
 '-------------------------------------------------------------------------------------------------------------
 
 iUBound = arrPolygonLatLon.Length - 1

    'Xor in this case is true if expression1 and expression2 are: True False or False True
    '-------------------------------------------------------------------------------------
    For i = 0 To iUBound - 1
        Boolean1 = arrPolygonLatLon(i).Latitude > dLatitude
        Boolean2 = arrPolygonLatLon(i + 1).Latitude > dLatitude
        If Boolean1 <> Boolean2 Then
            m = (arrPolygonLatLon(i + 1).Longitude - arrPolygonLatLon(i).Longitude) / (arrPolygonLatLon(i + 1).Latitude - arrPolygonLatLon(i).Latitude)
            b = (arrPolygonLatLon(i).Longitude * arrPolygonLatLon(i + 1).Latitude - arrPolygonLatLon(i).Latitude * arrPolygonLatLon(i + 1).Longitude) / (arrPolygonLatLon(i + 1).Latitude - arrPolygonLatLon(i).Latitude)
            If m * dLatitude + b > dLongitude Then
                NumSidesCrossed = NumSidesCrossed + 1
            End If
        End If
    Next

    Return (NumSidesCrossed = 1)

End Sub

'-------------------------------------------------------------------------------------------------

Sub CheckPointInPolygone(dPointLat As Double, dPointLon As Double, arrPolygonLatLon() As OSMDroid_GeoPoint) As Boolean
 Dim x1 As Double
 Dim y1 As Double
 Dim x2 As Double
 Dim y2 As Double
 Dim D As Double
 Dim i As Int
 Dim ni As Int
 Dim iPoints As Int
 'Klaus
 '-------------
 
 ni = 0
 x1 = arrPolygonLatLon(0).Latitude
 y1 = arrPolygonLatLon(0).Longitude
 iPoints = arrPolygonLatLon.Length
 For i = 1 To iPoints
  If i < iPoints Then
   x2 = arrPolygonLatLon(i).Latitude
   y2 = arrPolygonLatLon(i).Longitude
  Else
   x2 = arrPolygonLatLon(0).Latitude ' checks the last line
   y2 = arrPolygonLatLon(0).Longitude
  End If
       
  If dPointLon >= Min(y1, y2) Then
   If dPointLon <= Max(y1, y2) Then
    If dPointLat <= Max(x1, x2) Then
     If (dPointLat = x1 And dPointLon = y1) Or (dPointLat = x1 And dPointLat = x2) Then ' checks vertices and vertical lines
      Return True
     End If
     If y1 <> y2 Then
      D = (dPointLon - y1) * (x2 - x1) / (y2 - y1) + x1
      If x1 = x2 Or dPointLat <= D Then
       ni = ni + 1
      End If
     End If
    End If
   End If
  End If
  x1 = x2
  y1 = y2
 Next
 
 If ni Mod 2 = 0 Then
  Return False
 Else
  Return True
 End If
End Sub


RBS
 
Upvote 0
Top