  # Android Question area of polygon

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

1. Thanks for testing, emexes ! The java source for the Surfacearea2 method is following, if you find anything wrong just point it out:
Code:
`/**      * The methods calculates the area on the Earth's spherical enclosed inside a polygon.      * n is the polygon's size.      * GindxLat and GindxLong are the indices of the parameters Lat and Long in the array.      * The returned area is in square Km.      */     public double SphericalArea2(double [][] Poligon, int n,             int GindxLat, int GindxLong)        {            double meanlat = 0;            double [][] points =  new double[n] ; // new array because the content must change            for (int i = 0;i<n;i++)            {                points[i] = Poligon[i][GindxLat]* deg2rad ; // allocation by the index                points[i] = Poligon[i][GindxLong]* deg2rad ;                meanlat += points[i];            }            meanlat = meanlat/n ;            return Spherical(points , n, 0, 1, meanlat);        }          private double Spherical(double [][] points , int n, int UindxX, int UindxY, double meanlat)       {             double sum1,sum2,Rs,phis  ;                     double es = 0.0067394967423 ;  // eccentricity            double twoPi = 2 * Math.PI ;            phis = Math.atan(Math.tan(meanlat)/(1+es));            Rs = 6378.137 / Math.pow((1+es*Math.sin(phis)*Math.sin(phis)),0.5);  // km            // find angles for point 0            double f1 = Bearing(points,points,points,points);            double f2 = Bearing(points,points,points[n-1],points[n-1]);            sum1 = (f1 + twoPi - f2)%twoPi;            sum2 = (f2 + twoPi - f1)%twoPi;            // find angles for point n-1            f1 = Bearing(points[n-1],points[n-1],points,points);            f2 = Bearing(points[n-1],points[n-1],points[n-2],points[n-2]) ;            sum1 += (f1 + twoPi - f2)%twoPi;            sum2 += (f2 + twoPi - f1)%twoPi;            // find angles for the intermediate points            for (int i = 1;i<n-1;i++)            {            f1 = Bearing(points[i],points[i],points[i+1],points[i+1]);            f2 = Bearing(points[i],points[i],points[i-1],points[i-1]);            sum1 += (f1 + twoPi - f2)%twoPi;            sum2 += (f2 + twoPi - f1)%twoPi;            }            sum1 =  Math.min(sum1,sum2) ; // the minimun is the internal angles sum            return Math.abs((Math.abs(sum1)-(n-2)*Math.PI)) * Rs*Rs ;// in sqr km       }        // This method uses arguments in Radians and returns Radians !        private double Bearing(double LatR1,double LongR1,                               double LatR2,double LongR2)        {            double  range,bearing;            double L = LongR1 - LongR2 ;            double D = (Math.sin(LatR1) * Math.sin(LatR2)) +                        (Math.cos(LatR1) * Math.cos(LatR2) *Math.cos(L ));            range = Math.acos(D);            double C = (Math.sin(LatR2) - (Math.sin(LatR1) * D)) /                        (Math.cos(LatR1) * Math.sin(range));            if (Math.sin(L) < 0)                 {bearing = Math.acos(C) ;}            else                 {bearing = (2*Math.PI - Math.acos(C))  ;}            return  bearing ;        }`
The code can be easily translated to enable b4a debug.

Last edited: Jun 9, 2019
2. That's what I thought, too... I've translated it, but it's coming back with massively huge areas which I think has something to do with this bit:
Code:
`1            sum1 += (f1 + twoPi - f2)%twoPi;2            sum2 += (f2 + twoPi - f1)%twoPi;3            }4            sum1 =  Math.min(sum1,sum2) ; // the minimun is the internal angles sum5            return Math.abs((Math.abs(sum1)-(n-2)*Math.PI)) * Rs*Rs ;// in sqr km`
where, as I understand it, on lines 1 and 2: sum1 and sum2 are constrained to 0..twoPi by the % mod operation,
and thus on line 4, sum1 must thus also similarly be within 0..twoPi,
but then on line 5 when we subtract (n-2)*Pi, I end up with a large number is then multiplied by the Rs squared

But I've been using nav.GeoDistanceBearing to calculate the bearing; now that you've posted the Bearing function, I'll use that instead... like perhaps I got the units or zero base angle wrong.

The only thing I've worked out thus far is that it doesn't handle repeated points, which is fair enough because it would be pretty difficult to calculate a bearing from one position to the same position.

3. which is probably this line:
Code:
`double C = (Math.sin(LatR2) - (Math.sin(LatR1) * D)) / (Math.cos(LatR1) * Math.sin(range));`
where if either (or both!) of cos(LatR1) or sin(range) evaluate to zero, then it's game over.

4. I'm still debugging the translation... of course the NaN comes from division by range = 0

5. I wouldn't post a complex polygon as I know it can't handle that.
I think if the code in the library falls over repeated points then it should just check for these points and skip them.
I added some code to my app to clear repeated points.

Calculating radius from latitude is quite simple and fast, so added that. I take the latitude from the calculated centroid of the area.
With that the code to calculate the area, posted by Klaus seems fine and has the benefit that there is no library, so you can see what is happening.
Sofar all working fine.

Code:
`Sub GetPolygonAreaSimple(lstPoints As List) As Double     Dim i As Int    Dim Lat1, Lat2, Lng1, Lng2 As Double    Dim dArea, LatScale, LngScale As Double    Dim dRadius As Double    Dim gpCentroid As OSMDroid_GeoPoint    'Get the areas.    Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(0) 'dNavA = 6378137 If bCalculateRadius Then 'this looks the better way  gpCentroid = GetAreaCentroid(lstPoints)  dRadius = GetEarthRadius(gpCentroid.Latitude)  LatScale =     cPI * dRadius / 180        ' converts latitude degrees to m Else  LatScale =     cPI * dNavA / 180        ' converts latitude degrees to m End If    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  'General.RunLog("GetPolygonAreaSimple, " & i & ": " & oGeoPoint2.Latitude & "," & oGeoPoint2.Longitude)    Next     Return Round2(Abs(dArea / 1000000), 4) End SubSub GetEarthRadius(dLatitude As Double) As Double  Dim dRadius As Double Dim dRadiusP As Double = 6356752 'at poles Dim dRadiusE As Double = 6378137 'at equator   dRadius = Sqrt((Power(dRadiusE * dRadiusE * Cos(dLatitude), 2) + Power(dRadiusP * dRadiusP * Sin(dLatitude), 2)) / _    (Power(dRadiusE * Cos(dLatitude), 2) + Power(dRadiusP * Sin(dLatitude), 2)))     Return dRadius End SubSub GetAreaCentroid(lstPoints As List) As 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("GetAreaCentroid, polygon_area: " & polygon_area)     X = X / 6 / polygon_area    Y = Y / 6 / polygon_area  General.RunLog("GetAreaCentroid, Y: " & Y) General.RunLog("GetAreaCentroid, X: " & X)    '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 oGeoPointCenterEnd Sub'this is not the actual area in km2Sub GetPolygonArea(lstPoints As List) As 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`

RBS

6. My translation from Java to B4A thus far is below. I did trip over another little anomaly in the Bearing function, which could be of my own doing by caching intermediate Sin/Cos values into temporary variables, probably loses a few low-end bits of precision in the process.

But my main problem is still that I am still getting large numbers out of Abs(Abs(MinSum) - (N - 2) * cPI) on the Spherical() return line.

My other looming problem is that it is nighttime here and I am starting to go a bit fuzzy, so I'm letting go until tomorrow. Code:
`Sub Bearing(LatR1 As Double, LongR1 As Double, LatR2 As Double, LongR2 As Double) As Double     Dim SinLatR1 As Double = Sin(LatR1)    Dim SinLatR2 As Double = Sin(LatR2)    Dim CosLatR1 As Double = Cos(LatR1)    Dim CosLatR2 As Double = Cos(LatR2)     Dim L As Double = LongR1 - LongR2    Dim D As Double = (SinLatR1 * SinLatR2) + (CosLatR1 * CosLatR2 * Cos(L))    Dim Range As Double = ACos(D)    Dim C As Double = (SinLatR2 - (SinLatR1 * D)) / (CosLatR1 * Sin(Range))     'sometimes C falls very slightly outside of the range -1..1    'eg C = -1.0000000000000997 and thus ACos(C) = NAN     If C < -1 Then        C = -1    else if C > 1 Then        C = 1    End If                   Dim ReturnValue As Double     If Sin(L) < 0 Then        ReturnValue = ACos(C)    Else        ReturnValue = 2 * cPI - ACos(C)    End If     Return ReturnValueEnd Sub Sub Spherical(Points(,) As Double , N As Int, MeanLat As Double) As Double    Dim es As Double = 0.0067394967423    'eccentricity    Dim twoPi As Double = 2 * cPI    Dim phis As Double = ATan(Tan(MeanLat) / (1 + es))    Dim Sinphis As Double = Sin(phis)    Dim Rs As Double = 6378.137 / Sqrt(1 + es * Sinphis * Sinphis)    Rs = 6378.137    'close enough for time being     Dim Sum1 As Double = 0    Dim Sum2 As Double = 0     For I = 0 To N - 1        Dim J As Int = (I + 1) Mod N        Dim K As Int = (I + 2) Mod N           Dim F1 As Double = Bearing(Points(j, 0), Points(j, 1), Points(k, 0), Points(k, 1))        Dim F2 As Double = Bearing(Points(j, 0), Points(j, 1), Points(i, 0), Points(i, 1))           Sum1 = (Sum1 + F1 + twoPi - F2) Mod twoPi        Sum2 = (Sum2 + F2 + twoPi - F1) Mod twoPi    Next     Dim MinSum As Double = Min(Sum1, Sum2)    'the minimum is the internal angles sum     Return Abs(Abs(MinSum) - (N - 2) * cPI) * Rs * Rs    'in square kmEnd SubSub SphericalArea2(Poligon(,) As Double, N As Int) As Double    Dim MeanLat As Double = 0    Dim Points(N, 2) As Double    'new Array because the content must change     For I = 0 To N - 1        Points(I, 0) = Poligon(I, 0) * cPI / 180        Points(I, 1) = Poligon(I, 1) * cPI / 180        MeanLat = MeanLat + Points(I, 0)    Next     MeanLat = MeanLat / N     Return Spherical(Points, N, MeanLat) End Sub`

Last edited: Jun 9, 2019
7. <thumbs up> x 5

Yea, stick with the flat area calculations. It'll be nice to get the spherical calculations working, but it's more a labour of love than a degree of precision actually required in the real world. Even NASA suggested that straight lines are often close enough, and that a line had to be over a hundred kilometres long before it was out by a pixel on their display compared to an absolutely correct line.

If it's good enough for NASA, then it's good enough for me, and presumably good enough for your front garden landscape planning or whatever your mission is. :-/

8. Checking the rectangle case in SphericalArea I found that a slight difference is required between two points Lat or Long:
Code:
`Sub test2        Dim p(4,2) As Double    p(0,0) = 50    p(0,1) = 30.000001    p(1,0) = 51    p(1,1) = 30    p(2,0) = 51    p(2,1) = 31    p(3,0) = 50.00001    p(3,1) = 31    Dim lst As List    lst.Initialize    For i = 0 To 3        lst.Add(p(i,0) & "," & p(i,1))        Log(lst.Get(i))    Next    Log(nav.SphericalArea2(p,4,0,1))    Log(nav.SphericalArea(lst))End Sub`
I'll try a simplified bearing calculation and see if it makes a change.

9. emexes, your note in the code drew my attention and I found that it causes all the problem. I added two lines to eliminate these cases and everything works fine.
Here is the complete translated code, as an app in b4j:
Code:
`Sub Process_Globals'    Private fx As JFX'    Private MainForm As Form    Private deg2rad As Double = cPI/180End SubSub AppStart (Form1 As Form, Args() As String)'    MainForm = Form1    'MainForm.RootPane.LoadLayout("Layout1") 'Load the layout file.'    MainForm.Show    test2End SubSub test2    Dim points(4,2) 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    Log("Area = " & SphericalArea2(points,4,0,1))  End Sub  'Return true to allow the default exceptions handler to handle the uncaught exception.Sub Application_Error (Error As Exception, StackTrace As String) As Boolean    Return TrueEnd Sub'  The methods calculates the area on the Earth's spherical enclosed inside a polygon.'  n Is the polygon's size.'  GindxLat And GindxLong are the indices of the parameters Lat And Long in the Array.'  The returned area Is in square Km.Public Sub  SphericalArea2(Poligon(,) As Double,  n As Int, _    GindxLat As Int ,  GindxLong As Int) As Double    Dim  meanlat As Double = 0    Dim  points(n,2) As Double ' new Array because the content must change    For i = 0 To n-1        points(i,0) = Poligon(i,GindxLat)* deg2rad ' allocation by the index        points(i,1) = Poligon(i,GindxLong)* deg2rad        meanlat = meanlat + points(i,0)    Next    meanlat = meanlat/n    Return Spherical(points , n, meanlat)End SubSub  Spherical( points(,) As Double ,n As Int, meanlat As Double) As Double      Dim sum1,sum2,phis As Double     Dim es = 0.0067394967423  As Double  ' eccentricity    Dim twoPi = 2 * cPI As Double    phis = ATan(Tan(meanlat)/(1+es))    Dim Rs As Double = 6378.137 /Power((1 + es * Sin(phis) * Sin(phis)), 0.5) ' km    ' find angles For point 0    Dim f1 As Double = Bearing(points(0,0), points(0,1), points(1,0), points(1,1))    Dim f2 As Double = Bearing(points(0,0),points(0,1),points(n-1,0),points(n-1,1))    sum1 = (f1 + twoPi - f2) Mod twoPi    sum2 = (f2 + twoPi - f1) Mod twoPi    ' find angles For point n-1    f1 = Bearing(points(n-1,0),points(n-1,1),points(0,0),points(0,1))    f2 = Bearing(points(n-1,0),points(n-1,1),points(n-2,0),points(n-2,1))    sum1 = sum1 + (f1 + twoPi - f2) Mod twoPi    sum2 = sum2 + (f2 + twoPi - f1) Mod twoPi    ' find angles For the intermediate points    For i = 1 To n-2        f1 = Bearing(points(i,0),points(i,1),points(i+1,0),points(i+1,1))        f2 = Bearing(points(i,0),points(i,1),points(i-1,0),points(i-1,1))        sum1 = sum1 + (f1 + twoPi - f2) Mod twoPi        sum2 = sum2 + (f2 + twoPi - f1) Mod twoPi    Next    sum1 =  Min(sum1,sum2) ' the minimun Is the internal angles sum    Return Abs((Abs(sum1)-(n-2)*cPI)) * Rs*Rs ' in sqr kmEnd Sub             ' This method uses arguments in Radians And returns Radians ! Sub Bearing( LatR1 As Double, LongR1 As Double, LatR2 As Double, LongR2 As Double) As Double          Dim  range, Bearing1 As Double    Dim L As Double = LongR1 - LongR2    Dim D As Double= (Sin(LatR1) * Sin(LatR2)) + (Cos(LatR1) * Cos(LatR2) *Cos(L ))    D = Max(Min(D,1),-1)   ' Addition    range = ACos(D)    Dim C As Double = (Sin(LatR2) - (Sin(LatR1) * D)) /(Cos(LatR1) * Sin(range))    C = Max(Min(C,1),-1)   ' Addition    If (Sin(L) < 0) Then        Bearing1 = ACos(C)    Else        Bearing1 = (2*cPI - ACos(C))    End If    Return  Bearing1End Sub`
When I'll remember how to update a library in Eclipse I'll update the library...
The area of the small rectangle is much smaller than 1 kmr :
Area = 0.06396440387541677

10. Library updated to 1.2 !

emexes likes this.
11. Thanks, will try it out.
Won't have to deal with duplicate points anymore.

RBS

12. Just tested and seems to work well.
SphericalArea2 fine with small areas, eg 0.015 km2.
Did you look at the accuracy of the 3 mentioned algo's?

RBS

13. Please reload the library, I made some more corrections of the same problem in other methods.

14. Will do.

RBS

15. I had a similar discovery myself (a good sleep works wonders), in my translation of your code to:
Code:
`Sum1 = (Sum1 + F1 + twoPi - F2) Mod twoPiSum2 = (Sum2 + F2 + twoPi - F1) Mod twoPi`
the opening bracket needs to move along one step, ie:
Code:
`Sum1 = Sum1 + (F1 + twoPi - F2) Mod twoPiSum2 = Sum2 + (F2 + twoPi - F1) Mod twoPi`
and now I am getting results that very nearly match:

12380.233251685318 from B4A translation

which I'm almost happy with. Will try running same code with B4J (ie x86 FPU, JVM) instead of on Android device (ie ARM FPU, Android JVM) - I doubt that will fix it, but stranger things have happened. Also, I changed the order in which the polygon angles are summed, so I'll try reverting to your way (first angle, last angle, then the middle angles).

What could possibly go wrong? 16. 12380.233251685318 from B4J translation, so that's kinda reassuring.

17. Changng order of polygon angles summing doesn't change anything, which is very reassuring ;-)

12380.233251685318 from B4A translation using original summing order

18. Righto, final result:

12380.233251685318 from my translation to B4X (sums polygon angles in the order that they are presented by the array)
12380.233251685318 from your translation to B4X (sums polygon angles in first, last, remainder order)

I'm prepared to say that the library result difference sure looks like a single-vs-double floating point issue. If pressed, my backup excuse is that the output can only be as good as the input, so we would expect accuracy in the order of the 7 digits of the supplied earth radius precision.

Last edited: Jun 10, 2019
19. Postscript: yeah, I should just let this go, instead of chasing down a less-than-a-millionth-of-a-percent discrepancy, but... :-?

compiled and ran original Java code, still got very slightly different answer to that from Navigation library:

12380.233251685318 from my translation to B4X (sums polygon angles in the order that they are presented by the array)
12380.233251685318 from your translation to B4X (sums polygon angles in first, last, remainder order)
12380.233251685318 from compiling and running original Java code on a Windows computer

so... I don't know. The usual suspects would be:

1/ a floating point literal in the library code is being default cast to single precision rather than double precision
2/ the Math.pow(x. 0.5) comprises two precision-losing operations eg exp(log(x) * 0.5) rather than one square-root operation
3/ deg2rad is single, not double (although I doubt that, given that Math.PI is a double, and the compiler promptly points out the loss of precision)

20. BTW I tried this, and my problem test cases were magically fixed too. derez likes this.