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

Question

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

我想要完成的是,使用我的程序的用户可以创建一个多边形来标记他的属性.现在我想知道表面区域是什么,所以我可以告诉用户他的财产有多大.

单位为m²或km²或公顷.

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

我在WPF和GMap.NET上使用C#.地图是WindowsFormHost这样的,所以我可以使用GMap.Net的Winforms东西,因为这可以提供覆盖等.

我希望有人可以帮助我或给我一个帖子,其中解释我没有找到.

Answer 1

使用2D矢量空间近似(局部切线空间)

在本节中,我将详细介绍如何使用这些公式.

让我们注意Points多边形的点(Points[0] == Points[Points.Count - 1]关闭多边形的位置).

下一种方法背后的想法是将多边形分成三角形(该区域是所有三角形区域的总和).但是,为了通过简单的分解(不仅是星形多边形)来支持所有多边形类型,一些三角形贡献是负的(我们有一个"负"区域).我使用的三角形分解是:{(O, Points[i], Points[i + 1]}其中O是仿射空间的起源.

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

在2D中:

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中,给定normal,多边形的酉法线(平面):

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结构有Add,Subtract,Multiply,CrossProduct和DotProduct操作.

在你的情况下,你有一个经度和经度.然后,您不在2D欧几里德空间.它是一个球形空间,计算任何多边形的面积要复杂得多.但是,它与2D矢量空间局部同胚(使用切线空间).然后,如果您尝试测量的区域不是太宽(几公里),则上述公式应该有效.

现在,您只需要找到多边形的法线.为了做到这一点,并减少误差(因为我们接近该区域),我们使用多边形质心处的法线.质心由下式给出:

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