A method to compute the area of a spherical polygon

[Jason_Steven]

Given a spherical polygon described by its vertices,find its area.

/// <summary>
/// Haversine function : hav(x) = (1-cos(x))/2
/// </summary>
/// <param name="x"></param>
/// <returns>Returns the value of Haversine function</returns>
public double Haversine( double x )
{
        return ( 1.0 - Math.Cos( x ) ) / 2.0;
}

/// <summary>
/// Compute the Area of a Spherical Polygon
/// </summary>
/// <param name="lat">the latitudes of all vertices(in radian)</param>
/// <param name="lon">the longitudes of all vertices(in radian)</param>
/// <param name="r">spherical radius</param>
/// <returns>Returns the area of a spherical polygon</returns>
public double SphericalPolygonArea( double[ ] lat , double[ ] lon , double r )
{
       double lam1 = 0, lam2 = 0, beta1 =0, beta2 = 0, cosB1 =0, cosB2 = 0;
       double hav = 0;
       double sum = 0;

       for( int j = 0 ; j < lat.Length ; j++ )
       {
	int k = j + 1;
	if( j == 0 )
	{
	     lam1 = lon[j];
	     beta1 = lat[j];
	     lam2 = lon[j + 1];
	     beta2 = lat[j + 1];
	     cosB1 = Math.Cos( beta1 );
	     cosB2 = Math.Cos( beta2 );
             }
	else
	{
	     k = ( j + 1 ) % lat.Length;
	     lam1 = lam2;
	     beta1 = beta2;
	     lam2 = lon[k];
                  beta2 = lat[k];
	     cosB1 = cosB2;
	     cosB2 = Math.Cos( beta2 );
	}
	if( lam1 != lam2 )
	{
	     hav = Haversine( beta2 - beta1 ) + 
                          cosB1 * cosB2 * Haversine( lam2 - lam1 );
	     double a = 2 * Math.Asin( Math.Sqrt( hav ) );
	     double b = Math.PI / 2 - beta2;
	     double c = Math.PI / 2 - beta1;
	     double s = 0.5 * ( a + b + c );
	     double t = Math.Tan( s / 2 ) * Math.Tan( ( s - a ) / 2 ) *  
                                Math.Tan( ( s - b ) / 2 ) * Math.Tan( ( s - c ) / 2 );

	     double excess = Math.Abs( 4 * Math.Atan( Math.Sqrt( 
                                        Math.Abs( t ) ) ) );

	     if( lam2 < lam1 )
	     {
		excess = -excess;
	     }

	     sum += excess;
	}
      }
      return Math.Abs( sum ) * r * r;
}

The lat and lon parameters is a one-to-one relation .
For example
(lat[0],lon[0]) is a coordinate pair,and so on.

This theory for this algorithm is as follows:
http://mathworld.wolfram.com/SphericalPolygon.html

[patmurris]

Nice, thanks. I guess most of the code is dedicated to computing the azimuth angle from one vertice to the other so as to determine the spherical excess. Isn't there already such a method in MathEngine?

Edit: did you check with polygons spanning the date line or going over the poles?

[Jason_Steven]

Originally Posted by patmurris

Nice, thanks. I guess most of the code is dedicated to computing the azimuth angle from one vertice to the other so as to determine the spherical excess. Isn't there already such a method in MathEngine?

Edit: did you check with polygons spanning the date line or going over the poles?

There doesn't seem to be such a method in MathEngine.
The original algorithm was written in C from Internet.
It seemed OK when I tested it.

[patmurris]

Originally Posted by Jason_Steven

There doesn't seem to be such a method in MathEngine.

How about

/// <summary>
/// Calculates the azimuth from latA/lonA to latB/lonB
/// Borrowed from http://williams.best.vwh.net/avform.htm
/// </summary>
public static Angle Azimuth( Angle latA, Angle lonA, Angle latB, Angle lonB )
{
    double cosLatB = Math.Cos(latB.Radians);
    Angle tcA = Angle.FromRadians( Math.Atan2(
        Math.Sin(lonA.Radians - lonB.Radians) * cosLatB,
        Math.Cos(latA.Radians) * Math.Sin(latB.Radians) - 
        Math.Sin(latA.Radians) * cosLatB * 
        Math.Cos(lonA.Radians - lonB.Radians)));
    if(tcA.Radians < 0) 
        tcA.Radians = tcA.Radians + Math.PI*2;
    tcA.Radians = Math.PI*2 - tcA.Radians;
        return tcA;
}

Thanks!!This code was very useful for me in calculating the spherical polygon area for my project.

Regards,
Harika M

[long quote removed]

[Robert99]

Can I use the code above under an open-source (GPL) license? Thanks, Robert

[zkgis]

Hello I want to calculate the Surface area ,I mean I want to think over the hypsography of mountain
is anybody could share his code
thank you very mach !
I'm Waiting

[QCoder]

*sigh* While mathematically correct, the code is unusable in numerical work. Check Some Algorithms for Polygons on a Sphere to see why, and a better algorithm. I wouldn't feel too bad if you did use it, however, since apparently even the developers of the R package geosphere grabbed and used it, at least through the version available as of 10th June 2011. Also, it's probably not terribly far off. I'm hoping to update the R package in cooperation with the author.

what do you mean by "the code is unusable in numerical work"?

[QCoder]

Originally Posted by Edward

what do you mean by "the code is unusable in numerical work"?

My judgment is based upon failures of the code in specific tests. Naturally, each application is different, and the direct use of the theorem and the code may be appropriate for some applications. In the case of a project I'm involved with, it was not, because we needed accuracies at both extremely large and extremely small spatial scales. I tracked the code back to this source.

To see some examples where the code performs fair-to-poor, consider this case:



then this case:



and finally the failure of this case, where the L'Huiler's implementation is given by the value associated with the tag "geosphere":



The latter is the testcase advertised at the MATLAB site I previously referenced. The MATLAB algorithm calculates a line integral numerically, based upon Green's Theorem adapted to the purpose of finding the area of a polygonal patch on a sphere. Note they use a unit radius for testing purposes. (It's easier to see differences that way, instead of seeing just big numbers.) It's hard to find good algorithms to solve this clever case.

I am in the process of deriving a similar algorithm based upon Green's Theorem for R. The complicated "Christmas tree-like shape" used in some of this testing is available here.