Announcement

Collapse
No announcement yet.

A method to compute the area of a spherical polygon

Collapse
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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, 02:34 PM.

  • #2
    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, 04:30 PM.
    My World Wind Java Blog & WW.net Plugins page

    Comment


    • #3
      Spherical Area

      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.

      Comment


      • #4
        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

        Comment


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

          Regards,
          Harika M

          [long quote removed]

          Comment


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

            Comment


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

              Comment


              • #8
                *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.

                Comment


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

                  Comment


                  • #10
                    numerical failures of L'Huiler's Theorem

                    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.

                    Comment


                    • #11
                      Are you doing your calculations using a sphere or are you compensating for the ellipsoidal shape of the earth?
                      Neil
                      http://www.nlneilson.com

                      Comment


                      • #12
                        On a sphere ...

                        Originally posted by nlneilson View Post
                        Are you doing your calculations using a sphere or are you compensating for the ellipsoidal shape of the earth?
                        These are for a sphere.

                        The thing about L'Huiller's that started us on this path is that if you take a complicated shape like the tree shape described, on a small scale, and move that fixed shape about a spherical Earth, recalculating the area, L'Huiller's will give you a different answer. For example, if it's placed up at about 45N, it does okay. Move it to the Equator, and it gives a different answer.

                        We of course very carefully checked that our rotation of the points was correct. The rotation is done by first calculating unit vector positions for all the points in the tree shape, applying the same rotation matrix to them (done in R using the standard long precision), carefully renormalizing before and after. Moreover, the other algorithms, like a unit vector area algorithm we have, and the NASA JPL one, do give the same area for the tree at the Equator as they give up at 45N, and they agree with each other.

                        Alas, these algorithms also fail on the clever test case which the people at Mathworks cooked up for their areaint function in MATLAB. So, workin' on this ....

                        If there's interest, I can take on publishing a code for doing this here, in ANSI C, but I have a commitment to first do this and deliver it for within B. This is almost entirely on my own time, so I can't promise a schedule.

                        The Green's theorem technique I'm looking at would generalize to WGS-84, although that's not what I'm aiming for right now, and whatever code I might offer will only work for a sphere. Also, the radius of said sphere will be a parameter (since there are several ones used for spherical approximations of Earth), although it will default to something sensible.

                        Thanks for your question.

                        Comment


                        • #13
                          I am mostly concerned with geodesic distance and angles.
                          Have been using the Vincenty Formula/s using C++.
                          Was recently informed of the accuracy limitations by a mathematician who has done formulas at greater accuracy, more than I need.

                          I like to tinker with numbers.
                          Neil
                          http://www.nlneilson.com

                          Comment


                          • #14



                            I have made a c++ program to make icon and line/polygon feature layer.xml files. I have been using the haversine formula to calculate the distances around a point as above, distance between two points, bearing between two points, position along a route, etc. It is quite useful! I also used trig functions for the altitude and radius.
                            Last edited by danzig70; 02-13-2012, 03:08 AM.

                            Comment


                            • #15
                              double lam1 = 0, lam2 = 0, beta1 =0, beta2 = 0, cosB1 =0, cosB2 = 0;
                              use double type,There is an error;but user float,it is right
                              float lam1 = 0, lam2 = 0, beta1 =0, beta2 = 0, cosB1 =0, cosB2 = 0;

                              Comment

                              Working...
                              X