基于GDAL的标准图幅生成及数据批量裁剪方法*1

邹广黔，夏正清，吴孔江，王静宇

(贵州省第一测绘院，贵州 贵阳　550025)



基于GDAL的标准图幅生成及数据批量裁剪方法*1

邹广黔，夏正清，吴孔江，王静宇

(贵州省第一测绘院，贵州 贵阳550025)

摘要：目前，大多数GIS软件均能实现矢量数据和栅格数据的裁剪，但集标准图框生成及可定制的批量裁剪为一体的软件则较少。据此，文章通过研究矩形图框、梯形图框的计算方法以及高斯投影正、反算法，在VC6.0中使用C++及开源GDAL库实现了各种比例尺标准图框的批量生成，并可根据用户的各类需求实现矢量数据(.shp格式)和栅格数据(.tif/.img/.pix格式)的批量裁剪，从而实现了GIS数据裁剪的高效性及可扩展性。

关键词：GDAL；矢量数据；栅格数据；批量裁剪；C++；标准图幅

0引言

现今GIS软件已能在县级行政区范围甚至省级行政区范围内实现所有测绘地理信息数据的管理与分析，但对于实际工作来说，将项目区按一定比例尺划分是必不可少的：一方面便于图形的编绘，另一方面便于成果的管理。在测绘地理信息领域，数据主要由两部分构成，即矢量数据和栅格数据。矢量数据是通过记录坐标的方式将点、线和面等地理空间实体进行表示和存储，具有“位置明显、属性隐含”的特点；栅格数据是将地理实体以像元形式表示和存储，每个像元由行列号确定其具体位置，具有“属性明显、位置隐含”的特点[1]。地图有两种分幅方式，即矩形分幅和经纬线分幅：前者对应比例尺为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数据、地表覆盖及地理国情要素数据库等内容[4]。在项目实际开展过程中，贵州省将原始DOM进行1∶1万标准图幅裁剪，并以此为基础进行地表覆盖、地理国情要素等数据的生产，然后再进行接边得到全县的地表覆盖数据库、地理国情要求数据库等内容。在项目后期质量检查评定过程中，需要按图幅进行统计评定，为此，需将县级数据库按照1∶1万标准图幅进行裁剪以便统计相关矢量数据个数或面积，且在相关专题资料的使用过程中也要对其进行裁剪，如第一次全国水利普查数据。由上述内容可以看出，需要对矢量数据(地表覆盖数据、地理国情数据)及栅格数据(DOM)进行标准图幅裁剪(1∶1万、1∶2.5万及1∶5万)，对于大比例尺测图，可能会有1∶500、1∶2 000等矢量数据及栅格数据的裁剪。目前，市面上的GIS软件大都带有裁剪功能(如ArcGIS的clip及extract by mask分别用来对矢量数据及栅格数据进行裁剪)，但基于各类比例尺图框的批量裁剪则较少。为此，本文利用GDAL开源库，使用C++在VC6.0中编程实现了基于各类比例尺标准图框的矢量数据和栅格数据批量裁剪，且对裁剪的相关要求可实现定制(如外扩距离，左上角坐标等)。

1GDAL概述

GDAL全称为Geospatial Data Abstraction Library，是一个在X/MIT许可协议下读写栅格数据和矢量数据的开源库，它利用抽象数据模型来表达所支持的各种格式[5]。GDAL作为GIS领域的开源库，它提供了对各种格式的栅格数据及矢量数据的读写、转换、处理等功能，并公布了相关的源代码，为用户从底层进行功能扩展提供了便利[6]。GDAL自2007年发布1.1.0版本以来，至今已经到了1.11.1版本，本文使用的版本为2014年4月25日发布的1.11.0版本。

2数据按标准图幅批量裁剪流程

为便于分幅图形数据的再使用，数据按标准图幅批量裁剪分两大部分组成，即标准图幅生成(以图幅名称命名的ShapeFile格式面文件)和数据裁剪(矢量数据及栅格数据)，具体流程见图1。

图1　数据批量裁剪流程Fig.1　Batch clipping flow of data

3标准图幅批量生成

3.1矩形图框批量生成

本文所述矩形图框仅限于大比例尺，即1∶500、1∶1 000、1∶2 000，图幅大小为50 cm×50 cm。矩形图框以整公里、整百米或整五十米坐标进行分幅，图幅号以西南角坐标公里数编号[5-7]。

3.1.1批量计算图幅号

通过输入的起止投影直角坐标，批量计算相应比例尺图幅号，计算方法[7]如下：

1)计算东、北方向的图幅总数

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进行循环，得到所有图幅号，如当前row为i，col为j，则图幅号计算如下：

x=(MINX+i×DX)/1 000；y=(MINY+j×DY)/1 000；图幅号即为y-x。

3.1.2批量绘制图框

在VC中使用GDAL库提供的相关函数实现图框的批量绘制，保存格式为ShapeFile面文件，输出方式有两种：一种为一幅一个文件，并以图幅号作为其文件名；另一种为所有图框均生成到一个文件中，并为其增加一个TFH字段以便管理。现以第一种情况为例对批量绘制图框进行介绍：

1)首先通过图幅号获取四角坐标，并存入对应的数组，对图框有特殊需求的，可在此计算，如地理国情普查项目中DOM裁剪计算公式[8]为：

(1)

式中：R为DOM分辨率。通过计算得到新的起止坐标，并以此得到四角坐标，以西南角为起点顺时针将东坐标及北坐标分别存入X及Y数组中，通过编写代码即可实现图框文件的批量生成。

2)将数组中的坐标串转换为面状图形，起止点重合才能形成一个完整面，代码为：

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万图框向下扩展而得，其图幅总数计算方法较矩形框有所不同，不能通过坐标直接计算，只能通过图幅号间接计算。通过输入的起止经纬坐标，获取起止图幅号，计算公式为：

(2)

式中：a、b表示当前图幅在1∶100万比例尺中的行列号，前者用字母表示，后者用实际值表示；c、d表示当前图幅号的行列号；LONG和LAT表示待求图幅号的经纬度，用十进制度表示；DLONG和DLAT表示不同比例尺的经差及纬差(见表2)。

表2中小比例尺经差及纬差取值及代码

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

比例尺DLONGDLAT总行数总列数代码1∶50001'52.5″1'15″192192H1∶100003'45″2'30″9696G1∶250007'30″5'4848F1∶5000015'10'2424E

对于输入坐标为投影直角坐标的，需通过高斯反算公式计算其经纬度方可计算图幅号，相关公式见文献[9]。

2)获取起止点内的所有图幅号

依据起止点图幅号中的a、b、c、d四个值，批量计算范围内的所有图幅号，需考虑如下4种情况进行计算：首先是a、b均不变，即当前比例尺图幅都在一幅1∶100万图幅内；其次是a变b不变；再次是a不变b变；最后是a、b均改变。

3.2.2批量绘制图框

批量绘制图框步骤：

1)通过图幅计算西南角经纬度

通过3.2.1的图幅号计算公式反推计算LAT及LONG，该经纬度即为当前图幅西南角坐标。

2)获取四角经纬坐标

通过第一步中的经纬度及不同比例尺的经差和纬差，得到四角经纬坐标如下：

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

3)获取四角直角坐标

通过第二步中获取的经纬度，通过高斯正算公式计算获取其对应的投影直角坐标，在此可对图框的特殊要求进行定制，将最终的四角X坐标和四角Y坐标写入对应的X及Y数组中，高斯正算公式见文献[8]。

4)图形批量生成

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

4批量裁剪

4.1矢量数据批量裁剪

通过读取被裁剪数据文件夹及裁剪框文件夹，实现批量自动裁剪。保存结果有两种：一种为全部裁剪结果存放到一个文件夹中，并在被裁剪文件名前加入裁剪框名称，以示区别；另一种为以裁剪框名称新建文件夹，将裁剪结果按原名称保存到对应的文件夹下。现以第一种方法为例简要叙述实现过程[5，7]：

1)获取保存文件名，代码如下：

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

2)获取待裁文件，代码如下：

OGRDataSource *pSDS=OGRSFDriverRegistrar：：Open(BCWJFilePath，FALSE)；

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

3)获取裁剪框文件，代码如下：

OGRDataSource *pEDS=OGRSFDriverRegistrar：：Open(CJKFilePath，FALSE)；

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

4)创建裁剪结果图层，代码如下：

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

对点、线、面使用不同的参数进行创建，代码如下：

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)执行裁剪，代码如下：

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

裁剪结果中的属性名称既有裁剪框的也保留了被裁数据的，若裁剪框只有自带的FID及Shape两个字段，则裁剪结果属性名称仅为被裁数据的，属性值将复制原有内容，但面积、长度等与几何图形相关的属性在裁剪后并未自动更新，需手动重新计算以更新其值。

4.2栅格数据批量裁剪

栅格数据裁剪分两种：基于矩形的规则裁剪和基于多边形的不规则裁剪，GDAL库中对应的接口为RasterIO函数和GDALWarp工具。由于本文所述标准图幅均为规则裁剪，所以选择RasterIO函数来进行影像裁剪。裁剪框坐标系统应与被裁文件保持一致，如果被裁文件为地理坐标，则裁剪框坐标也必须为地理坐标。本文栅格数据为投影直角坐标，所以裁剪框坐标必须为投影直角坐标，具体裁剪流程如下：

1)获取待裁栅格数据的6参数及其它相关信息，代码如下：

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

double GeoTranform[6]；

pSrcDS->GetGeoTransform(GeoTranform)；

关于6参数的说明：

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)获取西北角直角坐标及东南角直角坐标

对于矩形框，通过图幅号即可计算西北角直角坐标；对于梯形框，需通过图幅号计算其西北角经纬度后使用高斯正算公式得到直角坐标。此处计算西北角直角坐标为像元起始点的中心点坐标，对于有外扩及其它要求的，在此进行计算即可。

算例：计算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)获取在被裁数据中的行列号，代码如下：

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)计算裁剪后的行列数，代码如下：

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

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

5)创建裁剪结果，代码如下：

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

6)执行裁剪，代码如下：

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)；

如被裁数据小于裁剪框范围，需正确计算Row、Col、iEndX、iEndY、nPixel、nLines等6个参数的值才能保证裁剪结果的正确性。分以下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中实现了各种比例尺的标准图框生成。图框生成方式有最小外接矩形、标准图框及最小外接矩形外扩等，结果为一个或若干个ShapeFile面文件；通过图框可实现矢量数据的批量裁剪，通过GDAL库提供的RasterIO函数，可实现基于各类比例尺标准图框的栅格数据批量裁剪，对于栅格数据裁剪有特殊需求的均可在程序中实现，如计算起始点像元中心点坐标、生成TFW文件、生成栅格数据的投影文件等。

[参考文献]

[1]牟乃夏，刘文宝，王海银，等.ArcGIS10地理信息系统教程[M].北京：测绘出版社，2012.

[2]祝国瑞，郭礼珍，尹贡白，等.地图设计与编绘[M].武汉：武汉大学出版社，2008.

[3]刘宏林.国家基本比例尺地形图新旧图幅编号变换公式及其应用[J].测绘通报，1998(8)：36-37.

[4]国务院第一次全国地理国情普查领导小组办公室.第一次全国地理国情普查实施方案[Z].北京：国务院第一次全国地理国情普查领导小组办公室，2013.

[5]李民录.GDAL源码剖析与开发指南[M].北京：人民邮电出版社，2014.

[6]赵辉辉.基于GDAL的农田信息系统研究[D].哈尔滨：东北农业大学，2011.

[7]姚领田.精通MFC[M].北京：人民邮电出版社，2007.

[8]第一次全国地理国情普查领导小组办公室.GDPJ 05—2013 数字正射影像生产技术规定[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～)，男，贵州平坝人，高级工程师，现主要从事测绘与地理信息应用研究方面的工作。

地矿测绘2016，32(2)：1～4

CN 53-1124/TDISSN 1007-9394

Surveying and Mapping of Geology and Mineral Resources