第11章-空间查询与分析
第11章：空间查询与分析
11.1 空间查询基础
11.1.1 范围查询
// 创建查询范围
var queryEnvelope = new Envelope(116.0, 117.0, 39.0, 40.0);
// 执行范围查询
var featureDataSet = new FeatureDataSet();
vectorLayer.DataSource.ExecuteIntersectionQuery(queryEnvelope, featureDataSet);
// 获取结果
if (featureDataSet.Tables.Count > 0)
{
var table = featureDataSet.Tables[0];
Console.WriteLine($"查询到 {table.Rows.Count} 个要素");
foreach (FeatureDataRow row in table.Rows)
{
Console.WriteLine($"ID: {row["ID"]}, Name: {row["NAME"]}");
Console.WriteLine($"几何类型: {row.Geometry.GeometryType}");
}
}
11.1.2 几何查询
// 创建查询几何
var factory = new GeometryFactory();
var queryPoint = factory.CreatePoint(new Coordinate(116.4, 39.9));
var queryBuffer = queryPoint.Buffer(0.1); // 0.1度缓冲区
// 执行几何查询
var featureDataSet = new FeatureDataSet();
vectorLayer.DataSource.ExecuteIntersectionQuery(queryBuffer, featureDataSet);
11.1.3 点击查询
public class PointQueryHelper
{
private readonly Map _map;
public PointQueryHelper(Map map)
{
_map = map;
}
public IEnumerable<FeatureDataRow> QueryAtPoint(PointF screenPoint, double tolerance = 5)
{
// 转换为地理坐标
var worldPoint = _map.ImageToWorld(screenPoint);
// 创建查询范围（考虑容差）
double worldTolerance = tolerance * _map.PixelWidth;
var queryEnvelope = new Envelope(
worldPoint.X - worldTolerance,
worldPoint.X + worldTolerance,
worldPoint.Y - worldTolerance,
worldPoint.Y + worldTolerance);
var results = new List<FeatureDataRow>();
// 查询所有矢量图层
foreach (var layer in _map.Layers.OfType<VectorLayer>().Where(l => l.Enabled))
{
var ds = new FeatureDataSet();
layer.DataSource.ExecuteIntersectionQuery(queryEnvelope, ds);
if (ds.Tables.Count > 0)
{
results.AddRange(ds.Tables[0].Rows.Cast<FeatureDataRow>());
}
}
return results;
}
}
11.2 属性查询
11.2.1 使用 LINQ 查询
// 获取所有要素
var allFeatures = GetAllFeatures(vectorLayer);
// LINQ 过滤
var filtered = allFeatures
.Where(f => Convert.ToDouble(f["POPULATION"]) > 1000000)
.OrderByDescending(f => Convert.ToDouble(f["POPULATION"]))
.Take(10);
foreach (var feature in filtered)
{
Console.WriteLine($"{feature["NAME"]}: {feature["POPULATION"]}");
}
// 辅助方法
private IEnumerable<FeatureDataRow> GetAllFeatures(VectorLayer layer)
{
var ds = new FeatureDataSet();
layer.DataSource.ExecuteIntersectionQuery(layer.Envelope, ds);
if (ds.Tables.Count > 0)
{
return ds.Tables[0].Rows.Cast<FeatureDataRow>();
}
return Enumerable.Empty<FeatureDataRow>();
}
11.2.2 SQL 风格查询
public class FeatureQuery
{
public static IEnumerable<FeatureDataRow> Where(IEnumerable<FeatureDataRow> features,
Func<FeatureDataRow, bool> predicate)
{
return features.Where(predicate);
}
public static IEnumerable<FeatureDataRow> Select(IEnumerable<FeatureDataRow> features,
string[] columns)
{
// 投影查询（只返回指定字段）
foreach (var feature in features)
{
var newRow = feature.Table.NewRow() as FeatureDataRow;
newRow.Geometry = feature.Geometry;
foreach (var col in columns)
{
if (feature.Table.Columns.Contains(col))
{
newRow[col] = feature[col];
}
}
yield return newRow;
}
}
}
// 使用
var results = FeatureQuery.Where(allFeatures, f =>
f["TYPE"]?.ToString() == "City" &&
Convert.ToDouble(f["POPULATION"] ?? 0) > 500000);
11.3 空间关系判断
11.3.1 基本空间关系
using NetTopologySuite.Geometries;
var factory = new GeometryFactory();
var polygon = factory.CreatePolygon(new Coordinate[]
{
new Coordinate(0, 0),
new Coordinate(10, 0),
new Coordinate(10, 10),
new Coordinate(0, 10),
new Coordinate(0, 0)
});
var point = factory.CreatePoint(new Coordinate(5, 5));
var outsidePoint = factory.CreatePoint(new Coordinate(15, 15));
// 空间关系判断
Console.WriteLine($"包含: {polygon.Contains(point)}"); // true
Console.WriteLine($"在内部: {point.Within(polygon)}"); // true
Console.WriteLine($"相交: {polygon.Intersects(point)}"); // true
Console.WriteLine($"分离: {polygon.Disjoint(outsidePoint)}"); // true
Console.WriteLine($"相邻: {polygon.Touches(point)}"); // false
11.3.2 空间关系查询
public class SpatialRelationQuery
{
public static IEnumerable<FeatureDataRow> Contains(
IEnumerable<FeatureDataRow> features, Geometry queryGeometry)
{
return features.Where(f => f.Geometry.Contains(queryGeometry));
}
public static IEnumerable<FeatureDataRow> Within(
IEnumerable<FeatureDataRow> features, Geometry queryGeometry)
{
return features.Where(f => f.Geometry.Within(queryGeometry));
}
public static IEnumerable<FeatureDataRow> Intersects(
IEnumerable<FeatureDataRow> features, Geometry queryGeometry)
{
return features.Where(f => f.Geometry.Intersects(queryGeometry));
}
public static IEnumerable<FeatureDataRow> WithinDistance(
IEnumerable<FeatureDataRow> features, Geometry queryGeometry, double distance)
{
return features.Where(f => f.Geometry.Distance(queryGeometry) <= distance);
}
}
11.4 空间分析
11.4.1 缓冲区分析
// 单一缓冲区
var point = factory.CreatePoint(new Coordinate(116.4, 39.9));
var buffer = point.Buffer(0.1); // 0.1 度缓冲区
// 不同端点样式
var lineBuffer1 = line.Buffer(10, EndCapStyle.Round); // 圆形端点
var lineBuffer2 = line.Buffer(10, EndCapStyle.Flat); // 平头端点
var lineBuffer3 = line.Buffer(10, EndCapStyle.Square); // 方形端点
// 多环缓冲区
public static List<Geometry> CreateMultiRingBuffer(Geometry geometry, double[] distances)
{
var buffers = new List<Geometry>();
Geometry previousBuffer = geometry;
foreach (var distance in distances.OrderBy(d => d))
{
var buffer = geometry.Buffer(distance);
var ring = buffer.Difference(previousBuffer);
buffers.Add(ring);
previousBuffer = buffer;
}
return buffers;
}
// 使用
var rings = CreateMultiRingBuffer(centerPoint, new[] { 1.0, 2.0, 3.0, 5.0 });
11.4.2 叠加分析
// 交集
var intersection = polygon1.Intersection(polygon2);
// 并集
var union = polygon1.Union(polygon2);
// 差集
var difference = polygon1.Difference(polygon2);
// 对称差
var symDifference = polygon1.SymmetricDifference(polygon2);
// 批量合并
public static Geometry UnionAll(IEnumerable<Geometry> geometries)
{
return geometries.Aggregate((current, next) => current.Union(next));
}
// 使用 CascadedPolygonUnion 提高性能
using NetTopologySuite.Operation.Union;
var geometryList = new List<Geometry> { polygon1, polygon2, polygon3 };
var unionResult = CascadedPolygonUnion.Union(geometryList);
11.4.3 距离计算
// 两点距离
var distance = point1.Distance(point2);
// 最近点
var nearestPoints = NetTopologySuite.Operation.Distance.DistanceOp
.NearestPoints(geometry1, geometry2);
// 查找最近要素
public static FeatureDataRow FindNearest(IEnumerable<FeatureDataRow> features,
Geometry queryGeometry)
{
return features
.OrderBy(f => f.Geometry.Distance(queryGeometry))
.FirstOrDefault();
}
// 查找距离范围内的要素
public static IEnumerable<FeatureDataRow> FindWithinDistance(
IEnumerable<FeatureDataRow> features, Geometry queryGeometry, double maxDistance)
{
return features
.Where(f => f.Geometry.Distance(queryGeometry) <= maxDistance)
.OrderBy(f => f.Geometry.Distance(queryGeometry));
}
11.4.4 面积和长度计算
// 面积（坐标系单位）
double area = polygon.Area;
// 长度（坐标系单位）
double length = lineString.Length;
// 周长
double perimeter = polygon.Length;
// 计算地理坐标系下的实际距离（使用 Haversine 公式）
public static double CalculateDistance(Coordinate point1, Coordinate point2)
{
const double R = 6371000; // 地球半径（米）
double lat1 = point1.Y * Math.PI / 180;
double lat2 = point2.Y * Math.PI / 180;
double deltaLat = (point2.Y - point1.Y) * Math.PI / 180;
double deltaLon = (point2.X - point1.X) * Math.PI / 180;
double a = Math.Sin(deltaLat / 2) * Math.Sin(deltaLat / 2) +
Math.Cos(lat1) * Math.Cos(lat2) *
Math.Sin(deltaLon / 2) * Math.Sin(deltaLon / 2);
double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
return R * c; // 米
}
// 计算地理坐标系下的面积
public static double CalculateGeographicArea(Polygon polygon)
{
// 使用投影转换到平面坐标系计算
// 或使用球面几何公式
throw new NotImplementedException("需要更复杂的球面计算");
}
11.5 空间索引
11.5.1 使用 STRtree
using NetTopologySuite.Index.Strtree;
// 创建空间索引
var index = new STRtree<FeatureDataRow>();
// 添加要素到索引
foreach (FeatureDataRow feature in featureTable.Rows)
{
index.Insert(feature.Geometry.EnvelopeInternal, feature);
}
// 构建索引
index.Build();
// 查询
var queryEnvelope = new Envelope(116.0, 117.0, 39.0, 40.0);
var results = index.Query(queryEnvelope);
// 精确过滤
var preciseResults = results
.Where(f => f.Geometry.Intersects(queryGeometry))
.ToList();
11.5.2 使用 Quadtree
using NetTopologySuite.Index.Quadtree;
var quadtree = new Quadtree<FeatureDataRow>();
// 添加要素
foreach (FeatureDataRow feature in featureTable.Rows)
{
quadtree.Insert(feature.Geometry.EnvelopeInternal, feature);
}
// 查询
var results = quadtree.Query(queryEnvelope);
11.6 空间分析示例
11.6.1 点在面内查询
public class PointInPolygonQuery
{
private readonly STRtree<FeatureDataRow> _polygonIndex;
public PointInPolygonQuery(IEnumerable<FeatureDataRow> polygons)
{
_polygonIndex = new STRtree<FeatureDataRow>();
foreach (var polygon in polygons)
{
_polygonIndex.Insert(polygon.Geometry.EnvelopeInternal, polygon);
}
_polygonIndex.Build();
}
public FeatureDataRow FindContainingPolygon(Point point)
{
var candidates = _polygonIndex.Query(point.EnvelopeInternal);
return candidates.FirstOrDefault(p => p.Geometry.Contains(point));
}
public Dictionary<string, int> CountPointsInPolygons(IEnumerable<Point> points,
string polygonIdField)
{
var counts = new Dictionary<string, int>();
foreach (var point in points)
{
var polygon = FindContainingPolygon(point);
if (polygon != null)
{
string id = polygon[polygonIdField]?.ToString() ?? "";
if (!counts.ContainsKey(id))
counts[id] = 0;
counts[id]++;
}
}
return counts;
}
}
11.6.2 最近邻分析
public class NearestNeighborAnalysis
{
public static Dictionary<FeatureDataRow, FeatureDataRow> FindNearestNeighbors(
IEnumerable<FeatureDataRow> sourceFeatures,
IEnumerable<FeatureDataRow> targetFeatures)
{
var results = new Dictionary<FeatureDataRow, FeatureDataRow>();
var targetIndex = new STRtree<FeatureDataRow>();
foreach (var target in targetFeatures)
{
targetIndex.Insert(target.Geometry.EnvelopeInternal, target);
}
targetIndex.Build();
foreach (var source in sourceFeatures)
{
var envelope = source.Geometry.EnvelopeInternal;
envelope.ExpandBy(1); // 初始搜索范围
FeatureDataRow nearest = null;
double minDistance = double.MaxValue;
while (nearest == null)
{
var candidates = targetIndex.Query(envelope);
foreach (var candidate in candidates)
{
var dist = source.Geometry.Distance(candidate.Geometry);
if (dist < minDistance)
{
minDistance = dist;
nearest = candidate;
}
}
if (nearest == null)
{
envelope.ExpandBy(envelope.Width); // 扩大搜索范围
}
}
results[source] = nearest;
}
return results;
}
}
11.7 本章小结
本章介绍了 SharpMap 的空间查询与分析功能：
- 空间查询：范围查询、几何查询、点击查询
- 属性查询：LINQ 查询、SQL 风格查询
- 空间关系：Contains、Within、Intersects 等判断
- 空间分析：缓冲区、叠加分析、距离计算
- 空间索引：STRtree、Quadtree 的使用
下一章预告：第12章将介绍 WinForms 桌面应用开发。
浙公网安备 33010602011771号