# Android Questionarea of polygon

#### RB Smissaert

##### Well-Known Member
Longtime 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.

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

Sub FindCentroid(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("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

#### derez

##### Expert
Longtime 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] ; // new array because the content must change
for (int i = 0;i<n;i++)
{
points[i] = Poligon[i][GindxLat]* deg2rad ; // allocation by the index
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
}

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

##### Expert
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

##### Expert
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
Longtime User
I'm still debugging the translation... of course the NaN comes from division by range = 0

#### RB Smissaert

##### Well-Known Member
Longtime 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 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)
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 dRadiusP As Double = 6356752 'at poles
Dim dRadiusE As Double = 6378137 'at equator

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

##### Expert
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

##### Expert
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. :-/

#### derez

##### Expert
Longtime 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
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
Longtime 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.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
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

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

range = ACos(D)
Dim C As Double = (Sin(LatR2) - (Sin(LatR1) * D)) /(Cos(LatR1) * Sin(range))
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
Longtime User
Library updated to 1.2 !

• emexes

#### RB Smissaert

##### Well-Known Member
Longtime User
Library updated to 1.2 !

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

RBS

#### RB Smissaert

##### Well-Known Member
Longtime User
Library updated to 1.2 !

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

Longtime User

Longtime User

Will do.

RBS

#### emexes

##### Expert
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.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

##### Expert
12380.233251685318 from B4J translation, so that's kinda reassuring.

#### emexes

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

12380.233251685318 from B4A translation using original summing order

#### emexes

##### Expert
Checking...
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:

#### emexes

##### Expert
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.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)

#### emexes

##### Expert
Library updated to 1.2 !
BTW I tried this, and my problem test cases were magically fixed too. • derez

Android Question Draw a Polygon on a Map
Replies
4
Views
1K
Replies
25
Views
10K
Replies
2
Views
1K