Calculating area of a polygon in geographic coordinates

[OrenGal]

Spent a lot of time finding the best function to calculate the area of a polygon in web mercator projection. Because the Silverlight API converts from wkid 102100 to wkid 4326 in the client side, It was easy to convert the point collection and then calculate the area.

Copy this function if you don't have it. It may save you time.

        public double get_area_of_polygon(PointCollection pc_all_4326)
        {
            double lam1 = 0, lam2 = 0, beta1 = 0, beta2 = 0, cosB1 = 0, cosB2 = 0;
            double hav = 0;
            double sum = 0;

            for (int j = 0; j < pc_all_4326.Count; j++)
            {
                int k = j + 1;
                if (j == 0)
                {
                    lam1 = pc_all_4326.X;
                    beta1 = pc_all_4326.Y;
                    lam2 = pc_all_4326[j + 1].X;
                    beta2 = pc_all_4326[j + 1].Y;
                    cosB1 = Math.Cos(beta1);
                    cosB2 = Math.Cos(beta2);
                }
                else
                {
                    k = (j + 1) % pc_all_4326.Count;
                    lam1 = lam2;
                    beta1 = beta2;
                    lam2 = pc_all_4326.X;
                    beta2 = pc_all_4326.Y;
                    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) * m_earth_radius * m_earth_radius;
        }

public static double Haversine(double x)
{
    return (1.0 - Math.Cos(x)) / 2.0;
}

A live demo of this measure can be seen here http://bit.ly/dh1Qfu
Measured some squres near the poles and it seems right.
Oren.