geos使用心得
缘由
因为项目需要用到 切割 与 合并 算法,百般搜索后,寻到geos。效果很好,学习经历曲折,为避免遗忘,便有此文。
采用的 goes3.5.0,使用的是 C语言版 的接口,非常标准,很好调用。
开源协议
geos 采用的是 LGPL 协议。LGPL允许商业软件通过类库引用(link)方式使用LGPL类库而不需要开源商业软件的代码。满足项目要求,可用!
五种开源协议的比较(BSD,Apache,GPL,LGPL,MIT)
http://www.ha97.com/833.html
本地编译
3.5.0 【编译通过 可用】,采用 VS2012 进行的编译。5分钟内 编译完成
主要参照这篇文章:
https://www.bbsmax.com/A/ke5jNGr75r/
-
下载geos-3.5.0,放在d:\geos350中
下载网站:http://trac.osgeo.org/geos/ -
进入Visual Studio Tools下的VS2012 开发人员命令提示,本例为
C:\Program Files (x86)\Microsoft Visual Studio 11.0>
3、依次执行如下命令
VCVARS32.BAT
cd d:\geos350
atuogen.bat
nmake /f makefile.vc
编译成功后,会在d:\geos350/src目录下生成geos.lib, geos_i.lib, geos_c_i.lib, geos.dll, geos_c.dll等五个文件
不想编译的伙伴,这里有编译好的文件,可进行下载:
https://download.csdn.net/download/lvye1221/16614981
分割功能实现
特别注意,分割前一定要进行相交的转换。【必须要求是相接的】
VS2015 GDAL c++ 开发—— GEOS 转面
https://www.bilibili.com/read/cv4426449/
QT 分割的主要参考代码
GEOSContextHandle_t handle = GEOS_init_r();
const GEOSGeometry* * geoms = new const GEOSGeometry*[polylineCount];
GEOSCoordSequence* seq = nullptr;
for (int i = 0; i < polylineCount; i++) {
seq = GEOSCoordSeq_create_r(handle, 2, 0);
GEOSCoordSeq_setX_r(handle, seq, 0, formatDouble(vecLine[i].p1().x()));
GEOSCoordSeq_setY_r(handle, seq, 0, formatDouble(vecLine[i].p1().y()));
GEOSCoordSeq_setX_r(handle, seq, 1, formatDouble(vecLine[i].p2().x()));
GEOSCoordSeq_setY_r(handle, seq, 1, formatDouble(vecLine[i].p2().y()));
// qDebug() << vecLine[i].p1() << vecLine[i].p2();
geoms[i] = (GEOSGeom_createLineString_r(handle, seq));
}
GEOSGeometry * outGeom = GEOSPolygonize_r(handle, geoms, polylineCount);
int num = GEOSGetNumGeometries_r(handle, outGeom);
// qDebug() << "num:" << num << endl;
char * s = GEOSGeomToWKT_r(handle, outGeom);
// qDebug() << "s:" << s << endl;
for (int i = 0; i < num; i++) {
const GEOSGeometry * gm = GEOSGetGeometryN_r(handle, outGeom, i);
const GEOSGeometry * line = GEOSGetExteriorRing_r(handle, gm);
int pointCount = GEOSGeomGetNumPoints_r(handle, line);
GSGraphicsPolygon* polygon = new GSGraphicsPolygonZone();
// pointCount-1 代表 不需要 后面那个相同的点
for (int j = 0; j < pointCount-1; j++) {
GEOSGeometry * point = GEOSGeomGetPointN_r(handle, line, j);
double x = 0.0;
double y = 0.0;
GEOSGeomGetX_r(handle, point, &x);
GEOSGeomGetY_r(handle, point, &y);
polygon->addPoint(QPointF(x,y));
}
vecPolygon.append(polygon);
}
delete[] geoms;
// 结束句柄
GEOS_finish_r(handle);
重点解析
GEOSPolygonize_r 要求的线段要相接,所以要对图形进行处理,满足输入要求。
GEOSGetNumGeometries_r 输出图形的数量
GEOSGeomToWKT_r 输出具体图形的信息
产生相接图形的流程
合并功能实现
QT 合并的主要参考代码
GEOSContextHandle_t handle = GEOS_init_r();
GEOSCoordSequence* seq = nullptr;
GEOSGeometry* boxLr = nullptr;
GEOSGeometry* boxP = nullptr;
GEOSGeometry* boxLr2 = nullptr;
GEOSGeometry* boxP2 = nullptr;
int n1 = polygon1->m_vecPoint.size();
seq = GEOSCoordSeq_create_r(handle, n1+1, 0);
for (int i = 0; i < n1; i++) {
QPointF point = polygon1->m_vecPoint[i];
GEOSCoordSeq_setX_r(handle, seq, i, formatDouble(point.x()));
GEOSCoordSeq_setY_r(handle, seq, i, formatDouble(point.y()));
}
GEOSCoordSeq_setX_r(handle, seq, n1, formatDouble(polygon1->m_vecPoint.first().x()));
GEOSCoordSeq_setY_r(handle, seq, n1, formatDouble(polygon1->m_vecPoint.first().y()));
// 创建环线
boxLr = GEOSGeom_createLinearRing_r(handle, seq);
//创建面
boxP = GEOSGeom_createPolygon_r(handle, boxLr, nullptr, 0);
int n2 = polygon2->m_vecPoint.size();
seq = GEOSCoordSeq_create_r(handle, n2+1, 0);
for (int i = 0; i < n2; i++) {
QPointF point = polygon2->m_vecPoint[i];
GEOSCoordSeq_setX_r(handle, seq, i, formatDouble(point.x()));
GEOSCoordSeq_setY_r(handle, seq, i, formatDouble(point.y()));
}
GEOSCoordSeq_setX_r(handle, seq, n2, formatDouble(polygon2->m_vecPoint.first().x()));
GEOSCoordSeq_setY_r(handle, seq, n2, formatDouble(polygon2->m_vecPoint.first().y()));
// 创建环线
boxLr2 = GEOSGeom_createLinearRing_r(handle, seq);
//创建面
boxP2 = GEOSGeom_createPolygon_r(handle, boxLr2, nullptr, 0);
GEOSGeometry * outGeom = GEOSUnion_r(handle, boxP, boxP2);
int num = GEOSGetNumGeometries_r(handle, outGeom);
char * s = GEOSGeomToWKT_r(handle, outGeom);
if (num > 1) {
// 说明生成失败
return nullptr;
}
const GEOSGeometry * gm = GEOSGetGeometryN_r(handle, outGeom, 0);
const GEOSGeometry * line = GEOSGetExteriorRing_r(handle, gm);
int pointCount = GEOSGeomGetNumPoints_r(handle, line);
GSGraphicsPolygon* polygon = new GSGraphicsPolygonZone();
// pointCount-1 代表 不需要 后面那个相同的点
for (int j = 0; j < pointCount-1; j++) {
GEOSGeometry * point = GEOSGeomGetPointN_r(handle, line, j);
double x = 0.0;
double y = 0.0;
GEOSGeomGetX_r(handle, point, &x);
GEOSGeomGetY_r(handle, point, &y);
polygon->addPoint(QPointF(x,y));
}
心得体会
扎实的基本功,不断尝试,探索研究,耐下性子,终有所成。
看到网上文章,有些没有输出图形信息和个数,这时候,就不断寻找资料,看到说要看示例代码程序,再大胆尝试,就有思路和想法了,最终解决了问题。
当输入参数时,要求相接,所以凭自己想法写了一个相接函数,流程如上图所示,还算解决了问题
其实作为私企的开发者,重点是解决问题,不再是重复造轮子,如何让代码更简洁、易维护、项目更稳定 是最为重要的。如有时间,深入研究实现算法也是有价值的,只是优先级不要那么高。总之,提升解决问题的能力是一直要注意的。
参考资料
GEOS库在windows中的编译和测试(vs2012)
https://www.cnblogs.com/denny402/p/4965213.html
关于GEOS库配置与安装
https://blog.csdn.net/qhqlnannan/article/details/110000955
JTS 相关问题解答
https://locationtech.github.io/jts/jts-faq.html#A1
叠置算法(5):拓扑构建算法
http://www.whudj.cn/?p=970
求交集的
要注意: buffer
https://blog.csdn.net/songyu0120/article/details/104489282
C 接口的测试
https://blog.csdn.net/dongyesang/article/details/78979287
多段线切割多边形
https://www.jianshu.com/p/f252f3eb69f8
geos 多边形分割
https://blog.totoroxiao.com/geo-polygon-split-union/
CGAL,Computational Geometry Algorithms Library,计算几何算法库
https://www.freesion.com/article/47271023281/
Boost.Geometry介绍
https://blog.csdn.net/yanfeng1022/article/details/100000741
geos 标准c接口文档
https://trac.osgeo.org/geos