GEOS库一些自己的见解和注意点

GEOS的一些自己的见解和注意点

GEOS的前身是JTS，JTS提供了全功能的，强大的空间操作和空间判断。
后来PostGIS缺少一套完整的空间查询操作，于是就将JTS移植成为C++版本，正式命名为GEOS。GEOS为开源库，它包括了完整的空间查询和一大部分空间操作，是从事图形操作和GIS行业开发人员经常接触的开发库。较为知名的使用GEOS的GIS软件就有QGIS，QGIS使用GEOS的c接口，c接口函数名称不会经常发生更改，具有更多的稳定性。

                          图一

图一为GEOS定义的数据类型：点 多点 线串 闭合线串 多线串 多边形 Multi多边形 几何集合

GEOS是基于九交模型来判断空间关系的，不了解九交模型的请移步http://www.cnblogs.com/denny402/p/4968201.html，其实九交模型名字很高大上，理解起来是很简单的，一个geometry分为三种空间关系，外部、内部、边界三种，空间相对关系是相对于多个几何体来说的，九交模型是对于两个几何体的空间关系来说的，也就是3*3中，例如边界对于边界相交则为多点集合，维度为0维（面在GEOS内定义为2维，线的维度则为1维，点的维度为0维）。
由于C++ 接口有API文档可以查阅，我不做太多介绍，想要查阅的请移步http://geos.osgeo.org/doxygen/
我主要介绍一下C接口，通过查看GEOS的源码发现，c接口内部调用的还是c++的接口函数，只是在外部封了一层，但GEOS的官网上为什么会说C的接口稳定？我也不是很理解，希望大神们和我探讨一下
c接口只有一个头文件geo_c.h介绍过于少，其中有需要注意的地方：
一、要注意C接口中空间操作GEOSUnion等空间分析API，这些API是具有深拷贝行为的，但create、get之类的接口，比如createPolygon（）只是获得了指向内外环的指针，是浅拷贝。
二、GEOSBoundary()是对几何体对象进行降维操作，例如获得一个多边形的内外环线，得到外环线使用GEOSGeomGetExteriorRing(),内环线是GEOSGeomGetInteriorRingN()
三、C接口的创建面的基本顺序是与C++创建顺序是相同的，首先创建一个点序列，再创建一个闭合线串（点序列首末点需相同），最后使用闭合线串创建面。
四、每一个API后都会有一个加_r的相同功能的API，通过查看源码发现，内部定义了一个空间操作上下文，用意是什么，不太确定，希望大神可以告知
五、c接口的运行时间是与几何图形的复杂度成正比关系的，具体的复杂度我没有查看源码，无法确定，但你的几何图形只要十分复杂且点数稍多，接口计算时间就会变长，大家可以参考DEM转换出的矢量图，点数在30w左右，几何图形复杂QGIS调用GEOS_C的API，合并 拓扑化 相交等控件操作时间都会变长。特别对于intersection操作，性能十分敏感。

问题：基于GEOS的空间操作的不稳定性，比如：使用QGIS操作矢量要素，多次分割一个面要素，然后合并要素，竟然合不上！查看源码基本定位到GEOS接口的不稳定性上，进一步研究发现，只要两个点的xy坐标到15位不相同，GEOS就认为这两个点不是同一个点！原因基本是因为GEOS的精度模型是很高的，同时使用双精度数据的空间判断很不稳定，但GEOS提供了精度模型类，通过demo程序设置，不起作用！最近我们小组一直在纠结这个问题，大神们有没有什么解决方法。。。。。
通过源码查看，GEOS内部有几种策略，如果需要设置精度模型，需要设置相关的策略，精度模型才能起效。
基本上可以确认，GEOS是使用double值来计算，关于空间计算的库，例如boost的geometry等，只要使用了双精度计算，都会存在数值不稳定型问题，根源再于四舍五入的问题！在高精度的计算中，这个问题显得尤为严重，可以引发这个geometry不合法，发生拓扑错误，GEOS内部也设置了许多策略来克服这个问题，但都没有很完美的解决，只能降低发生的概率，不能完全解决。
我了解的c接口需要注意的就这么多，有什么问题大家可以私信我或者在下面评论，大家一起讨论哈。
附上自己写的代码…

#include <iostream>
#include "geos_c.h"
int main()
{
    // 初始化
    GEOSContextHandle_t handle = GEOS_init_r();
    GEOSCoordSequence* seq = nullptr;
    GEOSGeometry* trian = nullptr;
    GEOSGeometry* triPolygon = nullptr;
    GEOSGeometry* boxLr = nullptr;
    GEOSGeometry* boxP = nullptr;

    // 创建三角形
    seq = GEOSCoordSeq_create_r(handle, 4, 0);
    // 赋值坐标
    // 0
    GEOSCoordSeq_setX_r(handle, seq, 0, 300.0);
    GEOSCoordSeq_setY_r(handle, seq, 0, 0.0);
    // 1
    GEOSCoordSeq_setX_r(handle, seq, 1, 400.0);
    GEOSCoordSeq_setY_r(handle, seq, 1, 0.0);
    // 2
    GEOSCoordSeq_setX_r(handle, seq, 2, 300.0);
    GEOSCoordSeq_setY_r(handle, seq, 2, 100.0);
    // 3
    GEOSCoordSeq_setX_r(handle, seq, 3, 300.0);
    GEOSCoordSeq_setY_r(handle, seq, 3, 0);
    //创建环线
    trian = GEOSGeom_createLinearRing_r(handle, seq);
    // 创建线串

    //创建三角面
    triPolygon = GEOSGeom_createPolygon_r(handle, trian, nullptr, 0);

    // 创建 面
    seq = GEOSCoordSeq_create_r(handle, 5, 0);
    // 赋值坐标
    // 0
    GEOSCoordSeq_setX_r(handle, seq, 0, -200.0);
    GEOSCoordSeq_setY_r(handle, seq, 0, -200.0);
    // 1
    GEOSCoordSeq_setX_r(handle, seq, 1, 200.0);
    GEOSCoordSeq_setY_r(handle, seq, 1, -200.0);
    // 2
    GEOSCoordSeq_setX_r(handle, seq, 2, 200.0);
    GEOSCoordSeq_setY_r(handle, seq, 2, 200.0);
    // 3
    GEOSCoordSeq_setX_r(handle, seq, 3, -200.0);
    GEOSCoordSeq_setY_r(handle, seq, 3, 200.0);
    // 4
    GEOSCoordSeq_setX_r(handle, seq, 4, -200);
    GEOSCoordSeq_setY_r(handle, seq, 4, -200);

    //创建环线
    boxLr = GEOSGeom_createLinearRing_r(handle, seq);
    //创建面
    boxP = GEOSGeom_createPolygon_r(handle, boxLr, nullptr, 0);

    // 测试空间关系
    //char contains = GEOSContains_r(handle, boxP, triPolygon);
    //char insert = GEOSIntersects_r(handle, boxP, triPolygon);
    //char within = GEOSWithin_r(handle, boxP, triPolygon);
    //std::cout << insert << std::endl;

    // 空间切割
    GEOSGeometry* geom = GEOSClipByRect_r(handle, triPolygon, 300, 0, 350, 50);
    char isValid = GEOSisValid_r(handle, geom);
    int num = GEOSGetNumCoordinates_r(handle, geom);

    //ID
    int id = GEOSGeomTypeId_r(handle, triPolygon);


    // 删除内存
    GEOSGeom_destroy_r(handle, triPolygon);
    GEOSGeom_destroy_r(handle, boxP);
    // 结束句柄
    GEOS_finish_r(handle);

    return 0;

}