我正在寻找一种方法来计算多边形下的表面。
我想要完成的事情是，使用我的程序的用户可以创建一个多边形来标记他的属性。现在我想知道表面积是多少，这样我就可以告诉用户他的财产有多大。
单位平方米或平方公里或公顷。
多边形的点有纬度和经度。
我将 C# 与 WPF 和 GMap.NET 一起使用。该地图位于WindowsFormHost 中，因此我可以使用 GMap.Net 中的 Winforms 东西，因为这会避免覆盖等。
我希望有人可以帮助我，或者给我看一篇解释了我没有找到的帖子。
【问题讨论】：
标签： c# wpf geometry-surface gmap.net
在本节中，我可以详细说明我是如何得出这些公式的。
让我们注意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);
}
}
【讨论】：
我使用来自 Cédric 的部分代码和来自互联网的代码修复了它。
我使用 poly.Points 和 poly.LocalPoints 修复了它。 poly.Points 是纬度和经度，而 LocalPoints 是屏幕上地图中心的点。
C# 库有一个计算距离 (km) 的函数，所以我计算了距离，然后我在 LocalPoints 中计算了距离。将 localPoints 潜入以 km 为单位的长度，然后您就知道 1 Localpoint 以 km 为单位有多长。
代码：
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 回答了我的问题并帮助我解决了我遇到的问题。
【讨论】：
只使用像 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();
}
}
}
【讨论】：