我正在寻找一种计算多边形面积的方法。
我想要做的事情是，让使用我的程序的用户可以创建一个多边形来标记他的房产。现在我想知道这个区域的面积大小，以便告诉用户他的房产有多大。
单位为平方米或平方千米或公顷。
多边形的点具有纬度和经度。
我正在使用带有WPF和GMap.NET的C#。地图位于一个WindowsFormHost中，因此我可以使用来自GMap.Net的Winforms工具，因为它提供了叠加等功能。
我希望有人能帮助我或向我展示解释了这个问题的文章。
C# GMap.Net 计算多边形面积
3
- Boeykes
3个回答
5
使用二维向量空间逼近（局部切线空间）
在本节中，我将详细介绍如何得出这些公式。
让我们把多边形的点称为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 操作。
在您的情况下，您有一个纬度和经度。然后，您不处于一个二维欧几里得空间中。这是一个球形空间，计算任何多边形的面积要复杂得多。 但是，它在局部上是同胚于一个二维向量空间（使用切空间）。 因此，如果您尝试测量的面积不太大（几公里），上述公式应该适用。
现在，您只需要找到多边形的法线。为了这样做，并减少误差（因为我们在估计面积），我们使用多边形重心处的法线。重心由以下公式给出：
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);
}
}
- Cédric Bignon
3
谢谢您的快速回复，我认为这可以帮助我很多。我会尝试今天实现它。再次感谢您的快速回复。 - Boeykes
请问您能给我一个 Vector3 类的例子吗？因为我找不到这个类有哪些属性以及如何在叉积计算中使用“normal”。 - Boeykes
我已经添加了Vector3的代码。我仍在考虑使用球面三角学来解决这个问题的确切方法。 - Cédric Bignon
1
我使用了来自Cédric的一部分代码和网络上的代码来修复它。
我通过使用poly.Points和poly.LocalPoints来修复它。 poly.Points是纬度和经度，而LocalPoints是屏幕上地图中心点的点。
C#库有一个函数可以计算距离（公里），所以我计算了距离，然后计算了LocalPoints中的距离。将LocalPoints除以公里数长度，然后你就知道1个Localpoint在公里中有多长。
代码：
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; }
}
因为我们是使用 2D 计算表面，所以会有一点误差，但对于我的应用来说这是可以接受的。
感谢 Cédric 回答我的问题并帮助我解决了我遇到的问题。
- Boeykes
0
使用后端数据库，如SQL Server 2008或MySql，将多边形的点发送到服务器并返回面积、长度、周长、交集等信息会更容易。
如果可行，可以搜索Sql Server地理/几何数据类型中的STArea()或STIntersect。
这是我一直在研究的一个示例。
using Microsoft.SqlServer.Types;
using System.Data.SqlClient;
GMap.NET.WindowsForms.GMapOverlay o = new GMapOverlay();
GMap.NET.WindowsForms.GMapOverlay o1 = new GMapOverlay();
List<PointLatLng> p = new List<PointLatLng>();
List<string> p1 = new List<string>();
//below adds a marker to the map upon each users click. At the same time adding that Lat/Long to a <PointLatLng> list
private void gMapControl1_MouseClick(object sender, MouseEventArgs e)
{
if (e.Button == MouseButtons.Right )
{
p.Add(new PointLatLng(Convert.ToDouble(gMapControl2.FromLocalToLatLng(e.X, e.Y).Lat), Convert.ToDouble(gMapControl2.FromLocalToLatLng(e.X, e.Y).Lng)));
p1.Add(gMapControl2.FromLocalToLatLng(e.X, e.Y).Lng + " " + gMapControl2.FromLocalToLatLng(e.X, e.Y).Lat);
GMarkerGoogle marker = new GMarkerGoogle(gMapControl2.FromLocalToLatLng(e.X, e.Y), GMarkerGoogleType.red_small);
// marker.Tag =(gMapControl1.FromLocalToLatLng(e.X, e.Y).Lng + " " + gMapControl1.FromLocalToLatLng(e.X, e.Y).Lat);
o.Markers.Add(marker);
gMapControl2.Overlays.Add(o);
}
}
//Then the below code will add that <PointLatLng> List to a SQL query and return Area and Centoid of polygon. Area is returned in Acres
private void gMapControl1_MouseDoubleClick(object sender, MouseEventArgs e)
{
if (e.Button == MouseButtons.Left)
{
try
{
o.Clear();
n = new GMapPolygon(p, "polygon");
n.Fill = new SolidBrush(Color.Transparent);
n.Stroke = new Pen(Color.Red, 1);
o.Polygons.Add(n);
gMapControl2.Overlays.Add(o);
StringBuilder a = new StringBuilder();
StringBuilder b = new StringBuilder();
p1.ToArray();
for (int i = 0; i != p1.Count; i++)
{
a.Append(p1[i].ToString() + ",");
}
a.Append(p1[0].ToString());
cs.Open();
SqlCommand cmd = new SqlCommand("Declare @g geography; set @g = 'Polygon((" + a + "))'; Select Round((@g.STArea() *0.00024711),3) As Area", cs);
SqlCommand cmd1 = new SqlCommand("Declare @c geometry; set @c =geometry::STGeomFromText('Polygon((" + a + "))',0); Select Replace(Replace(@c.STCentroid().ToString(),'POINT (',''),')','')AS Center", cs);
poly = "Polygon((" + a + "))";
SqlDataReader sdr = cmd.ExecuteReader();
while (sdr.Read())
{
txtArea.Text = sdr["Area"].ToString();
}
sdr.Close();
SqlDataReader sdr1 = cmd1.ExecuteReader();
while (sdr1.Read())
{
center = sdr1["Center"].ToString();
lat = center.Substring(center.IndexOf(" ") + 1);
lat = lat.Remove(9);
lon = center.Substring(0, (center.IndexOf(" ")));
lon = lon.Remove(10);
txtCenter.Text = lat + ", " + lon;
}
sdr1.Close();
}
catch (Exception ex)
{
MessageBox.Show("Please start the polygon over, you must create polygon in a counter-clockwise fasion","Counter-Clockwise Only!",MessageBoxButtons.OK,MessageBoxIcon.Error);
}
finally
{
p.Clear();
p1.Clear();
cs.Close();
o.Markers.Clear();
}
}
}
- Nate S.
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接