geos 学习笔记-1

最近负责地形系统的数据处理，空间几何的关联操作需求较多，得知geos比较成熟，于是决定学习吸收。

下面是我刚编译完geos源码后写的测试代码。并且用2个简单多边形的相交测试对比了geos的操作效率，代码如下：

#include <vector>
#include <string>
#include <memory>
#include <list>
#include "geos.h"
#include <stdio.h>
#include <functional>
#include<iomanip>
#include <time.h>
#include "tMath.h"
using namespace std;
using namespace geos;
using namespace geom;

#include <functional>
class ScopeGuard
{
public:
    explicit ScopeGuard(std::function<void()> onExit): __onExit(onExit),__dismissed(false)
    {
    }
    ~ScopeGuard()
    {
        if(!__dismissed)
            __onExit();
    }
    void Dismiss(bool dismissed=true)
    {
        __dismissed=dismissed;
    }
    ScopeGuard(ScopeGuard const&)=delete;
    ScopeGuard& operator=(ScopeGuard const&)=delete;
private:
    std::function<void()> __onExit;
    bool __dismissed;
};

class Test0
{
public:
    void init()
    {
        _geoFactory.reset(new GeometryFactory);
    }
private:
    Polygon* _createPolygon(const geom::CoordinateArraySequence& points)
    {
        LinearRing *shell=_createRing(points);
        return _createPolygon(shell);
    }

    Polygon* _createPolygon(LinearRing* shell,std::vector<Geometry*> *holes = NULL)
    {
        return _geoFactory->createPolygon(shell,holes);
    }

    LinearRing* _createRing(const geom::CoordinateArraySequence& points)const
    {
        return _geoFactory->createLinearRing(points);
    }

    LineString* _createLine(const geom::CoordinateArraySequence& points) const
    {
        return _geoFactory->createLineString(points);
    }

    double _distGeo(const Geometry *geoA,const Geometry *geoB)const
    {
        return geoA->distance(geoB);
    }

    Geometry* _intersectGeo(const Geometry *geoA,const Geometry *geoB)const
    {
        return geoA->intersection(geoB);
    }

    Geometry* _unionGeo(const Geometry *geoA,const Geometry *geoB)const
    {
        return geoA->Union(geoB);
    }

    Geometry* _diffGeo(const Geometry *geoA,const Geometry *geoB)const
    {
        return geoA->difference(geoB);
    }

    void _coutGeo(const Geometry *geo)
    {
        if(geo->isEmpty())
        {
            cout<<"empty"<<std::endl;
            return;
        }
            
        CoordinateArraySequence *vCoords(dynamic_cast<CoordinateArraySequence*>(geo->getCoordinates()));
        const std::vector<Coordinate>* points=vCoords->toVector();
        std::vector<Coordinate>::const_iterator itPt=points->begin();
        int iFlag=1;
        while(itPt!=points->end())
        {
            Coordinate pt=*itPt++;
            cout<<setprecision(10)<<pt.x<<" "<<pt.y;
            if((iFlag++%4)==0)
                cout<<endl;
            else
                cout<<"  ";
        }

        const Polygon *po=dynamic_cast<const Polygon*>(geo);
        if (po != NULL)
        {
            cout<<setprecision(10)<<"area:"<<po->getArea();
        }

        cout<<std::endl<<endl;
    }

    void _coutSharePt(const std::vector<std::pair<double,double>>& sharedPoints)
    {
        int iFlag=1;
        std::vector<std::pair<double,double>>::const_iterator itPt=sharedPoints.begin();
        while (itPt != sharedPoints.end())
        {
            const std::pair<double,double>& pt=*itPt++;
            cout<<setprecision(10)<<pt.first<<" "<<pt.second;

            if((iFlag++%4)==0)
                cout<<endl;
            else
                cout<<"  ";
        }

        cout<<std::endl<<endl;
    }

    void _time(int use)
    {
        int hour=use/(CLOCKS_PER_SEC*3600);
        int minute=(use%(CLOCKS_PER_SEC*3600))/(CLOCKS_PER_SEC*60);
        int second=(use%(CLOCKS_PER_SEC*60))/CLOCKS_PER_SEC;
        int millSecond=use%1000;

        cout<<"----------------------------------------------------------"<<std::endl;
        cout<<"总耗时："<<hour<<"小时"<<minute<<"分钟"<<second<<"秒"<<millSecond<<"微秒"<<std::endl;
    }
public:
    void mainExe()
    {
        ScopeGuard onExit([&]
        {
            char cTemp=getchar();
        });
        //---------------------------------
        //用geos作相交并输出交点
        CoordinateArraySequence vCoords;
        vCoords.clear();
        vCoords.add(Coordinate(118.0050000000,32.0060000000));
        vCoords.add(Coordinate(118.0060000000,32.0060000000));
        vCoords.add(Coordinate(118.0060000000,32.0050000000));
        vCoords.add(Coordinate(118.0050000000,32.0050000000));
        vCoords.add(Coordinate(118.0050000000,32.0060000000));
        LinearRing *tileRing=_createRing(vCoords);
        if (tileRing->isClosed()==false)
        {
            cout<<"tileRing is not closed"<<std::endl;
            return;
        }
        Polygon *tile=_createPolygon(tileRing);
        _coutGeo(tile);

        vCoords.clear();
        vCoords.add(Coordinate(118.0050000000,32.0055000000));
        vCoords.add(Coordinate(118.0055000000,32.0057500000));
        vCoords.add(Coordinate(118.0062500000,32.0050000000));
        vCoords.add(Coordinate(118.0055000000,32.0047500000));
        vCoords.add(Coordinate(118.0050000000,32.0055000000));
        LinearRing *holeRing=_createRing(vCoords);
        if(holeRing->isClosed()==false)
        {
            cout<<"tileRing is not closed"<<std::endl;
            return;
        }
        Polygon *hole=_createPolygon(vCoords);
        _coutGeo(hole);

        int hour(0),minute(0),second(0);
        clock_t begin,end,use;
        
        begin=clock();
        Geometry *geo=tile->intersection(hole);
        end=clock();
        use=end-begin;
        _time(use);
        _coutGeo(geo);

        //-------------------------------------
        //用逐边裁剪处理同样的2个简单多边形
        Rect tileRec(118.0050000000,118.0060000000,32.0050000000,32.0060000000);
        std::vector<std::pair<double,double>> holePts;
        holePts.push_back(make_pair(118.0050000000,32.0055000000));
        holePts.push_back(make_pair(118.0055000000,32.0057500000));
        holePts.push_back(make_pair(118.0062500000,32.0050000000));
        holePts.push_back(make_pair(118.0055000000,32.0047500000));
        std::vector<std::pair<double,double>> sharedPoints;
        begin=clock();
        if(false==sutherlangCut(tileRec,holePts,sharedPoints))
        {
            cout<<"sutherlangCut-error"<<std::endl;
            return;
        }
        for(int i=0;i<19;i++)
            sutherlangCut(tileRec,holePts,sharedPoints);
        
        end=clock();
        use=end-begin;
        _time(use);
        _coutSharePt(sharedPoints);
 
    }
private:
    shared_ptr<GeometryFactory> _geoFactory;
};

void main()
{
    Test0 test0;
    test0.init();
    test0.mainExe();
}

结果如下：

我的相交测试因为只针对不带洞的凸多边形，所以实际效率会比geos的操作要快很多，我循环了同样的操作20次。

不过考虑到geos的范用性，该对比并不能作为判断geos效率的依据，只是个学习代码而已。

我的工程代码就2个文件，下面提供下载地址，对geos不是和在下一样刚接触的就不用看了额。

 

工程文件2个: 文件.rar