柚子快报激活码778899分享：geos使用心得

http://www.51969.com/

缘由

因为项目需要用到 切割 与 合并 算法，百般搜索后，寻到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

柚子快报激活码778899分享：geos使用心得

http://www.51969.com/

查看原文