Android Question area of polygon

Discussion in 'Android Questions' started by RB Smissaert, Jun 6, 2019.

  1. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Needed a way to get the centre of a polygon, defined by a list of osmdroid_geopoints.
    Part of this algorithm is to get the area of the polygon. All this works fine and I get the
    centre of the polygon (defined as making the polygon balance on a pin if put at the centre)
    Having this area I would like to know how to get the real are in square km or whatever unit.
    The actual obtained area value is very small, so there need to be some multiplication factor, but
    not sure how to get this factor.

    Code:
    Sub GetPolygonArea(lstPoints As ListAs Double

        
    Dim i As Int
        
    Dim dArea As Double

        
    'Get the areas.
        For i = 0 To lstPoints.Size - 2
      
    Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(i)
      
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i + 1)
            dArea = dArea + _
                    (oGeoPoint2.Longitude - oGeoPoint1.Longitude) * _
                    (oGeoPoint2.Latitude + oGeoPoint1.Latitude) / 
    2
     
    Next

     
    Return Abs(dArea)

    End Sub

    Sub FindCentroid(lstPoints As ListAs OSMDroid_GeoPoint

        
    Dim i As Long
        
    Dim second_factor As Double
        
    Dim polygon_area As Double
        
    Dim X As Double
        
    Dim Y As Double
     
    Dim oGeoPointCenter As OSMDroid_GeoPoint

        
    For i = 0 To lstPoints.Size - 2
      
    Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(i)
      
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i + 1)
            second_factor = _
            oGeoPoint1.Longitude * oGeoPoint2.Latitude - _
                            oGeoPoint2.Longitude * oGeoPoint1.Latitude
            X = X + (oGeoPoint1.Longitude + oGeoPoint2.Longitude) * _
                second_factor
            Y = Y + (oGeoPoint1.Latitude + oGeoPoint2.Latitude) * _
                second_factor
        
    Next

        
    'Divide by 6 times the polygon's area.
        polygon_area = GetPolygonArea(lstPoints)
     General.RunLog(
    "FindCentroid, polygon_area: " & polygon_area)
     
        X = X / 
    6 / polygon_area
        Y = Y / 
    6 / polygon_area
     

        
    'If the values are negative, the polygon is
        'oriented counterclockwise. Reverse the signs.
     '(this works for me, but maybe wrong for other locations)
     If Y < 0 Then
      X = -X
      Y = -Y
     
    End If
     
     oGeoPointCenter.Initialize(Y, X)

     
    Return oGeoPointCenter

    End Sub

    RBS
     
  2. klaus

    klaus Expert Licensed User

    Routine to calculate the area in squaremeters.

    Code:
    Sub GetPolygonArea(lstPoints As ListAs Double
        
    Dim i As Int
        
    Dim Lat1, Lat2, Lng1, Lng2 As Double
        
    Dim dArea, LatScale, LngScale As Double

        
    'Get the areas.
        Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(0)
        LatScale =     
    cPI * 6371000 / 180        ' converts latitude degrees to m
        LngScale = LatScale * CosD(oGeoPoint1.Latitude) ' converts longitude degrees to m
        Lat1 = oGeoPoint1.Latitude * LatScale
        Lng1 = oGeoPoint1.Longitude * LngScale
        
    For i = 1 To lstPoints.Size - 1
            
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i)
            Lat2 = oGeoPoint2.Latitude * LatScale
            Lng2 = oGeoPoint2.Longitude * LngScale
            dArea = dArea + (Lng2 - Lng1) * (Lat2 + Lat1) / 
    2
            Lat1 = Lat2
            Lng1 = Lng2
        
    Next
       
        
    Return Abs(dArea)
    End Sub
     
    Last edited: Jun 7, 2019
    rboeck, DonManfred and Erel like this.
  3. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Thanks, will try that and report back.

    RBS
     
  4. derez

    derez Expert Licensed User

  5. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

  6. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Shouldn't your 6317000 be 6378137, according to:
    https://en.wikipedia.org/wiki/World_Geodetic_System ?

    RBS
     
  7. klaus

    klaus Expert Licensed User

    Sorry, it should be 6371000, I permuted the 1 and 7.
    Amended in post#2.
    I got the value from HERE, 6371 km is a kind of mean radius.
    6378137 is the semi-major axis and 6356752 is the semi-minor axis.
     
  8. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    The area I am interested in is only small, latitude from about 52.45 to 52.65, so I will take the mean radius
    of that and use that for all locations.

    RBS
     
  9. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    nav.SphericalArea2 doesn't seem to work with squares or rectangles, gives me zero:

    Code:
    Sub GetPolygonAreaSpherical(lstPoints As ListAs Double

     
    Dim i As Int
     
    Dim iPoints As Int
     
    Dim dArea As Double
     
     iPoints = lstPoints.Size - 
    1
     
     
    Dim arrPoints(iPoints, 2As Double
     
     
    For i = 0 To iPoints - 1 'so don't include last point, which is same as first
      Dim oGeoPoint As OSMDroid_GeoPoint = lstPoints.Get(i)
      arrPoints(i, 
    0) = oGeoPoint.Latitude
      arrPoints(i, 
    1) = oGeoPoint.Longitude
     
    Next

     dArea = nav.SphericalArea2(arrPoints, iPoints, 
    01)
     
     
    Return Round2(Abs(dArea), 4)

    End Sub
    Running the above Sub on a non-rectangle works fine.
    Am I overlooking something?


    RBS
     
  10. derez

    derez Expert Licensed User

    Please check your data. I run this sub and get 150.0697679885843

    Code:
    Sub test2
    Dim points(4,2As Double
        points(
    0,0) = 52.45
        points(
    0,1) = 30.10
        points(
    1,0) = 52.65
        points(
    1,1) = 30.10
        points(
    2,0) = 52.65
        points(
    2,1) = 30.20
        points(
    3,0) = 52.45
        points(
    3,1) = 30.20
        
    Log(nav.SphericalArea2(points,4,0,1))
    End Sub
     
  11. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Try this:

    Code:
    Sub test2
     
    Dim points(4,2As Double
     
    Dim dArea As Double
     points(
    0,0) = 52.597316
     points(
    0,1) = -2.150112
     points(
    1,0) = 52.597316
     points(
    1,1) = -2.147188
     points(
    2,0) = 52.594396
     points(
    2,1) = -2.147188
     points(
    3,0) = 52.594396
     points(
    3,1) = -2.150112
     dArea = nav.SphericalArea2(points,
    4,0,1)
     
    Log(dArea)
     
    Log(Round2(dArea, 4))
    End Sub
    Gives me NaN and 0


    RBS
     
  12. derez

    derez Expert Licensed User

    The area enclosed is less than 1 square km, so you should not use this method here. I decreased the distance slowly and got down to area of 0.2 , I guess the bearing calculation in the algorithm cannot handle it.
     
  13. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    So, I should then use FlatArea2 first and see if the area is > 1 km2 and then use SpericalArea2?
    Although it might a bit less accurate, in that case it makes sense perhaps to use the algo posted by Klaus.

    There also seems to be a problem with SphericalArea2 if there are a lots of points:

    Code:
    Sub test2
     
    Dim points(72,2As Double
     
    Dim dArea As Double
     points(
    0,0) = 52.613914
     points(
    0,1) = -2.20366
     points(
    1,0) = 52.614383
     points(
    1,1) = -2.2019
     points(
    2,0) = 52.615191
     points(
    2,1) = -2.200913
     points(
    3,0) = 52.618057
     points(
    3,1) = -2.200012
     points(
    4,0) = 52.619073
     points(
    4,1) = -2.198381
     points(
    5,0) = 52.620037
     points(
    5,1) = -2.197265
     points(
    6,0) = 52.619933
     points(
    6,1) = -2.181859
     points(
    7,0) = 52.617719
     points(
    7,1) = -2.179069
     points(
    8,0) = 52.616911
     points(
    8,1) = -2.177224
     points(
    9,0) = 52.616442
     points(
    9,1) = -2.174949
     points(
    10,0) = 52.616259
     points(
    10,1) = -2.173275
     points(
    11,0) = 52.617588
     points(
    11,1) = -2.171688
     points(
    12,0) = 52.619047
     points(
    12,1) = -2.170572
     points(
    13,0) = 52.620506
     points(
    13,1) = -2.174091
     points(
    14,0) = 52.623633
     points(
    14,1) = -2.176408
     points(
    15,0) = 52.624362
     points(
    15,1) = -2.174434
     points(
    16,0) = 52.626967
     points(
    16,1) = -2.174649
     points(
    17,0) = 52.62741
     points(
    17,1) = -2.171344
     points(
    18,0) = 52.628686
     points(
    18,1) = -2.170958
     points(
    19,0) = 52.631161
     points(
    19,1) = -2.171044
     points(
    20,0) = 52.631812
     points(
    20,1) = -2.172203
     points(
    21,0) = 52.631578
     points(
    21,1) = -2.174992
     points(
    22,0) = 52.630796
     points(
    22,1) = -2.176194
     points(
    23,0) = 52.630796
     points(
    23,1) = -2.176194
     points(
    24,0) = 52.631526
     points(
    24,1) = -2.181429
     points(
    25,0) = 52.632698
     points(
    25,1) = -2.181816
     points(
    26,0) = 52.635511
     points(
    26,1) = -2.184047
     points(
    27,0) = 52.635849
     points(
    27,1) = -2.18585
     points(
    28,0) = 52.635068
     points(
    28,1) = -2.187695
     points(
    29,0) = 52.634339
     points(
    29,1) = -2.191171
     points(
    30,0) = 52.634339
     points(
    30,1) = -2.191171
     points(
    31,0) = 52.634573
     points(
    31,1) = -2.193188
     points(
    32,0) = 52.636448
     points(
    32,1) = -2.193918
     points(
    33,0) = 52.637829
     points(
    33,1) = -2.194862
     points(
    34,0) = 52.637959
     points(
    34,1) = -2.197008
     points(
    35,0) = 52.637438
     points(
    35,1) = -2.199625
     points(
    36,0) = 52.636136
     points(
    36,1) = -2.200269
     points(
    37,0) = 52.634339
     points(
    37,1) = -2.200312
     points(
    38,0) = 52.633193
     points(
    38,1) = -2.201428
     points(
    39,0) = 52.632724
     points(
    39,1) = -2.202887
     points(
    40,0) = 52.633193
     points(
    40,1) = -2.205934
     points(
    41,0) = 52.633375
     points(
    41,1) = -2.208251
     points(
    42,0) = 52.634964
     points(
    42,1) = -2.211728
     points(
    43,0) = 52.635328
     points(
    43,1) = -2.213487
     points(
    44,0) = 52.633974
     points(
    44,1) = -2.214903
     points(
    45,0) = 52.633036
     points(
    45,1) = -2.215032
     points(
    46,0) = 52.63288
     points(
    46,1) = -2.213873
     points(
    47,0) = 52.631838
     points(
    47,1) = -2.209753
     points(
    48,0) = 52.631421
     points(
    48,1) = -2.209067
     points(
    49,0) = 52.629911
     points(
    49,1) = -2.209024
     points(
    50,0) = 52.629129
     points(
    50,1) = -2.209367
     points(
    51,0) = 52.628556
     points(
    51,1) = -2.208681
     points(
    52,0) = 52.627306
     points(
    52,1) = -2.204303
     points(
    53,0) = 52.626394
     points(
    53,1) = -2.206707
     points(
    54,0) = 52.626915
     points(
    54,1) = -2.209239
     points(
    55,0) = 52.626759
     points(
    55,1) = -2.21044
     points(
    56,0) = 52.62556
     points(
    56,1) = -2.210569
     points(
    57,0) = 52.625665
     points(
    57,1) = -2.211899
     points(
    58,0) = 52.625665
     points(
    58,1) = -2.211899
     points(
    59,0) = 52.626576
     points(
    59,1) = -2.21323
     points(
    60,0) = 52.626238
     points(
    60,1) = -2.215075
     points(
    61,0) = 52.624883
     points(
    61,1) = -2.214946
     points(
    62,0) = 52.622617
     points(
    62,1) = -2.215161
     points(
    63,0) = 52.621887
     points(
    63,1) = -2.21486
     points(
    64,0) = 52.620532
     points(
    64,1) = -2.214517
     points(
    65,0) = 52.619464
     points(
    65,1) = -2.213916
     points(
    66,0) = 52.619829
     points(
    66,1) = -2.210226
     points(
    67,0) = 52.620975
     points(
    67,1) = -2.206792
     points(
    68,0) = 52.621575
     points(
    68,1) = -2.202973
     points(
    69,0) = 52.61949
     points(
    69,1) = -2.202286
     points(
    70,0) = 52.617041
     points(
    70,1) = -2.202973
     points(
    71,0) = 52.615113
     points(
    71,1) = -2.204217
     dArea = nav.SphericalArea2(points,
    72,0,1)
     
    Log(dArea)
     
    Log(Round2(dArea, 4))
    End Sub

    RBS
     
  14. emexes

    emexes Well-Known Member Licensed User

    What areas are being measuring, and for what purpose? If the data above is typical, and the areas are less than a square degree, it's probable that flat areas (and straight-line distances) are accurate enough for your purpose.

    If they're not accurate enough, then perhaps revisit the conversion of angles (longitude and latitude) to linear distances, because using an average radius is going to introduce errors up to the order of +/- 0.3% (for distance, and twice that for areas).

    On the bright side, if you are just comparing areas or distances to see if one is smaller or closer than the other, then most of that error will cancel out anyway. Bonus!
     
  15. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    It doesn't have to be seriously accurate and I am not comparing area's.
    My range of latitude is very small, so one radius for all locations is fine.

    To me it looks there is a bug in nav.SpericalArea2 if there are a lot of points.

    RBS
     
  16. emexes

    emexes Well-Known Member Licensed User

    I think you're right. My first guess would be a float (single-precision) operation being used with small-difference-of-large-number maths.

    Bonus!

    I stayed out of the radius-of-earth discussion, but one useful factoid I use all the time when close-enough is good-enough, is that the metre was originally defined as 10,000,000 metres = 90 degrees of earth circumference (from north pole to equator, through Paris). Thus for converting degrees latitude to distance you can use:

    metres = 10000000 * latitude / 90

    no Pi or hard-to-remember earth radius required :)
     
  17. emexes

    emexes Well-Known Member Licensed User

    I did some tests, using a 1x1 degree patch near the equator, and:
    1/ same problem occurs for just 4 points
    2/ problem randomly comes and goes as patch is slowly enlarged from 1x1 to 1x2 degrees
    3/ problem seems to be with the latitude
    4/ when it does return results, the area should be increasing monotonically as the NS latitude size of the patch increases, but... check this out:
    Code:
    latitude size         area

    2.000000000152  NaN
    2.000000000153  24749.007912058205
    2.000000000154  24749.00791209434
    2.000000000155  24749.007912166602
    2.000000000156  24749.00791209434
    2.000000000157  24749.007912166602
    2.000000000158  24749.00791213047
    2.000000000159  24749.00791213047
    2.000000000160  24749.007912166602
    2.000000000161  24749.007912166602
    2.000000000162  24749.00791209434
    2.000000000163  24749.00791213047
    2.000000000164  24749.007912058205
    2.000000000165  24749.00791213047
    2.000000000166  24749.00791209434
    2.000000000167  24749.00791209434
    2.000000000168  24749.007912166602
    2.000000000169  NaN
    2.000000000170  24749.00791220273
    2.000000000171  NaN
    2.000000000172  24749.007912274992
    2.000000000173  NaN
    Hmm. Handy to know. Library author/creator Derez is still listed as a licenced user and active forum member; I'll let you have the honour of giving him the good news (if you haven't already).

    :)
     
  18. emexes

    emexes Well-Known Member Licensed User

    Had a bit of an epiphany over lunch, modified your code to randomly shift the latitudes up/down by ~1/10,000,000th of a degree (ie < 1 cm), figured that if each of those nudged 72 points has a 50% chance of landing in a "good place" then it should only take about 2^72 tests to find a set where all points are good and the spherical area is returned... at say 100,000 tests per second, that should only take... (grabs calculator) 1,496,427,638 years... hmm...

    No worries, tried it anyway. That lead me to discovering that your polygon has three duplicate points, which I removed and then... voila! comes back with area that seems about right. I still think there's a problem in the library though, given that my simple 4-point patches did not contain duplicate points.
    Code:
    ** Activity (main) Create, isFirst = true **
    4.409755228643426
    4.4098
    ** 
    Activity (main) Resume **
    Code:
    Dim points(722As Double
    Dim dArea As Double

    points(
    0,0) = 52.613914
    points(
    0,1) = -2.20366
    points(
    1,0) = 52.614383
    points(
    1,1) = -2.2019
    points(
    2,0) = 52.615191
    points(
    2,1) = -2.200913
    points(
    3,0) = 52.618057
    points(
    3,1) = -2.200012
    points(
    4,0) = 52.619073
    points(
    4,1) = -2.198381
    points(
    5,0) = 52.620037
    points(
    5,1) = -2.197265
    points(
    6,0) = 52.619933
    points(
    6,1) = -2.181859
    points(
    7,0) = 52.617719
    points(
    7,1) = -2.179069
    points(
    8,0) = 52.616911
    points(
    8,1) = -2.177224
    points(
    9,0) = 52.616442
    points(
    9,1) = -2.174949
    points(
    10,0) = 52.616259
    points(
    10,1) = -2.173275
    points(
    11,0) = 52.617588
    points(
    11,1) = -2.171688
    points(
    12,0) = 52.619047
    points(
    12,1) = -2.170572
    points(
    13,0) = 52.620506
    points(
    13,1) = -2.174091
    points(
    14,0) = 52.623633
    points(
    14,1) = -2.176408
    points(
    15,0) = 52.624362
    points(
    15,1) = -2.174434
    points(
    16,0) = 52.626967
    points(
    16,1) = -2.174649
    points(
    17,0) = 52.62741
    points(
    17,1) = -2.171344
    points(
    18,0) = 52.628686
    points(
    18,1) = -2.170958
    points(
    19,0) = 52.631161
    points(
    19,1) = -2.171044
    points(
    20,0) = 52.631812
    points(
    20,1) = -2.172203
    points(
    21,0) = 52.631578
    points(
    21,1) = -2.174992
    points(
    22,0) = 52.630796
    points(
    22,1) = -2.176194
    'same as previous point
    'points(23,0) = 52.630796
    'points(23,1) = -2.176194
    points(24-1,0)= 52.631526
    points(
    24-1,1)= -2.181429
    points(
    25-1,0)= 52.632698
    points(
    25-1,1)= -2.181816
    points(
    26-1,0)= 52.635511
    points(
    26-1,1)= -2.184047
    points(
    27-1,0)= 52.635849
    points(
    27-1,1)= -2.18585
    points(
    28-1,0)= 52.635068
    points(
    28-1,1)= -2.187695
    points(
    29-1,0)= 52.634339
    points(
    29-1,1)= -2.191171
    'same as previous point
    'points(30,0) = 52.634339
    'points(30,1) = -2.191171
    points(31-2,0)= 52.634573
    points(
    31-2,1)= -2.193188
    points(
    32-2,0)= 52.636448
    points(
    32-2,1)= -2.193918
    points(
    33-2,0)= 52.637829
    points(
    33-2,1)= -2.194862
    points(
    34-2,0)= 52.637959
    points(
    34-2,1)= -2.197008
    points(
    35-2,0)= 52.637438
    points(
    35-2,1)= -2.199625
    points(
    36-2,0)= 52.636136
    points(
    36-2,1)= -2.200269
    points(
    37-2,0)= 52.634339
    points(
    37-2,1)= -2.200312
    points(
    38-2,0)= 52.633193
    points(
    38-2,1)= -2.201428
    points(
    39-2,0)= 52.632724
    points(
    39-2,1)= -2.202887
    points(
    40-2,0)= 52.633193
    points(
    40-2,1)= -2.205934
    points(
    41-2,0)= 52.633375
    points(
    41-2,1)= -2.208251
    points(
    42-2,0)= 52.634964
    points(
    42-2,1)= -2.211728
    points(
    43-2,0)= 52.635328
    points(
    43-2,1)= -2.213487
    points(
    44-2,0)= 52.633974
    points(
    44-2,1)= -2.214903
    points(
    45-2,0)= 52.633036
    points(
    45-2,1)= -2.215032
    points(
    46-2,0)= 52.63288
    points(
    46-2,1)= -2.213873
    points(
    47-2,0)= 52.631838
    points(
    47-2,1)= -2.209753
    points(
    48-2,0)= 52.631421
    points(
    48-2,1)= -2.209067
    points(
    49-2,0)= 52.629911
    points(
    49-2,1)= -2.209024
    points(
    50-2,0)= 52.629129
    points(
    50-2,1)= -2.209367
    points(
    51-2,0)= 52.628556
    points(
    51-2,1)= -2.208681
    points(
    52-2,0)= 52.627306
    points(
    52-2,1)= -2.204303
    points(
    53-2,0)= 52.626394
    points(
    53-2,1)= -2.206707
    points(
    54-2,0)= 52.626915
    points(
    54-2,1)= -2.209239
    points(
    55-2,0)= 52.626759
    points(
    55-2,1) = -2.21044
    points(
    56-2,0) = 52.62556
    points(
    56-2,1) = -2.210569
    points(
    57-2,0) = 52.625665
    points(
    57-2,1) = -2.211899
    'same as previous point
    'points(58,0) = 52.625665
    'points(58,1) = -2.211899
    points(59-3,0) = 52.626576
    points(
    59-3,1) = -2.21323
    points(
    60-3,0) = 52.626238
    points(
    60-3,1) = -2.215075
    points(
    61-3,0) = 52.624883
    points(
    61-3,1) = -2.214946
    points(
    62-3,0) = 52.622617
    points(
    62-3,1) = -2.215161
    points(
    63-3,0) = 52.621887
    points(
    63-3,1) = -2.21486
    points(
    64-3,0) = 52.620532
    points(
    64-3,1) = -2.214517
    points(
    65-3,0) = 52.619464
    points(
    65-3,1) = -2.213916
    points(
    66-3,0) = 52.619829
    points(
    66-3,1) = -2.210226
    points(
    67-3,0) = 52.620975
    points(
    67-3,1) = -2.206792
    points(
    68-3,0) = 52.621575
    points(
    68-3,1) = -2.202973
    points(
    69-3,0) = 52.61949
    points(
    69-3,1) = -2.202286
    points(
    70-3,0) = 52.617041
    points(
    70-3,1) = -2.202973
    points(
    71-3,0) = 52.615113
    points(
    71-3,1) = -2.204217
       
    dArea = nav.SphericalArea2(points, 
    72-301)
    Log(dArea)
    Log(Round2(dArea, 4))
     
  19. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Thanks for spotting those duplicates.
    Should have checked.

    RBS
     
  20. emexes

    emexes Well-Known Member Licensed User

    No worries. Old Einstein here completely missed them too, even though I was on-the-ball enough to plot the polygon to make sure it didn't cross over itself. Spewin! upload_2019-6-9_14-4-50.png
     
Loading...
  1. This site uses cookies to help personalise content, tailor your experience and to keep you logged in if you register.
    By continuing to use this site, you are consenting to our use of cookies.
    Dismiss Notice