Page 1 of 2 12 LastLast
Results 1 to 10 of 15

Thread: A method to compute the area of a spherical polygon

  1. #1
    Junior Member
    Join Date
    Nov 2008
    Location
    United Kingdom
    Posts
    29

    Smile A method to compute the area of a spherical polygon

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

    Code:
    
    /// <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
    Last edited by Jason_Steven; 11-24-2008 at 01:34 PM.

  2. #2
    WWJ Consultant patmurris's Avatar
    Join Date
    Jun 2005
    Location
    Saint-Paul de Vence, Alpes Maritimes, France
    Posts
    3,383

    Default

    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?
    Last edited by patmurris; 11-24-2008 at 03:30 PM.
    My World Wind Java Blog & WW.net Plugins page

  3. #3
    Junior Member
    Join Date
    Nov 2008
    Location
    United Kingdom
    Posts
    29

    Smile Spherical Area

    Quote Originally Posted by patmurris View Post
    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.

  4. #4
    WWJ Consultant patmurris's Avatar
    Join Date
    Jun 2005
    Location
    Saint-Paul de Vence, Alpes Maritimes, France
    Posts
    3,383

    Default

    Quote Originally Posted by Jason_Steven View Post
    There doesn't seem to be such a method in MathEngine.
    How about
    Code:
    /// <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;
    }
    My World Wind Java Blog & WW.net Plugins page

  5. #5
    Unregistered
    Guest

    Smile

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

    Regards,
    Harika M

    [long quote removed]

  6. #6
    Junior Member
    Join Date
    Apr 2010
    Posts
    1

    Default

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

  7. #7
    Junior Member
    Join Date
    Oct 2010
    Posts
    1

    Default

    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

  8. #8
    Junior Member
    Join Date
    Jul 2011
    Location
    Massachusetts
    Posts
    3

    Default

    *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.

  9. #9
    Edward
    Guest

    Default

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

  10. #10
    Junior Member
    Join Date
    Jul 2011
    Location
    Massachusetts
    Posts
    3

    Smile numerical failures of L'Huiler's Theorem

    Quote Originally Posted by Edward View Post
    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.

Thread Information

Users Browsing this Thread

There are currently 1 users browsing this thread. (0 members and 1 guests)

Similar Threads

  1. WWJ SDK Alpha 4.1 - 0.4.1 available
    By tag in forum WWJ Release Announcements
    Replies: 56
    Last Post: 11-14-2009, 12:51 PM
  2. How to calculate area of an irregular polygon by 3D coordinates?
    By Jason_Steven in forum Developers' Corner
    Replies: 9
    Last Post: 11-13-2008, 07:11 AM
  3. add Urban Area mode as an icon on the toolbar!
    By Phil T in forum Add-ons & Scripts
    Replies: 10
    Last Post: 02-25-2005, 04:48 PM

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •