鄒廣黔，夏正清，吳孔江，王靜宇

(貴州省第一測繪院，貴州 貴陽　550025)

?

基于GDAL的標準圖幅生成及數(shù)據(jù)批量裁剪方法*1

鄒廣黔，夏正清，吳孔江，王靜宇

(貴州省第一測繪院，貴州 貴陽550025)

摘要：目前，大多數(shù)GIS軟件均能實現(xiàn)矢量數(shù)據(jù)和柵格數(shù)據(jù)的裁剪，但集標準圖框生成及可定制的批量裁剪為一體的軟件則較少。據(jù)此，文章通過研究矩形圖框、梯形圖框的計算方法以及高斯投影正、反算法，在VC6.0中使用C++及開源GDAL庫實現(xiàn)了各種比例尺標準圖框的批量生成，并可根據(jù)用戶的各類需求實現(xiàn)矢量數(shù)據(jù)(.shp格式)和柵格數(shù)據(jù)(.tif/.img/.pix格式)的批量裁剪，從而實現(xiàn)了GIS數(shù)據(jù)裁剪的高效性及可擴展性。

關鍵詞：GDAL；矢量數(shù)據(jù)；柵格數(shù)據(jù)；批量裁剪；C++；標準圖幅

0引言

現(xiàn)今GIS軟件已能在縣級行政區(qū)范圍甚至省級行政區(qū)范圍內實現(xiàn)所有測繪地理信息數(shù)據(jù)的管理與分析，但對于實際工作來說，將項目區(qū)按一定比例尺劃分是必不可少的：一方面便于圖形的編繪，另一方面便于成果的管理。在測繪地理信息領域，數(shù)據(jù)主要由兩部分構成，即矢量數(shù)據(jù)和柵格數(shù)據(jù)。矢量數(shù)據(jù)是通過記錄坐標的方式將點、線和面等地理空間實體進行表示和存儲，具有“位置明顯、屬性隱含”的特點；柵格數(shù)據(jù)是將地理實體以像元形式表示和存儲，每個像元由行列號確定其具體位置，具有“屬性明顯、位置隱含”的特點[1]。地圖有兩種分幅方式，即矩形分幅和經(jīng)緯線分幅：前者對應比例尺為1∶2 000、1∶1 000及1∶500三種；后者則對應為國家基本比例，分別為1∶100萬、1∶50萬、1∶25萬、1∶10萬、1∶5萬、1∶2.5萬、1∶1萬及1∶5 000八種[2-3]。

以2014年全國范圍內開展的第一次地理國情普查項目為例，該項目要求上交1∶2.5萬或1∶5萬標準圖幅的DOM數(shù)據(jù)、地表覆蓋及地理國情要素數(shù)據(jù)庫等內容[4]。在項目實際開展過程中，貴州省將原始DOM進行1∶1萬標準圖幅裁剪，并以此為基礎進行地表覆蓋、地理國情要素等數(shù)據(jù)的生產，然后再進行接邊得到全縣的地表覆蓋數(shù)據(jù)庫、地理國情要求數(shù)據(jù)庫等內容。在項目后期質量檢查評定過程中，需要按圖幅進行統(tǒng)計評定，為此，需將縣級數(shù)據(jù)庫按照1∶1萬標準圖幅進行裁剪以便統(tǒng)計相關矢量數(shù)據(jù)個數(shù)或面積，且在相關專題資料的使用過程中也要對其進行裁剪，如第一次全國水利普查數(shù)據(jù)。由上述內容可以看出，需要對矢量數(shù)據(jù)(地表覆蓋數(shù)據(jù)、地理國情數(shù)據(jù))及柵格數(shù)據(jù)(DOM)進行標準圖幅裁剪(1∶1萬、1∶2.5萬及1∶5萬)，對于大比例尺測圖，可能會有1∶500、1∶2 000等矢量數(shù)據(jù)及柵格數(shù)據(jù)的裁剪。目前，市面上的GIS軟件大都帶有裁剪功能(如ArcGIS的clip及extract by mask分別用來對矢量數(shù)據(jù)及柵格數(shù)據(jù)進行裁剪)，但基于各類比例尺圖框的批量裁剪則較少。為此，本文利用GDAL開源庫，使用C++在VC6.0中編程實現(xiàn)了基于各類比例尺標準圖框的矢量數(shù)據(jù)和柵格數(shù)據(jù)批量裁剪，且對裁剪的相關要求可實現(xiàn)定制(如外擴距離，左上角坐標等)。

1GDAL概述

GDAL全稱為Geospatial Data Abstraction Library，是一個在X/MIT許可協(xié)議下讀寫柵格數(shù)據(jù)和矢量數(shù)據(jù)的開源庫，它利用抽象數(shù)據(jù)模型來表達所支持的各種格式[5]。GDAL作為GIS領域的開源庫，它提供了對各種格式的柵格數(shù)據(jù)及矢量數(shù)據(jù)的讀寫、轉換、處理等功能，并公布了相關的源代碼，為用戶從底層進行功能擴展提供了便利[6]。GDAL自2007年發(fā)布1.1.0版本以來，至今已經(jīng)到了1.11.1版本，本文使用的版本為2014年4月25日發(fā)布的1.11.0版本。

2數(shù)據(jù)按標準圖幅批量裁剪流程

為便于分幅圖形數(shù)據(jù)的再使用，數(shù)據(jù)按標準圖幅批量裁剪分兩大部分組成，即標準圖幅生成(以圖幅名稱命名的ShapeFile格式面文件)和數(shù)據(jù)裁剪(矢量數(shù)據(jù)及柵格數(shù)據(jù))，具體流程見圖1。

圖1　數(shù)據(jù)批量裁剪流程Fig.1　Batch clipping flow of data

3標準圖幅批量生成

3.1矩形圖框批量生成

本文所述矩形圖框僅限于大比例尺，即1∶500、1∶1 000、1∶2 000，圖幅大小為50 cm×50 cm。矩形圖框以整公里、整百米或整五十米坐標進行分幅，圖幅號以西南角坐標公里數(shù)編號[5-7]。

3.1.1批量計算圖幅號

通過輸入的起止投影直角坐標，批量計算相應比例尺圖幅號，計算方法[7]如下：

1)計算東、北方向的圖幅總數(shù)

row=(MAXY-MINY)/DY；col=(MAXX-MINX)/DX。其中，(MINX，MINY)和(MAXX，MAXY)分別為起止點的坐標，DX和DY為不同比例尺圖幅的東方向和北方向實際距離，見表1。

表1大比例尺DX及DY取值(50 cm×50 cm標準圖幅)

Tab.1Large-scale DX and DY values(50 cm×50 cm standard mapsheet)

比例尺DXDY1∶5002502501∶10005005001∶200010001000

2)獲取所有圖幅的圖幅號

以上一步中求取的row及col進行循環(huán)，得到所有圖幅號，如當前row為i，col為j，則圖幅號計算如下：

x=(MINX+i×DX)/1 000；y=(MINY+j×DY)/1 000；圖幅號即為y-x。

3.1.2批量繪制圖框

在VC中使用GDAL庫提供的相關函數(shù)實現(xiàn)圖框的批量繪制，保存格式為ShapeFile面文件，輸出方式有兩種：一種為一幅一個文件，并以圖幅號作為其文件名；另一種為所有圖框均生成到一個文件中，并為其增加一個TFH字段以便管理?，F(xiàn)以第一種情況為例對批量繪制圖框進行介紹：

1)首先通過圖幅號獲取四角坐標，并存入對應的數(shù)組，對圖框有特殊需求的，可在此計算，如地理國情普查項目中DOM裁剪計算公式[8]為：

(1)

式中：R為DOM分辨率。通過計算得到新的起止坐標，并以此得到四角坐標，以西南角為起點順時針將東坐標及北坐標分別存入X及Y數(shù)組中，通過編寫代碼即可實現(xiàn)圖框文件的批量生成。

2)將數(shù)組中的坐標串轉換為面狀圖形，起止點重合才能形成一個完整面，代碼為：

for(int i=0；i<4；i++)

{ pt.setX(X[i])；pt.setY(Y[i])；pRing->addPoint(&pt)；}

pt.setX(X[0])；pt.setY(Y[0])；

pRing->addPoint(&pt)；

pRing->closeRings()；

pPolygon->addRing(pRing)；

OGRGeometry *pGeo=(OGRGeometry*)pPolygon；

OGRFeature *PF=OGRFeature：：CreateFeature(pDefn)；

PF->SetGeometry(pGeo)；

pLayer->CreateFeature(PF)；

3.2梯形圖框批量生成

梯形圖框主要用于國家基本比例尺，即1∶5 000-1∶50萬，按照最新編號方法對其進行編號。

3.2.1批量計算圖幅號

1)計算起止圖幅號

梯形圖框均以1∶100萬圖框向下擴展而得，其圖幅總數(shù)計算方法較矩形框有所不同，不能通過坐標直接計算，只能通過圖幅號間接計算。通過輸入的起止經(jīng)緯坐標，獲取起止圖幅號，計算公式為：

(2)

式中：a、b表示當前圖幅在1∶100萬比例尺中的行列號，前者用字母表示，后者用實際值表示；c、d表示當前圖幅號的行列號；LONG和LAT表示待求圖幅號的經(jīng)緯度，用十進制度表示；DLONG和DLAT表示不同比例尺的經(jīng)差及緯差(見表2)。

表2中小比例尺經(jīng)差及緯差取值及代碼

Tab.2Latitude and longitude difference values and codes in small and medium scale

比例尺DLONGDLAT總行數(shù)總列數(shù)代碼1∶50001'52.5″1'15″192192H1∶100003'45″2'30″9696G1∶250007'30″5'4848F1∶5000015'10'2424E

對于輸入坐標為投影直角坐標的，需通過高斯反算公式計算其經(jīng)緯度方可計算圖幅號，相關公式見文獻[9]。

2)獲取起止點內的所有圖幅號

依據(jù)起止點圖幅號中的a、b、c、d四個值，批量計算范圍內的所有圖幅號，需考慮如下4種情況進行計算：首先是a、b均不變，即當前比例尺圖幅都在一幅1∶100萬圖幅內；其次是a變b不變；再次是a不變b變；最后是a、b均改變。

3.2.2批量繪制圖框

批量繪制圖框步驟：

1)通過圖幅計算西南角經(jīng)緯度

通過3.2.1的圖幅號計算公式反推計算LAT及LONG，該經(jīng)緯度即為當前圖幅西南角坐標。

2)獲取四角經(jīng)緯坐標

通過第一步中的經(jīng)緯度及不同比例尺的經(jīng)差和緯差，得到四角經(jīng)緯坐標如下：

西南(LONG，LAT)；東南(LONG+DLONG，LAT)；東北(LONG+DLONG，LAT+DLAT)；西北(LONG，LAT+DLAT)。

3)獲取四角直角坐標

通過第二步中獲取的經(jīng)緯度，通過高斯正算公式計算獲取其對應的投影直角坐標，在此可對圖框的特殊要求進行定制，將最終的四角X坐標和四角Y坐標寫入對應的X及Y數(shù)組中，高斯正算公式見文獻[8]。

4)圖形批量生成

按照3.1.2中的方法批量生成ShapeFile格式的面圖框。

4批量裁剪

4.1矢量數(shù)據(jù)批量裁剪

通過讀取被裁剪數(shù)據(jù)文件夾及裁剪框文件夾，實現(xiàn)批量自動裁剪。保存結果有兩種：一種為全部裁剪結果存放到一個文件夾中，并在被裁剪文件名前加入裁剪框名稱，以示區(qū)別；另一種為以裁剪框名稱新建文件夾，將裁剪結果按原名稱保存到對應的文件夾下?，F(xiàn)以第一種方法為例簡要敘述實現(xiàn)過程[5，7]：

1)獲取保存文件名，代碼如下：

EndFileName=DesFol+"\"+CJKShpName+ "-"+BCWJShpName+".shp"；

2)獲取待裁文件，代碼如下：

OGRDataSource *pSDS=OGRSFDriverRegistrar：：Open(BCWJFilePath，F(xiàn)ALSE)；

OGRLayer *pSLayer=pSDS->GetLayer(0)；

3)獲取裁剪框文件，代碼如下：

OGRDataSource *pEDS=OGRSFDriverRegistrar：：Open(CJKFilePath，F(xiàn)ALSE)；

OGRLayer *pELayer=pEDS->GetLayer(0)；

4)創(chuàng)建裁剪結果圖層，代碼如下：

OGRDataSource *pDstDS=pDriver->CreateDataSource(EndFileName，NULL)；

對點、線、面使用不同的參數(shù)進行創(chuàng)建，代碼如下：

OGRSpatialReference *pSRS=pSLayer->GetSpatialRef()；

if(pELayer->GetGeomType()==wkbPolygon)

pDstLayer=pDstDS->CreateLayer("polygon"，pSRS，wkbPolygon，NULL)；

if(pELayer->GetGeomType()==wkbLineString)

pDstLayer=pDstDS->CreateLayer("polygon"，pSRS，wkbLineString，NULL)；

if(pELayer->GetGeomType()==wkbPoint)

pDstLayer=pDstDS->CreateLayer("polygon"，pSRS，wkbPoint，NULL)；

5)執(zhí)行裁剪，代碼如下：

pSLayer->Intersection(pELayer，pDstLayer，NULL，GDALTermProgress，NULL)；

裁剪結果中的屬性名稱既有裁剪框的也保留了被裁數(shù)據(jù)的，若裁剪框只有自帶的FID及Shape兩個字段，則裁剪結果屬性名稱僅為被裁數(shù)據(jù)的，屬性值將復制原有內容，但面積、長度等與幾何圖形相關的屬性在裁剪后并未自動更新，需手動重新計算以更新其值。

4.2柵格數(shù)據(jù)批量裁剪

柵格數(shù)據(jù)裁剪分兩種：基于矩形的規(guī)則裁剪和基于多邊形的不規(guī)則裁剪，GDAL庫中對應的接口為RasterIO函數(shù)和GDALWarp工具。由于本文所述標準圖幅均為規(guī)則裁剪，所以選擇RasterIO函數(shù)來進行影像裁剪。裁剪框坐標系統(tǒng)應與被裁文件保持一致，如果被裁文件為地理坐標，則裁剪框坐標也必須為地理坐標。本文柵格數(shù)據(jù)為投影直角坐標，所以裁剪框坐標必須為投影直角坐標，具體裁剪流程如下：

1)獲取待裁柵格數(shù)據(jù)的6參數(shù)及其它相關信息，代碼如下：

GDALDataset *pSrcDS=(GDALDataset*)GDALOpen(m_RSPATH，GA_ReadOnly)；

double GeoTranform[6]；

pSrcDS->GetGeoTransform(GeoTranform)；

關于6參數(shù)的說明：

0—像素1，1(左上方)的Y坐標；

1—Y方向上的像素分辨率；

2—平移量；

3—像素1，1(左上方)的X坐標；

4—旋轉量；

5—X方向上的像素分辨率。

其中，0、3分別為東方向最小和北方向最大值(即西北方向起始像元中心點坐標)，1、5分別為東方向和北方向分辨率，2、4則為東方向和北方向的旋轉值(一般為0)。

獲取行列值，代碼如下：

int rXSize=pSrcDS->GetRasterXSize()；

int rYSize=pSrcDS->GetRasterYSize()；

獲取波段：

int nDstBC=pSrcDS->GetRasterCount()；

2)獲取西北角直角坐標及東南角直角坐標

對于矩形框，通過圖幅號即可計算西北角直角坐標；對于梯形框，需通過圖幅號計算其西北角經(jīng)緯度后使用高斯正算公式得到直角坐標。此處計算西北角直角坐標為像元起始點的中心點坐標，對于有外擴及其它要求的，在此進行計算即可。

算例：計算1∶1 000比例尺標準圖框外擴100個像元的西北角直角坐標及東南角坐標

MINX=atof(TFH.Right(5))*1000- GeoTranform[1]*100；

MAXY=atof(TFH.Left(6))*1000+500+ ABS(GeoTranform[5])*100；

MAXX=atof(TFH.Right(5))*1000+500+ GeoTranform[1]*100；

MINY=atof(TFH.Left(6))*1000- ABS(GeoTranform[5])*100；

3)獲取在被裁數(shù)據(jù)中的行列號，代碼如下：

double dTemp = GeoTransform[1]*GeoTransform[5] - GeoTransform[2]*GeoTransform[4]；

Col=int((GeoTransform[5]*(X-adfGeoTransform[0])-GeoTransform[2]*(Y-GeoTransform[3]))/dTemp+0.5)；

Row=int((GeoTransform[1]*(Y-GeoTransform[3])-GeoTransform[4]*(X-GeoTransform[0]))/dTemp+0.5)；

4)計算裁剪后的行列數(shù)，代碼如下：

int nPixel=(int)(((MAXX-MINX)/ adfGeoTranform[1])+0.5)；

int nLines=(int)(((MAXY-MINY)/ABS(adfGeoTranform[5]))+0.5)；

5)創(chuàng)建裁剪結果，代碼如下：

GDALDataset *hDstDS=hDriver->Create(pszDstFile，nPixel，nLines，nDstBandCount，eDT，NULL)；

6)執(zhí)行裁剪，代碼如下：

pSrcDS->RasterIO(GF_Read，Col，Row，nPixel，nLines，pDataBuff，nPixel，nLines，eDT，nDstBC，pBand，0，0，0)；

hDstDS->RasterIO(GF_Write，iEndX，iEndY，nPixel，nLines，pDataBuff，nPixel，nLines，eDT，nDstBC，pBand，0，0，0)；

如被裁數(shù)據(jù)小于裁剪框范圍，需正確計算Row、Col、iEndX、iEndY、nPixel、nLines等6個參數(shù)的值才能保證裁剪結果的正確性。分以下4種情況進行設置：

①裁剪框東邊大于被裁文件，代碼如下：

nPixel= rXSize-Col；iEndX=0；iEndY=0；

②裁剪框南邊大于被裁文件，代碼如下：

nLines= rYSize-Row；iEndX=0；iEndY=0；

③裁剪框西邊大于被裁文件，代碼如下：

iEndX=ABS(Col)；nPixel=nPixel-iEndX；iStartX=0；

④裁剪框北邊大于被裁文件，代碼如下：

iEndY=ABS(Row)；nLines=nLines-iEndY；iStartY=0；

5結論

本文通過C++語言及開源GDAL庫，在VC中實現(xiàn)了各種比例尺的標準圖框生成。圖框生成方式有最小外接矩形、標準圖框及最小外接矩形外擴等，結果為一個或若干個ShapeFile面文件；通過圖框可實現(xiàn)矢量數(shù)據(jù)的批量裁剪，通過GDAL庫提供的RasterIO函數(shù)，可實現(xiàn)基于各類比例尺標準圖框的柵格數(shù)據(jù)批量裁剪，對于柵格數(shù)據(jù)裁剪有特殊需求的均可在程序中實現(xiàn)，如計算起始點像元中心點坐標、生成TFW文件、生成柵格數(shù)據(jù)的投影文件等。

[參考文獻]

[1]牟乃夏，劉文寶，王海銀，等.ArcGIS10地理信息系統(tǒng)教程[M].北京：測繪出版社，2012.

[2]祝國瑞，郭禮珍，尹貢白，等.地圖設計與編繪[M].武漢：武漢大學出版社，2008.

[3]劉宏林.國家基本比例尺地形圖新舊圖幅編號變換公式及其應用[J].測繪通報，1998(8)：36-37.

[4]國務院第一次全國地理國情普查領導小組辦公室.第一次全國地理國情普查實施方案[Z].北京：國務院第一次全國地理國情普查領導小組辦公室，2013.

[5]李民錄.GDAL源碼剖析與開發(fā)指南[M].北京：人民郵電出版社，2014.

[6]趙輝輝.基于GDAL的農田信息系統(tǒng)研究[D].哈爾濱：東北農業(yè)大學，2011.

[7]姚領田.精通MFC[M].北京：人民郵電出版社，2007.

[8]第一次全國地理國情普查領導小組辦公室.GDPJ 05—2013 數(shù)字正射影像生產技術規(guī)定[S].北京：第一次全國地理國情普查領導小組辦公室，2013.

[9]孔祥元，郭際明.控制測量學：下冊[M].武漢：武漢大學出版社，2006.

Method of Batch Clipping GIS Data and Generating Standard Mapsheet Based on GDAL

ZOU Guang-qian，XIA Zheng-qing，WU Kong-jiang，WANG Jing-yu

(FirstSurveyingandMappingInstituteofGuizhouProvince，GuiyangGuizhou550025，China)

Abstract：At present，most GIS software can realize the clipping of the vector and raster data，but there has been less softwares which conbines the generation of standard frame and customizable batch clipping.And than this paper studies the calculation method of rectangular frame and a ladder frame，and Gauss projection positive and negative algorithm，using of C++ and the open source GDAL library in VC6.0 to generate the various scale standard frame，and can realize the batch clipping of the vector data(.shp format)and raster data(.tif/.img/.pix format)according to the needs of users，so as to achieve the clipping efficiency and scalability of GIS data.

Key words：GDAL；vector data；raster data；batch clipping；C++；standard mapsheet

* 收稿日期：2016-03-17

中圖分類號：P 208

文獻標識碼：A

文章編號：1007-9394(2016)02-0001-04

作者簡介：鄒廣黔(1974～)，男，貴州平壩人，高級工程師，現(xiàn)主要從事測繪與地理信息應用研究方面的工作。

地礦測繪2016，32(2)：1～4

CN 53-1124/TDISSN 1007-9394

Surveying and Mapping of Geology and Mineral Resources