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的一些簡單用法，其他幾個函數也基本相同，不妥之處，請批評指正！