在GDAL中如果判断矢量中两个多边形的关系，比如说取交集、取并集、相差等操作，是需要GEOS库支持的，然而GDAL默认是不编译GEOS的，这篇就是介绍如何在GDAL中支持GEOS，并有一个例子代码。
参考链接：http://www.cnblogs.com/denny402/p/4966558.html
http://www.cnblogs.com/sansan/p/3394636.html
1 下载GEOS源码
下载地址：http://trac.osgeo.org/geos/
2 打开VS201*的命令提示工具（Visual Studio 命令提示201x）
cd 到解压的根目录：
3 执行命令：
Nmake /f makefile.vc
稍等片刻即可编译完成。
编译成功后，会在*/src目录下生成geos.lib, geos_i.lib, geos_c_i.lib, geos.dll, geos_c.dll等五个文件
4 编译GDAL
4.1 下载GDAL
4.2 修改nmake.opt
GDAL2.2.3下
64行:修改GDAL目录
599行：修改GEOS_DIR:
由于GEOS是最新版，所以目录跟默认的可能不太一样，所以修改如下：
GEOS_DIR=G:\GDAL223\geos-3.6.2.tar\geos-3.6.2
GEOS_CFLAGS = -I$(GEOS_DIR)/capi -I$(GEOS_DIR)/include -DHAVE_GEOS
GEOS_LIB = $(GEOS_DIR)/src/geos_c_i.lib
4.3 2 打开VS201*的命令提示工具（Visual Studio 命令提示201x）
cd 到GDAL源码根目录：
执行命令：
nmake /f makefile.vc
nmake /f makefile.vc install（生成bin/html/data文件夹）
namek /f makefile.vc devinstall（生成lib/include文件夹）
至此编译完成
5 测试代码
这段代码摘自：https://blog.csdn.net/u013230291/article/details/78124785
如何在VS中引入GDAL请查阅其他教程(比如这个：
http://www.cnblogs.com/sansan/p/3394636.html)。注意的是，需要将生成的geos_c.dll复制到debuge目录下。
关于GDAL1.x与2.x对于矢量处理的差别，请看这篇：
https://blog.csdn.net/liminlu0314/article/details/72961663
主要功能是求两个shp的交集：
if (OGRGeometryFactory::haveGEOS() == false)
{
cout << "GDAL库未包含GEOS库" << endl;
return 0;
}
// 打开栅格文件
GDALDataset* poSrcDS1 = (GDALDataset*)GDALOpenEx(pszSrc1File, GDAL_OF_VECTOR, NULL, NULL, NULL);
if (poSrcDS1 == NULL)
{
return 0;
}
OGRLayer* poLayer1 = poSrcDS1->GetLayer(0);
poLayer1->ResetReading();
OGRFeature *poFeature1 = poLayer1->GetNextFeature();
OGRGeometry *poGeometry1 = poFeature1->GetGeometryRef();
GDALDataset* poSrcDS2 = (GDALDataset*)GDALOpenEx(pszSrc21File, GDAL_OF_VECTOR, NULL, NULL, NULL);
if (poSrcDS2 == NULL)
{
return 0;
}
OGRLayer* poLayer2 = poSrcDS2->GetLayer(0);
poLayer2->ResetReading();
OGRFeature *poFeature2 = poLayer2->GetNextFeature();
OGRGeometry *poGeometry2= poFeature2->GetGeometryRef();
OGRGeometry *poGeometry3 = poGeometry1->Difference(poGeometry2);
// 创建输出矢量文件
GDALDriver *poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");
if (poDriver == NULL)
{
printf("%s驱动不可用！\n", "ESRI Shapefile");
GDALClose((GDALDatasetH)poSrcDS1);
GDALClose((GDALDatasetH)poSrcDS2);
return 0;
}
//根据文件名创建输出矢量文件
GDALDataset* poDstDS = poDriver->Create(pszDstFile, 0, 0, 0, GDT_Unknown, NULL);
if (poDstDS == NULL)
{
GDALClose((GDALDatasetH)poSrcDS1);
GDALClose((GDALDatasetH)poSrcDS2);
return 0;
}
// 定义空间参考，与输入图像相同
OGRSpatialReference *poSpatialRef = new OGRSpatialReference(poSrcDS1->GetProjectionRef());
if (poSpatialRef==NULL)
{
return 0;
}
OGRLayer* poLayer = poDstDS->CreateLayer("Result", poSpatialRef, wkbPolygon, NULL);
if (poDstDS == NULL)
{
GDALClose((GDALDatasetH)poSrcDS1);
GDALClose((GDALDatasetH)poSrcDS2);
delete poSpatialRef;
poSpatialRef = NULL;
return 0;
}
OGRFeatureDefn *poDefn = poLayer->GetLayerDefn();
OGRFeature *poFeatureIntersection = OGRFeature::CreateFeature(poDefn);
poFeatureIntersection->SetGeometry(poGeometry3);
poLayer->CreateFeature(poFeatureIntersection);
OGRFeature::DestroyFeature(poFeatureIntersection);
OGRFeature::DestroyFeature(poFeature1);
OGRFeature::DestroyFeature(poFeature2);
GDALClose(poSrcDS1);
GDALClose(poSrcDS2);
GDALClose(poDstDS);
return 0;
这里面需要注意的是，获取空间参考应该改为：
OGRSpatialReference *poSpatialRef = poLayer1->GetSpatialRef();
本文结束！
