如何用C#快速筛选多面体内外的地理坐标？

问题背景 #

使用C#结合NetTopologySuite筛选大型国家精确边界内的135万条点数据，现有代码功能正常但耗时长达一天，已尝试Parallel.ForEach优化但效果不佳，寻求性能提升方案。

原始代码 #

var path = "fr.shp"; // 需要精确匹配大型国家边界

var reader = new ShapeDataReader(path);
var mbr = reader.ShapefileBounds;
var result = reader.ReadByMBRFilter(mbr);

var polygons = new List<Geometry>();

using (var coll = result.GetEnumerator())
{
    while (coll.MoveNext())
    {
        var item = coll.Current;

        if (item == null)
        {
            continue;
        }

        polygons.Add(item.Geometry);
    }
}

var polygon = new GeometryCombiner(polygons).Combine();

var points = new List<Point>();
List<Point> pointsToFilterWithBorders; // 从数据库加载，共1,350,000条点数据

Parallel.ForEach(pointsToFilterWithBorders, point =>
{
    if (polygon.Contains(point))
        points.Add(point);
});

优化方案 #

1. 给多边形构建空间索引 #

NetTopologySuite的STRtree可通过MBR快速过滤掉不可能在边界内的点，大幅减少Contains的计算次数，这是核心优化手段：

// 构建空间索引
var spatialIndex = new STRtree<Geometry>();
if (polygon is MultiPolygon multiPoly)
{
    foreach (var singlePoly in multiPoly.Geometries)
    {
        spatialIndex.Insert(singlePoly.EnvelopeInternal, singlePoly);
    }
}
else if (polygon is Polygon singlePoly)
{
    spatialIndex.Insert(singlePoly.EnvelopeInternal, singlePoly);
}

// 并行筛选时先通过索引缩小判断范围
var concurrentPoints = new ConcurrentBag<Point>();
var parallelOptions = new ParallelOptions { MaxDegreeOfParallelism = Environment.ProcessorCount };

Parallel.ForEach(pointsToFilterWithBorders, parallelOptions, point =>
{
    // 查询可能包含该点的多边形子集
    var candidatePolygons = spatialIndex.Query(point.EnvelopeInternal).Cast<Geometry>();
    if (candidatePolygons.Any(p => p.Contains(point)))
    {
        concurrentPoints.Add(point);
    }
});

var filteredPoints = concurrentPoints.ToList();

2. 修复并行处理的线程安全问题 #

原始代码直接操作非线程安全的List<Point>，会引发锁竞争和数据错误。改用ConcurrentBag<Point>或手动分区合并结果，避免不必要的性能损耗。

3. 优化几何合并结果类型 #

GeometryCombiner.Combine()返回的可能是GeometryCollection，转换为MultiPolygon可避免后续判断时的类型转换开销：

var combinedGeometry = new GeometryCombiner(polygons).Combine();
var countryPolygon = combinedGeometry switch
{
    MultiPolygon mp => mp,
    Polygon p => new MultiPolygon(new[] { p }),
    _ => throw new InvalidOperationException("无法转换为多边形类型")
};

4. 数据库端预过滤（若适用） #

如果点数据存储在支持空间查询的数据库（如PostGIS、SQL Server Spatial），先在数据库层面用国家边界的MBR过滤掉大部分无关点，再将剩余点导入客户端做精确判断，能大幅减少客户端处理的数据量。

5. 可选：几何简化（允许精度损失时） #

若业务场景允许轻微精度损失，用DouglasPeuckerSimplifier减少多边形顶点数，降低Contains的计算复杂度：

// 0.001为容差，值越大简化程度越高，精度损失越大
var simplifiedPolygon = DouglasPeuckerSimplifier.Simplify(countryPolygon, 0.001);

内容的提问来源于stack exchange，提问作者Denis H