gdal–矢量求交

前言

最近在帮同学写几个小程序，包括矢量数据求交（Intersection）。网上相关文章很少，所以我参考GDAL官网写了一个完整的程序（裁剪、合并、擦除、更新等功能只要修改代码的求交函数部分，下面会有说明）。

代码

#include "gdal_priv.h"  
#include "ogrsf_frmts.h" //for ogr 

bool VectorIntersection(const char *pszSrcShp, const char *pszMethodShp, const char *pszDstShp, const char* pszFormat);

int main()
{
    const char *pszSrcFile = "C:\\Users\\liuwei\\Desktop\\CreateIm\\shpManage\\SrcLayer.shp";//原始文件
    const char *pszMethodFile = "C:\\Users\\liuwei\\Desktop\\CreateIm\\shpManage\\MethodLayer.shp";//用来求交的文件
    const char *pszOutFile = "C:\\Users\\liuwei\\Desktop\\shp\\Intersection.shp";//求交后的结果文件
    VectorIntersection(pszSrcFile, pszMethodFile, pszOutFile, "ESRI Shapefile");//矢量数据求交
    return 0;
}

bool VectorIntersection(const char *pszSrcShp, const char *pszMethodShp, const char *pszDstShp, const char* pszFormat)
{
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

    //打开原始数据
    GDALDataset *poSrcDS;
    poSrcDS = (GDALDataset*)GDALOpenEx(pszSrcShp, GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (poSrcDS == NULL)
    {
        printf("Open failed.\n");
        exit(1);
    }

    //打开用来求交的数据
    GDALDataset *poMethodDS;
    poMethodDS = (GDALDataset*)GDALOpenEx(pszMethodShp, GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (poMethodDS == NULL)
    {
        printf("Open failed.\n");
        exit(1);
    }

    // 根据Shapefile驱动创建输出矢量文件 
    GDALDriver *poDriver;
    poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);\
    if (poDriver == NULL)
    {
        printf("%s driver not available.\n", pszFormat);
        exit(1);
    }
    //根据文件名创建输出矢量数据集 
    GDALDataset* poDstDS = poDriver->Create(pszDstShp, 0, 0, 0, GDT_Unknown, NULL);
    if (poDstDS == NULL)
    {
        printf("Creation of output file failed.\n");
        exit(1);
    }

    OGRLayer *poSrcLayer = poSrcDS->GetLayer(0);//得到原始数据第一个图层
    if (poSrcLayer == NULL)
    {
        printf("Creation of layer failed.\n");
        GDALClose(poSrcDS); //关闭文件  
        GDALClose(poMethodDS);
        GDALClose(poDstDS);
        exit(1);
    }

    OGRLayer *poMethodLayer = poMethodDS->GetLayer(0);
    if (poMethodLayer == NULL)
    {
        printf("Creation of layer failed.\n");
        GDALClose(poSrcDS); //关闭文件  
        GDALClose(poMethodDS);
        GDALClose(poDstDS);
        exit(1);
    }

    // 定义空间参考，与原始矢量数据相同 
    OGRSpatialReference *pSRS = poSrcLayer->GetSpatialRef();
    OGRLayer *poDstLayer = poDstDS->CreateLayer("NewLayer", pSRS, wkbPolygon, NULL);//创建图层
    if (poDstLayer == NULL)
    {
        printf("Creation of layer failed.\n");
        GDALClose(poSrcDS); //关闭文件  
        GDALClose(poMethodDS);
        GDALClose(poDstDS);
        exit(1);
    }

    poSrcLayer->Intersection(poMethodLayer, poDstLayer, NULL, GDALTermProgress, NULL);//图层求交

    GDALClose(poSrcDS); //关闭文件  
    GDALClose(poMethodDS);
    GDALClose(poDstDS);

    return true;
}

说明

1.实现矢量数据求交需要用到GDAL+Geos，关于环境配置，网上内容太杂，可参考以下博客：GDAL+GEOS(注意：先编译GEOS再编译GDAL，如果之前编译过GDAL了，建议重新解压GDAL再编译，不要在原来的基础上修改nmake.opt后再编译，那样可能还会报错：ERROR 6: GEOS support not enabled.我开始就是因为这报错)。
2.OGR库要添加GEOS库的支持，需要添加geos库的部分头文件的路径，主要是两个文件，一个是capi/geos_c.h，一个是include/geos/version.h。因此添加两个路径capi/和source/headers/就可以了。
3.GDAL集成GEOS（具体步骤参考我上面推荐的博客）需要找到GDAL源代码目录下的nmake.opt文件，用编辑器打开nmake.opt搜索goes，将

# Uncomment for GEOS support (GEOS >= 3.1.0 required)
#GEOS_DIR=C:/warmerda/geos
#GEOS_CFLAGS=-I$(GEOS_DIR)/capi -I$(GEOS_DIR)/source/headers -DHAVE_GEOS
#GEOS_LIB    =$(GEOS_DIR)/source/geos_c_i.lib

修改为：

# Uncomment for GEOS support (GEOS >= 3.1.0 required)
GEOS_DIR=D:\develop\geos-3.4.3
GEOS_CFLAGS = -I$(GEOS_DIR)/capi -I$(GEOS_DIR)/include -DHAVE_GEOS
GEOS_LIB     = $(GEOS_DIR)/src/geos_c_i.lib

4.注意：编译GEOS和GDAL后，记得将GEOS库文件夹下src文件夹中geos_c.dll文件，拷贝到GDAL编译后存放的目录下的bin文件夹中gdal202.dll（202是版本号）的同级目录下，否则会提示你找不到geos_c.dll文件。
5.矢量数据裁剪、合并、擦除、更新等功能只需要修改用到的图层操作函数，相应的函数定义如下：
可参考GDAL官网OGRLayer类中函数说明。将代码中poSrcLayer->Intersection(poMethodLayer, poDstLayer, NULL, GDALTermProgress, NULL)直接修改即可。

转载自：https://blog.csdn.net/secyb/article/details/80246105