本文提出的线线吸附算法的核心思路是:以参考线为基准,将待吸附线的线段投影到参考线上,并截取对应长度的一段作为结果。
具体流程如下:
数据预处理与坐标统一:将所有输入图层统一转换到EPSG:4498,该坐标系是上海地区常用的投影坐标系,能有效控制长度变形。
多部件几何体合并与拆分:将图层中所有要素的几何部件提取出来,合并到一个
MultiLineString中,再通过mergeLines()将首尾相连的线段合并为完整的线,最后拆分为单部件LineString集合。空间索引构建:为参考线创建20米缓冲区和空间索引(R-tree),加速空间查询。
候选参考线筛选:通过空间索引快速找出与待吸附线相交(缓冲区相交)的候选参考线,再通过精确的
intersects()判断筛选出真正相交的参考线。最佳匹配选择:选择与待吸附线交叠长度最长且超过待吸附线长度60%的参考线作为匹配目标。
线段吸附:
计算待吸附线的中点
将中点投影到选定的参考线上,得到投影位置
以待吸附线长度的一半为半径,在参考线上截取对应长度的线段
截取结果即为吸附后的几何体
#include "qgsmultilineString.h" QVector<QgsGeometry> mergeLayerFeatures(QgsVectorLayer* inputLayer) { QgsCoordinateReferenceSystem srcCrs = inputLayer->crs(); //上海常用4498 auto destCrs = QgsCoordinateReferenceSystem("EPSG:4498"); QgsCoordinateTransform transform(srcCrs, destCrs, QgsProject::instance()); auto refFeatureIter = inputLayer->getFeatures(); std::unique_ptr<QgsMultiLineString> multiLineString(new QgsMultiLineString()); QgsFeature refFeature; while (refFeatureIter.nextFeature(refFeature)) { if (!refFeature.hasGeometry()) { continue; } auto refFeatureGeo = refFeature.geometry(); if (refFeatureGeo.isEmpty() || !refFeatureGeo.isGeosValid()) { continue; } refFeatureGeo.transform(transform); auto refFeatureGeoParts = refFeatureGeo.constParts(); while (refFeatureGeoParts.hasNext()) { multiLineString->addGeometry(refFeatureGeoParts.next()->clone()); } } if (multiLineString->isEmpty()) { return QVector<QgsGeometry>();; } const QgsGeometry in(multiLineString.release()); //拆成QgsLineString集合 return in.mergeLines().asGeometryCollection(); } void lineAdsorption(QgsVectorLayer* adsLayer, QgsVectorLayer* refLayer) { auto adsGeos = mergeLayerFeatures(adsLayer); auto refGeos = mergeLayerFeatures(refLayer); QVector<QgsGeometry> adsResultGeos; QVector<QgsGeometry> halfPtGeos; adsResultGeos.resize(adsGeos.size()); halfPtGeos.resize(adsGeos.size()); for (int i = 0; i < adsGeos.size();i++) { adsResultGeos[i] = adsGeos[i]; halfPtGeos[i] = adsGeos[i].interpolate(adsGeos[i].length() / 2); } //吸附阈值 double adsThreshold = 20; QgsSpatialIndex refSpatialIndex; QVector<QgsGeometry> refBufferGeos; for (int i = 0; i < refGeos.size(); i++) { auto refBufferGeo = refGeos[i].buffer(adsThreshold, 20); refBufferGeos.append(refBufferGeo); refSpatialIndex.addFeature(i, refBufferGeo.boundingBox()); } #pragma omp parallel for for (int i = 0; i < adsGeos.size(); i++) { QList<QgsFeatureId> candidateIds = refSpatialIndex.intersects(adsGeos[i].boundingBox()); if (candidateIds.size() == 0) { continue; } int intersectIndex = -1; double intersectLen = -1.0; for (auto& candidateId: candidateIds) { if (!adsGeos[i].intersects(refBufferGeos[candidateId])) { continue; } auto intersectGeo = adsGeos[i].intersection(refBufferGeos[candidateId]); if (intersectGeo.length() > (adsGeos[i].length() * 0.6) && intersectGeo.length() > intersectLen) { intersectLen = intersectGeo.length(); intersectIndex = candidateId; } } if (intersectLen < 0) { continue; } //根据中点向两边延伸 auto resLocate = refGeos[intersectIndex].lineLocatePoint(refGeos[intersectIndex].nearestPoint(halfPtGeos[i])); auto sD = std::max(0.0, resLocate - adsGeos[i].length()/2); auto eD = std::min(refGeos[intersectIndex].length(), resLocate + adsGeos[i].length() / 2); const QgsLineString* lineString = qgsgeometry_cast<const QgsLineString*>(refGeos[intersectIndex].constGet()); adsResultGeos[i] = QgsGeometry(lineString->curveSubstring(sD, eD)->clone()); } //adsResultGeos[i]为吸附后结果 }本文代码在逻辑层面进行了完整的推演与验证。由于开发环境的限制,代码尚未在完整数据集上进行系统性测试。文中提供的算法流程和代码实现可作为技术参考,读者在实际部署时请根据具体业务场景进行适配和充分测试。