C# GMap.Net 计算多边形面积

float GetArea(List<Vector2> points)
{
    float area2 = 0;
    for (int numPoint = 0; numPoint < points.Count - 1; numPoint++)
    {
        MyPoint point = points[numPoint];
        MyPoint nextPoint = points[numPoint + 1];
        area2 += point.x * nextPoint.y - point.y * nextPoint.x;
    }
    return area2 / 2f;
}

float GetArea(List<Vector3> points, Vector3 normal)
{
    Vector3 vector = Vector3.Zero;
    for (int numPoint = 0; numPoint < points.Count - 1; numPoint++)
    {
        MyPoint point = points[numPoint];
        MyPoint nextPoint = points[numPoint + 1];
        vector += Vector3.CrossProduct(point, nextPoint);
    }
    return (1f / 2f) * Math.Abs(Vector3.DotProduct(vector, normal));
}

Vector3 GetCentroid(List<Vector3> points)
{
    Vector3 vector = Vector3.Zero;
    Vector3 normal = Vector3.CrossProduct(points[0], points[1]);  // Gets the normal of the first triangle (it is used to know if the contribution of the triangle is positive or negative)
    normal = (1f / normal.Length) * normal;  // Makes the vector unitary
    float sumProjectedAreas = 0;
    for (int numPoint = 0; numPoint < points.Count - 1; numPoint++)
    {
        MyPoint point = points[numPoint];
        MyPoint nextPoint = points[numPoint + 1];
        float triangleProjectedArea = Vector3.DotProduct(Vector3.CrossProduct(point, nextPoint), normal);
        sumProjectedAreas += triangleProjectedArea;
        vector += triangleProjectedArea  * (point + nextPoint);
    }
    return (1f / (6f * sumProjectedAreas)) * vector;
}

Vector3 GeographicCoordinatesToPoint(float latitude, float longitude)
{
    return EarthRadius * new Vector3(Math.Cos(latitude) * Math.Cos(longitude), Math.Cos(latitude) * Math.Sin(longitude), Math.Sin(latitude));
}

// Converts the latitude/longitude coordinates to 3D coordinates
List<Vector3> pointsIn3D = (from point in points
                            select GeographicCoordinatesToPoint(point.Latitude, point.Longitude))
                           .ToList();

// Gets the centroid (to have the normal of the vector space)
Vector3 centroid = GetCentroid(pointsIn3D );

// As we are on a sphere, the normal at a given point is the colinear to the vector going from the center of the sphere to the point.
Vector3 normal = (1f / centroid.Length) * centroid;  // We want a unitary normal.

// Finally the area is computed using:
float area = GetArea(pointsIn3D, normal);

public struct Vector3
{
    public static readonly Vector3 Zero = new Vector3(0, 0, 0);

    public readonly float X;
    public readonly float Y;
    public readonly float Z;

    public float Length { return Math.Sqrt(X * X + Y * Y + Z * Z); }

    public Vector3(float x, float y, float z)
    {
         X = x;
         Y = y;
         Z = z;
    }

    public static Vector3 operator +(Vector3 vector1, Vector3 vector2)
    {
        return new Vector3(vector1.X + vector2.X, vector1.Y + vector2.Y, vector1.Z + vector2.Z);
    }
    public static Vector3 operator -(Vector3 vector1, Vector3 vector2)
    {
        return new Vector3(vector1.X - vector2.X, vector1.Y - vector2.Y, vector1.Z - vector2.Z);
    }
    public static Vector3 operator *(float scalar, Vector3 vector)
    {
        return new Vector3(scalar * vector.X, scalar * vector.Y, scalar * vector.Z);
    }

    public static float DotProduct(Vector3 vector1, Vector3 vector2)
    {
        return vector1.X * vector2.X + vector1.Y * vector2.Y + vector1.Z * vector2.Z;
    }
    public static Vector3 CrossProduct(Vector3 vector1, Vector3 vector2)
    {
        return return new Vector3(vector1.Y * vector2.Z - vector1.Z * vector2.Y,
                                  vector1.Z * vector2.X - vector1.X * vector2.Z,
                                  vector1.X * vector2.Y - vector1.Y * vector2.X);
    }
}

List<PointLatLng> firstTwoPoints = new List<PointLatLng>();
         firstTwoPoints.Add(poly.Points[0]);
         firstTwoPoints.Add(poly.Points[1]);

         GMapPolygon oneLine = new GMapPolygon(firstTwoPoints,"testesddfsdsd"); //Create new polygone from messuring the distance.
         double lengteLocalLine =
            Math.Sqrt(((poly.LocalPoints[1].X - poly.LocalPoints[0].X)*(poly.LocalPoints[1].X - poly.LocalPoints[0].X)) +
                      ((poly.LocalPoints[1].Y - poly.LocalPoints[0].Y)*(poly.LocalPoints[1].Y - poly.LocalPoints[0].Y))); //This calculates the length of the line in LocalPoints. 
         double pointInKm = oneLine.Distance / lengteLocalLine; //This gives me the length of 1 LocalPoint in km.

         List<Carthesian> waarden = new List<Carthesian>(); 

//Here we fill the list "waarden" with the points.
//NOTE: the last value is NOT copied because this is handled in calculation method.
         foreach (GPoint localPoint in poly.LocalPoints)
         {
            waarden.Add(new Carthesian(){X = (localPoint.X * pointInKm), Y = (localPoint.Y * pointInKm)});
         }

         MessageBox.Show("" + GetArea(waarden)*1000000);

    }

//Method for calculating area
      private double GetArea(IList<Carthesian> points)
      {
         if (points.Count < 3)
         {
            return 0;
         }
         double area = GetDeterminant(points[points.Count - 1].X , points[points.Count - 1].Y, points[0].X, points[0].Y);

         for (int i = 1; i < points.Count; i++)
         {
            //Debug.WriteLine("Lng: " + points[i].Lng + "   Lat:" + points[i].Lat);
            area += GetDeterminant(points[i - 1].X, points[i - 1].Y , points[i].X,  points[i].Y);
         }
         return Math.Abs(area / 2);
      }

//Methode for getting the Determinant

      private double GetDeterminant(double x1, double y1, double x2, double y2)
      {
         return x1 * y2 - x2 * y1;
      }



//This class is just to make it nicer to show in code and it also was from previous tries
    class Carthesian
       {
          public double X { get; set; }
          public double Y { get; set; }
          public double Z { get; set; }
       }
Because we calculate the surface using 2D there