我正在寻找一种计算多边形下表面的方法。
我想要完成的事情是使用我的程序的用户可以创建一个多边形来标记他的属性。现在我想知道表面积是多少，这样我就可以告诉用户他的财产有多大。
单位 m² 或 km² 或公顷。
多边形的点具有纬度和经度。
我将 C# 与 WPF 和 GMap.NET 一起使用。该地图位于其中，WindowsFormHost因此我可以使用 GMap.Net 中的 Winforms 东西，因为这会避免叠加等。
我希望有人可以帮助我或向我展示一个解释我没有找到的帖子。
问问题
4600 次
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和操作。MultiplyCrossProductDotProduct
就您而言，您有纬度和经度。然后，您不在 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);
}
}
于 2013-07-21T19:20:34.550 回答
1
我用 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 回答了我的问题并帮助我解决了我遇到的问题。
于 2013-07-22T16:09:39.900 回答
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();
}
}
}
于 2013-09-16T14:08:52.723 回答