Android Question area of polygon

derez

Expert
Licensed User
Thanks for testing, emexes ! The java source for the Surfacearea2 method is following, if you find anything wrong just point it out:
B4X:
     /**
      * 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][2] ; // new array because the content must change
            for (int i = 0;i<n;i++)
            {
                points[i][0] = Poligon[i][GindxLat]* deg2rad ; // allocation by the index
                points[i][1] = Poligon[i][GindxLong]* deg2rad ;
                meanlat += points[i][0];
            }
            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[0][0],points[0][1],points[1][0],points[1][1]);
            double f2 = Bearing(points[0][0],points[0][1],points[n-1][0],points[n-1][1]);
            sum1 = (f1 + twoPi - f2)%twoPi;
            sum2 = (f2 + twoPi - f1)%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 += (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][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 += (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:

emexes

Well-Known Member
Licensed User
The code can be easily translated to enable b4a debug.
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:
B4X:
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 sum
5            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.

if you find anything wrong just point it out
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.
 

emexes

Well-Known Member
Licensed User
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.
which is probably this line:
B4X:
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.
 

derez

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

RB Smissaert

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!View attachment 81151
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.

B4X:
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 Sub

Sub 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 Sub

Sub 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 oGeoPointCenter

End Sub

'this is not the actual area in km2
Sub 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
 

emexes

Well-Known Member
Licensed User
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.

:)

B4X:
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 ReturnValue

End 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 km

End Sub

Sub 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:

emexes

Well-Known Member
Licensed User
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.
So far all working fine.
<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.

:)

Except for O-rings, let's not follow NASA's lead on those...

:-/
 

derez

Expert
Licensed User
Checking the rectangle case in SphericalArea I found that a slight difference is required between two points Lat or Long:
B4X:
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.
 

derez

Expert
Licensed User
'sometimes C falls very slightly outside of the range -1..1
'eg C = -1.0000000000000997 and thus ACos(C) = NAN
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:
B4X:
Sub Process_Globals
'    Private fx As JFX
'    Private MainForm As Form
    Private deg2rad As Double = cPI/180

End Sub

Sub AppStart (Form1 As Form, Args() As String)
'    MainForm = Form1
    'MainForm.RootPane.LoadLayout("Layout1") 'Load the layout file.
'    MainForm.Show
    test2
End Sub

Sub 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 True
End 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 Sub


Sub  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 km
End 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  Bearing1
End 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
 

derez

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

emexes

Well-Known Member
Licensed User
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.
I had a similar discovery myself (a good sleep works wonders), in my translation of your code to:
B4X:
Sum1 = (Sum1 + F1 + twoPi - F2) Mod twoPi
Sum2 = (Sum2 + F2 + twoPi - F1) Mod twoPi
the opening bracket needs to move along one step, ie:
B4X:
Sum1 = Sum1 + (F1 + twoPi - F2) Mod twoPi
Sum2 = Sum2 + (F2 + twoPi - F1) Mod twoPi
and now I am getting results that very nearly match:

12380.233223394947 from Navigation library
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?

:)
 

emexes

Well-Known Member
Licensed User
Changng order of polygon angles summing doesn't change anything, which is very reassuring ;-)

12380.233251685318 from B4A translation using original summing order
 

emexes

Well-Known Member
Licensed User
Checking...
Righto, final result:

12380.233223394947 from Navigation library 1.1
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:

emexes

Well-Known Member
Licensed User
Checking...
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.233223394947 from Navigation library 1.1
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)
 
Top