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

Question

我正在寻找一种计算多边形下表面的方法。

我想要完成的事情是，使用我的程序的用户可以创建一个多边形来标记他的属性。现在我想知道表面积是多少，这样我就可以告诉用户他的属性(property)有多大。

单位 m²、km² 或公顷。

多边形的点具有纬度和经度。

我正在将 C# 与 WPF 和 GMap.NET 结合使用。 map 位于 WindowsFormHost 中，因此我可以使用 GMap.Net 中的 Winforms 东西，因为这可以避免覆盖等。

我希望有人可以帮助我，或者给我看一篇帖子，其中解释了我没有找到的内容。

Answer 1

使用二维向量空间近似(局部切线空间)

在本节中，我可以详细说明我是如何得出这些公式的。

让我们记下Points多边形的点(其中Points[0] == Points[Points.Count - 1]以闭合多边形)。

接下来的方法背后的想法是将多边形分割成三角形(面积是所有三角形面积的总和)。但是，为了通过简单分解支持所有多边形类型(不仅仅是星形多边形)，一些三角形的贡献是负的(我们有一个“负”面积)。我使用的三角形分解是:{(O, Points[i], Points[i + 1]}，其中O是仿射空间的原点。

非自相交多边形(在欧几里德几何中)的面积由下式给出:

二维:

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;
}

在 3D 中，给定法线，多边形的单一法线(平面):

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));
}

在前面的代码中，我假设您有一个带有 Add、Subtract、Multiply、CrossProduct 的 Vector3 结构和 DotProduct 操作。

就您而言，您有纬度和经度。那么，你就不再处于二维欧几里得空间中。它是一个球形空间，计算任何多边形的面积要复杂得多。然而，它与二维向量空间局部同胚(使用切线空间)。然后，如果您尝试测量的区域不太宽(几公里)，则上述公式应该有效。

现在，您只需找到多边形的法线即可。为此，并减少误差(因为我们正在近似面积)，我们在多边形的质心处使用法线。质心由下式给出:

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 添加了一个新属性:Vector3.Length

最后，将纬度和经度转换为 Vector3:

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);

Vector3 结构

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);
    }
}