OGR1.10中空间叠加函数Union初探

    最近发现GDAL/OGR库发布了1.10版本，新版本中增加了矢量图层的空间叠加分析功能，在OGRLayer类中增加了Union、Intersection、Clip、Erase、Identify等相关函数。简单尝试了其中的Union函数，感觉还是极好的，用起来相对ArcGIS来说更方便。根据自己的理解和OGR提供的帮助，对Union函数的用法简单做一下介绍，其他几个函数的用法跟Union函数应该大同小异，不妥之处，见谅。

    这几个函数的实现首先需要GEOS库的支持，也就是需要编译GDAL的时候打开GEOS库的开关，关于GDAL/GEOS库联合编译的方法网上很多，具体可以参照链接http://blog.csdn.net/wangqinghao/article/details/8279672 。编译好后，便可在OGRLayer的对象中调用Union函数了。先看一下Union函数的原型：Union (OGRLayer *pLayerMethod, OGRLayer *pLayerResult, char **papszOptions=NULL, GDALProgressFunc pfnProgress=NULL, void *pProgressArg=NULL)，根据OGR提供的解释来看，前三个参数是比较重要的：

    pLayerMethod：要跟当前图层做Union的那个图层，如果我要对layer1、layer2做Union操作，通过layer1对象调用Union函数，那么layer2就是pLayerMethod，不能为空。

    pLayerResult：很明显，存储Union结果的图层，也不能为空。这里可以是用ArcGIS创建好的没有要素的图层，也可以是在程序中用OGR创建的新图层，个人倾向与在程序  中用代码创建一个，简单的代码就能实现。

    papszOptions：这个二维指针中有其实包括了四个参数——SKIP_FAILURES=YES/NO、PROMOTE_TO_MULTI=YES/NO、INPUT_PREFIX=string、METHOD_PREFIX=string。当SKIP_FAILURES=YES的时候，如果某个要素的合并出现错误，则会跳过该要素继续后续要素的合并。PROMOTE_TO_MULTI=YES时，可以将单空间对象转换为多空间对象，如将Polygons转换为MultiPolygons， LineStrings转换为MultiLineStrings。INPUT_PREFIX=string，输入图层的属性字段在结果图层中的前置标记。METHOD_PREFIX=string，操作图层的属性字段在结果图层中的前置标记。还是拿layer1和layer2来举例，layer1、layer2合并后的结果图层中包含有layer1和layer2的所有属性字段，那么如果我们设置INPUT_PREFIX=1，METHOD_PREFIX=2，那么在结果图层的字段中，来自layer1的所有字段将在字段名称前加上前缀"1"，来自layer2的将加前缀"2"。

    后面的两个参数没做尝试。

   下面看一个具体的代码示例

  

char *filePath = "D:\\CeShi_Data\\CESHI_NEW";char *layerName1 = "layer1";char *layerName2 = "layer2";OGRLayer *pLayer1 = NULL;OGRLayer *pLayer2 = NULL;OGRDataSource *pODS = NULL;OGRRegisterAll();pODS = OGRSFDriverRegistrar::Open(filePath,TRUE);//读取取要进行Union的两个图层pLayer1 = pODS->GetLayerByName(layerName1);pLayer2 = pODS->GetLayerByName(layerName2);//创建结果图层OGRLayer *pResultLayer = NULL;pResultLayer = pODS->CreateLayer("result",pLayer2->GetSpatialRef(),wkbMultiPolygon,NULL);//配置Union函数中的第三个参数char **p = new char *[4];p[0] = "SKIP_FAILURES=YES";p[1] = "PROMOTE_TO_MULTI=YES";p[2] = "INPUT_PREFIX=1";p[3] = "METHOD_PREFIX=2";pLayer2->Union(pLayer1,pResultLayer,p,NULL,NULL);  //将对pResultLayer的编辑写入文件，如果不加这句,result文件中将没有记录pResultLayer->SyncToDisk();OGRDataSource::DestroyDataSource(pODS);

如下图：

layer1

layer2:

结果图层result:

如果layer1、layer2中的属性字段相同，如示例中，在结果图层中想保存一份属性字段，则在创建结果图层的时给结果图层创建同layer1一样的属性字段集即可，在代码中做如下修改：

//创建结果图层OGRLayer *pResultLayer = NULL;pResultLayer = pODS->CreateLayer("result",pLayer2->GetSpatialRef(),wkbMultiPolygon,NULL);//为结果图层创建跟输入图层layer1相同的属性字段OGRFeatureDefn *pOGRFeatureDefn = pLayer1->GetLayerDefn();for (int i = 0 ; i < pOGRFeatureDefn->GetFieldCount(); i++)pResultLayer->CreateField( pOGRFeatureDefn->GetFieldDefn(i));pLayer2->Union(pLayer1,pResultLayer,p,NULL,NULL);  

生成的结果图层的属性表如下：

    以上便是Union的一些简单用法，其他几个函数也基本相同，不妥之处，请批评指正！