GIS开源库GEOS库学习教程(一)：编译及示例代码

1、介绍

  GEOS库是一个集合形状的拓扑关系操作实用库，简单得说，就是判断两个几何形状之间关系和对两个几何形状进行操作以形成新的几何形状的库。GEOS是仿照JTS库做的，是JTS的C++实现。下面是JTS Topology Suite (JTS) 拓扑运算函数库的介绍，同样适合于GEOS库的介绍。

  JTS 是一套用于处理几何要素拓扑关系的函数库。它提供了完整、稳定、可靠的基本二位平面线形图形运算算法实现。

1.1 简介

• JTS 包含了”JTS Test Builder“ 用于通过交互界面测试各种功能。大量的测试保证了 JTS 的可靠性。
• JTS 被大量 Java 开源地理信息软件使用，并有 GEOS 这一 C++ 移植版用于大量 C++ 软件。另有 C# 和avaScript 移植版本。
• 它尽可能实现了OpenGIS Simple Features Specification规范，所以在于GIS相关的开源世界里如Geotools、Udig等，JTS Topology Suite都得到了大量的应用，甚至可以说没有JTS Topology Suite的话，Geotools等的实现会相当的复杂。从这些角度而言，JTS Topology Suite可以看作一个几何对象中间件。

1.2. 核心功能

• 点、线、面二维几何对象
• 交、并、包含等几何运算
• 精度控制
• 几何变换操作
• 度量操作（如 2D 笛卡尔距离、豪斯道夫距离等）
• 空间运算（如缓冲、线偏移、线简化等）
• 几何数值测量（角度、长度等）
• 空间数据结构（四分树、R 树、德洛内三角化等）
• WKT、GML 等输入和输出
• 高精度运算

2、下载及编译

  下载地址：https://libgeos.org/usage/download/下载地址

  编译方式不再详述，可以自己搜索或者参照上面网址中的方法介绍编译链接地址。

注意：源码的D:\geos3.7.5\include\geos\geom\GeometryFactory.h 中geos::geom::GeometryFactory类设定为了protect保护类型，需要修改源码改为public类型再编译，否则使用起来提示错误。

  编译完成后debug/release如下所示：

3、示例代码

新建VS项目：
引用头文件和lib
VS属性lib输入选项：geos_c_i.lib;

运行后效果如下：

main.cpp

// geosProj1.cpp : 定义控制台应用程序的入口点。//#include"stdafx.h"#include<iostream>usingnamespace std;#include"geos.h"#include"geos/geom/GeometryFactory.inl"

string intersectResult(bool flag){
    string result =(flag ==true)?"相交":"不相交";return result;}voidtestIntersect(){
    cout <<"GEOS库版本为："<< GEOS_VERSION << endl;typedef Coordinate PT;
    geos::geom::GeometryFactory factory;
    CoordinateArraySequenceFactory csf;//构建第一个矩形p1
    CoordinateSequence* cs1 = csf.create(5,2);//五个2维点，第三维度z始终为0
    cs1->setAt(PT(0,0),0);
    cs1->setAt(PT(3,0),1);
    cs1->setAt(PT(3,3),2);
    cs1->setAt(PT(0,3),3);
    cs1->setAt(PT(0,0),4);//与第一个点相等，构成闭合
    LinearRing* ring1 = factory.createLinearRing(cs1);//点构成线
    Geometry* p1 = factory.createPolygon(ring1,NULL);//线构成面

    CoordinateSequence* cs2 = csf.create(5,2);//构建一个四边形p2
    cs2->setAt(PT(2,2),0);
    cs2->setAt(PT(4,5),1);
    cs2->setAt(PT(5,5),2);
    cs2->setAt(PT(5,4),3);
    cs2->setAt(PT(2,2),4);
    LinearRing * ring2 = factory.createLinearRing(cs2);
    Geometry* p2 =(factory.createPolygon(ring2,NULL));

    CoordinateSequence *cs3 =newCoordinateArraySequence();//构建一个三角形p3int xoffset =4, yoffset =4, side =2;
    cs3->add(PT(xoffset, yoffset));
    cs3->add(PT(xoffset, yoffset + side));
    cs3->add(PT(xoffset + side, yoffset + side));
    cs3->add(PT(xoffset, yoffset));
    LinearRing * ring3 = factory.createLinearRing(cs3);
    Geometry* p3 =(factory.createPolygon(ring3,NULL));bool flag12 = p1->intersects(p2);bool flag13 = p1->intersects(p3);bool flag23 = p2->intersects(p3);
    cout <<"图1与图2:"<<intersectResult(flag12)<< endl;
    cout <<"图1与图3:"<<intersectResult(flag13)<< endl;
    cout <<"图2与图3:"<<intersectResult(flag23)<< endl;}intmain(){testIntersect();}