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/

1. 下载geos-3.5.0，放在d:\geos350中
   下载网站：http://trac.osgeo.org/geos/

2. 进入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