查东平 ,林 辉 ,孙 华 ,刘年元 ,申 展
(1. 中南林业科技大学 a.林业遥感信息工程研究中心;b.林学院,湖南 长沙 410004;2. 湘西自治州林业局,湖南 吉首 416000)
基于GDAL的遥感影像数据快速读取与显示方法的研究
查东平1a,林 辉1a,孙 华1a,刘年元2,申 展1b
(1. 中南林业科技大学 a.林业遥感信息工程研究中心;b.林学院,湖南 长沙 410004;2. 湘西自治州林业局,湖南 吉首 416000)
利用GDAL开源类库相关函数可以读取多源遥感影像的物理大小、坐标投影系统、传感器波段和影像范围等信息。介绍了多源遥感影像在GDAL框架下的解析模式与方法,从底层开发的关键技术上做了相应的改进和优化,扩展了GDAL读取和快速显示遥感影像的方法。研究结果显示:获取2.13 GB遥感影像像元值,利用RasterIO函数所需要的时间为4.3 s左右,按照论文所用遥感影像按照屏幕缩放窗口大小进行缓存和采样的方法仅需要2.1 s;传统的采样方法无法申请超过2 GB大小的内存空间,论文提出的分块处理技术,能够满足大数据量的应用需求,并且能够提高系统内存应用效率和加快遥感影像缩放过程中的加载速度。
遥感;3S技术;GDAL;Worldview2
显示与处理遥感影像的软件种类繁多,但是从开发的角度上主要可以分为两类[1-2]:一是如ENVI、ERDAS、PCI、ArcGIS、Super Map 等 专业公司开发的软件,这些软件功能强大,安全稳定,但系统要求高,使用成本昂贵,功能复杂,与项目结合程度不高[3];二是软件开发人员对专业软件提供的开放接口进行的二次开发,如ESRI公司提供的Arc Object、Arc Engine[4],二次开发软件紧密结合项目,使用起来比专业软件更加方便,但是也存在不少弊端,如不易扩展、开发与运行时需要安装宿主软件、软件成本高、系统要求高等[5]。开源GIS系统具有自主创新、合作创
新的特点成为国内推动地理空间信息产业发展的一个热点[6]。GDAL凭借其开源、免费、强大的功能、易扩展以及大量热心开发人员长期对其扩展维护等优点成为开源GIS的佼佼者[7]。国内外开发人员借助GDAL类库和相关算法开发了上百个专业软件,例 如:FME,GRASS GIS,IDRISI,Map Server,Open EV,Quantum GIS(QGIS)等[8-9], 甚 至 ESRI的ArcGIS[10]和Google Earth[9,11]等著名软件的很多算法与模块都来源于GDAL开源库。论文在GDAL开源类库的基础上分析了传统的影像金字塔构建方法和重采样方法,研究一个遥感影像数据分块读取和高速缓存、快速显示的方法,增强了遥感影像图像显示效果和提升了计算机内存运用效率。
GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下用于栅格空间数据格式读写的类库[12-13],它使用单一的抽象数据模型,满足多种栅格数据的应用需求。GDAL抽象数据模型如图1所示[14-15]。
图1 抽象数据模型Fig.1 Model of abstract data
GDAL数据模型主要包括:数据集(Dataset)、坐标系统(Coordinate System)、仿射地理转换(Affine Geo Transform)、地理控制点(GCPS)、元数据(Metadata)、波段信息(Raster Band) 和 颜 色 表(Color Table) 等[15-16]。 遥感影像数据集包含了一组影像波段信息和影像属性信息,以及波段的地理参考转换和坐标系统定义。坐标系统采用Open GIS WKT(Well Known Text)描述,包含的内容有:全局坐标系统名称、地理坐标系统名称、大地参照系、参考椭球体信息、投影类型、投影参数等。GDAL数据集采用两种方式描述栅格位置和地理参考之间的关系:仿射变换和大地控制点。元数据用于表达特定数据的辅助信息,大小控制在100 kb内。
目前GDAL的最新版本为1.9.0,它支持超过120种栅格类型数据,下一个版本支持的栅格数据类型将能达到200种[13],几乎囊括了市面上所有的空间数据格式[17-18]。
C#调用C++开发的动态链接库属于非动态链接库,非动态链接库的调用需要使用DllImport属性。DllImport属性指示该属性化方法由非托管动态链接库作为静态入口点公开,并提供对非托管DLL导出的函数进行调用所有的信息。作为最低要求,方法必须提供包含入口的DLL名称。下面的实例说明如何使用DllImport属性调用非托管的DLL:
[DllImport("gdal19.dll", EntryPoint = "GDALO penShared", SetLastError = true)]
public static extern IntPtr GDALOpenShared(byte [] fi leName, GDALAccess access);
打开遥感影像主要调用了GDAL开源类库中的GDALOpenShared方法,通过传入文件路径和数据访问方式能够返回获取遥感影像数据集(dataset)。注意:在获取数据集之前要先进行文件格式驱动注册,GDALAllRegister()函数是针对GDAL支持的所有数据格式的驱动,若读取失败,该函数返回空值。
通过遥感影像数据集能够获得的属性信息包括:数据集名称、大小 、长度、宽度、坐标投影信息和波段信息等。这些信息存放在影像数据集的元数据中[15]。遥感影像的宽度和高度信息需要通过对GDALGetRasterXSize()和GDALGetRasterYSize()函 数 的 调 用,GetGeoTransform函数则能够返回仿射地理坐标系数。仿射变换共由6个参数构成,用式(1)的方法将坐标点映射到地理坐标:
式(1)中:Xgeo、Ygeo为地理坐标(变量);GT(0)~GT(5)为放射变换系数(在一幅影像中是常量);Xpixel为栅格影像上X像素坐标(变量);Yline为栅格影像上列位置(变量)。遥感影像左下角坐标以及分辨率的计算公式分别为:
式(2)中:Xmin、Ymin为影像上左下角坐标;HrasterSize为影像高度;Reswidth为像元宽;Resheight为像元高。
一景遥感影像数据具有若干个波段,从数据集中获取波段数目可以通过调用函数GDALGetBandNumber()。屏幕显示遥感影像需要从其数据集中选取3个波段进行有效组合。从遥感影像中获取一个波段需要通过调用GDALGetRasterBand()函数传入数据集对象和波段索引值。为了便于程序设计实现,做如下设定:如果波段数大于3,读取其前3个波段,按照band1,band2,band3组合方式显示到屏幕。波段另外一个重要属性信息是波段像素范围的最大值和最小值,在程序中可以通过调用GDAL函数GDALGetRasterStatistics()获得。该函数返回值信息中还包括波段像素的平均值和及标准差。 GDALGetRasterBandXSize()和 GDALGetRasterBandYSize() 能分别获取波段以像素为单位的X、Y的最大值和最小值。
遥感影像显示到屏幕之前需将遥感影像波段信息读取存入内存缓冲区[17]。GDAL开源类库中将遥感影像读入内存缓冲区的方法有多种,其中最简单的是调用GDALRasterBand()类中的GDALRasterIO函数。该函数可以实现自动转换数据类型、采样以及影像的裁剪。使用时传入行列偏移量以及x、y方向上的像素大小(nXoff/nYoff/nXsize/nYsize),能够对遥感影像数据全部读取或者分块读取。为了提升遥感影像的读取效率和缩放的平滑性,应该控制目标块在x方向和y方向上的缓冲区大小(nBufXSize/nBufYsize)。当遥感影像文件大于2 GB时,GDAL提供的函数GDALRasterIO()不能将文件全部映射到内存,同时内存映射文件到进程的地址空间时大小有限,可用的空间为2 GB[17],因此需要多遥感影像按照屏幕窗口大小进行分块映射。其步骤如下:
通过缩放视窗的大小重新计算显示在屏幕上的影像像素长度与宽度,其计算公式为:
式(3)中:W为屏幕上影像像素长度;H为屏幕上影像像素宽度;xmax为屏幕上影像x最大值;ymax为屏幕上影像y最大值;Wclient为屏幕宽度;Hclient为屏幕高度。
遥感影像在缩放过程中,需要重新计算遥感影像数据的映射和采样范围,为了使遥感影像缩放过程中不超出影像范围的边界大小,其缩放范围取影像缩放窗口大小和遥感影像大小之间的相交部分的范围(Clip)。窗口x的最小值(Xminclip)取值为缩放窗口与遥感影像边界x的较大值,窗口y的最小值(Yminclip)取缩放窗口和遥感影像边界y的较大值,窗口x的最大值(Xmaxclip)取缩放窗口与遥感影像边界的x的较小值,窗口y的最大值(Ymaxclip)取视图窗口和缩放窗口y的较小值。屏幕上缩放范围的大小和位置通过左下角和右上角两个坐标值来标识:
式(4)中:Xminclip、Yminclip为缩放窗口左下角坐标;Xmaxclip、Ymaxclip为缩放窗口右上角坐标;ptxmin、ptymin为影像对应左下角坐标;ptxmax、ptymax为影像对应右上角坐标。
内存缓冲区的大小和像素空间的大小根据缩放窗口的大小进行确定。在GDAL函数基础上编写的ComputeBufferSizeAndPixelRect()函 数用于实现确定像素空间大小和内存缓冲区的大小。计算大小之前首先判断遥感影像数据在屏幕范围内的数据时候存在金字塔数据,若存在直接从金字塔数据中获取,按照所计算范围进行构建。遥感影像空间像素范围大小和缓冲区大小的计算方式分别如下。
影像像素空间大小:
式(5)中:XmaxrasterExt、YmaxrasterExt为一幅影像右上角坐标;XminasterExt、YminrasterExt为一幅影像左下角坐标;WbandSize为一个波段的宽度;HbandSize为一个波段的高度;WpixelSize、HpixelSize为影像像元大小。转换参数:
式(6)中:dx、dy为x和y方向上转换参数;WmapCtrlPixelSize、HmapCtrlPixelSize为地图控件上像元大小。内存缓冲区大小:
式(7)中:WbufferSize、HbufferSize为内存中对应的宽度和高度;XmaxPixelRect、YmaxPixelRect为缩放影像右上角坐标;XminPixelRect、YminPixelRect为缩放影像左下角坐标。
遥感影像重采样按照BID数据格式颜色模式要求,内存缓冲区中每行数据的字节数必须为4的倍数:
使用Windows API函数memset为新申请的内存做初始化工作。内存初始化之后得到一个指向该内存的指针。根据内存指针使用转换参数对影像重采样之后对内存逐行逐列进行填充。填充之前分别定义一个行列参数用于记录数据库填充的进度。为提高系统运行效率,遥感影像块大小的获取和影像数据范围的计算只需要从数据集中取出一个波段进行计算。波段块的实际组织大小(blockSize)通过GDALGetBlockSize()函 数 获 取。GDALGetRasterBandXSize()和GDALGetRasterBandYSize()能分别获取影像像素宽度和高度。根据影像数据集和波段信息确定缩放视图大小范围确定在影像中的范围,其计算公式为:
式(9)中:xminblock、yminblock为数据块左下角坐标;xmaxblock、ymaxblock为数据块中右上角坐标;WblockSize、HblockSize为数据块宽度和高度。
遥感影像数据块对象在影像中按照行列顺序获取,数据块在影像中的范围获取表达式为:
式(10)中:XminblockRect、YminblockRect为数据块在遥感影像左下角坐标;XmaxblockRect、YmaxblockRect为数据块在遥感影像右上、角坐标;xBlockOff、yBlockOff为行列偏移量。
GDAL提供了GDALReadBlock函数高效读取遥感影像块对象。遥感影像块对象应该与屏幕上的位置——对应,因此数据块存放在内存中后还需要计算所读取的块对象左上角对应在屏幕上的坐标位置。
式(11)中:XclientPixel、YclientPixel为屏幕上对应的坐标。
内存中的影像数据显示在屏幕上,其方式是调用StretchDIBits将DIB影像数据块内像素通过颜色数据拷贝到指定的屏幕目标位置中。若目标矩形比原矩形范围要大,该函数会对数据的行和列进行拉伸,与目标矩形匹配,若目标矩形比原矩形要小,那么该函数会通过指定的光栅操作对其进行压缩。由于本方法严格计算屏幕缩放范围和目标影像范围,影像光栅数据不会进行任何拉伸和压缩。缓冲区中读取的数据逐行列顺序采样显示在屏幕上的正确位置,最后以实参方式返回像素填充位置。
为了说明本算法具有很好的操作性,能够对海量数据支持和快速读取与显示遥感影像,提升系统内存的应用效率,利用本算法和GDAL提供的RasterIO方法进行了对比实验。硬件环境:Dell INSPIRON N4120笔记本;处理器:Intel(R)Core(TM)i5-2450M CPU @ 2.5GHz 2.50GHz;内存:4.00GB;操作系统:64位 Windows 7;开发环境:Visual Studio 2008。影像尺寸:29362*19437;数据大小:2.13 GB;数据格式:TIF;像元数据类型:Byte。
实验结果:获取单个像元值,利用RasterIO()函数代码运行时间为4.3 s左右,采用本文算法分块读取的方式,运行时间为2.1 s左右。遥感影像重采样采用RasterIO()函数对整幅影像进行采样需要申请能容纳整个影像的内存空间,系统通常不允许一次性申请如此大的内存空间[17,19]。研究采用的算法将遥感影像分块处理和采样,对于大于2 GB的遥感影像也能有效读取采样和显示,而且效率也很高,算法运行程序界面如图2所示。
图2 遥感影像显示界面Fig.2 Display interface of remote sensing image
论文在分析了GDAL对遥感影像数据读取的基础上,研究了遥感影像快速显示的算法,主要结论如下:
(1)提出了一种新的基于GDAL快速读取与显示遥感影像的方法。该方法在GDAL开源类库的基础上进行研究和改进,实现了快速构建遥感影像数据集获取、波段信息获取、遥感影像分块采样和数据快速显示等功能。
(2)提升了海量遥感影像的采样和显示效率。按照屏幕显示范围将遥感影像进行分块读取与采样,改进了影像采样技术,提升了遥感影像显示效率,解决了大于2 GB遥感影像的读取与显示问题。
(3)论文提出的遥感影像快速显示方法具有很强的实用性和理论研究意义。利用本方法结合GDAL对遥感影像进行读取和转换,可以脱离商业遥感软件平台运行环境的限制,有利于软件的进一步推广使用。
[1] 刘建国,陆金桂.组件式GIS的矢量图形系统框架研究与实现[J].计算机时代,2004,5:5-7.
[2] 扈 洋.基于GDI+的矢量图形原型系统的设计与实现[D].成都:电子科技大学,2011.
[3] 陈谷良.基于WebGIS的城市电子地图框架设计研究[J].中国西部科技,2010,25(9):32-34.
[4] 苏开君,佃袁勇.流溪河流域数字水温信息管理系统的设计与开发[J].中南林业科技大学学报,2010,30(3):53-65.
[5] 戴吾蛟,邹峥嵘.小型集成地理信息系统建设中的若干问题[J].测绘工程,2002,11(1):18-21.
[6] 郝玉凤.利用开源推动国内地理空间信息产业发展[J].中州大学学报,2010,27(4):118-120.
[7] 刘亚东,李青元.开源库GDAL及其在影像拼接中的应用[J].数字技术与应用,2010,2:88-89.
[8] Open Source Geospatial Foundation .Software Using Gdal –GDAL [EB/OL] (2012-04)[2012-10].http://trac.osgeo.org/gdal/wiki/SoftwareUsingGdal.
[9] 张学宝,包富华.基于开源架构的GIS原型系统的设计与开发研究[J].测绘科学,2010,35(4):210-211.
[10] ArcGIS 9.2 Desktop Help: Supported raster dataset file formats[M]. ESRI. 2007-08-15.
[11] 孟婵媛.基于GDAL和NetCDF的影像金字塔构建方法[J].海洋测绘,2012,2:49-51.
[12] 百度百科. GDAL/OGR [EB/OL].(2012-08)[2012-10]http://baike.baidu.com/view/4087251.htm.
[13] 张 剑,何亚东.基于PDA的遥感影像读取与显示[J].测绘,2010, 33(2):69-71.
[14] 赵辉辉.基于GDAL的农田信息系统研究[D].哈尔滨:东北农业大学,2011.
[15] 刘昌明,陈 荤.GDAL多源空间数据访问中间件[J].地理空间信息,2011,9(5):58-61.
[16] 张宏伟,童恒建.基于GDAL大于2G遥感图像的快速浏览[J].计算机工程与应用,2012,48(13):159-162.
[17] 李 芳,邬群勇.基于GeoRaster的多源遥感数据存储研究[J].测绘科学,2009,4(3):150-151.
[18] 余盼盼,钟志农.基于GDAL的月球空间数据转换服务[J].兵工自动化,2010,29(12):45-48.
[19] 谭庆全,毕建涛.一种灵活高效的遥感影像金字塔构建算法[J].计算机系统应用,2008,4:124-127.
Fast reading and displaying of remote sensing images based on GDAL
ZHA Dong-ping1a, LIN Hui1a, SUN Hua1a, LIU Yuan-nian2, SHEN Zhan1b
(1a. Research Center of Forest Remote Sensing & Information Engineering; b. School of Forestry, Central South University of Forestry &Technology, Changsha 410004 , Hunan , China; 2. Forestry Bureau of Autonomous Prefecture of Xiangxi, Jishou 416000, Hunan, China)
GDAL is a powerful open source library that supports reading and writing of more than 120 types of raster and GIS vector data. By using related functions of GDAL open source library, the physical size of the multi-source remote sensing image coordinate projection information, band information and image range information were obtained. The mode and method to resolve multi-source remote sensing images under GDAL framework were introduced. Some key technologies from the bottom of GDAL methods were improved and optimized and the methods of quickly reading and displaying Geo-images were expanded. The results show that getting a pixel from a remote sensing image of 2.13GB took 4.3 seconds by the RasterIO function and the new method took only 2.1 seconds; The traditional method can’t apply to sample an image larger than 2GB, while the block processing technology method described in this paper can meet the requirement of large-size image, and improved the memory application eff i ciency of the system and also speed up the loading and scaling of the remote sensing image.
remote sensing;3S technology;GDAL;Worldview 2
S771.8
A
1673-923X(2013)01-0058-05
2012-10-16
“十二五”国家高技术研究发展计划(863计划)课题(2012AA102001):“数字化森林资源监测关键技术研究”;林业公益性行业科研专项(201104028):“林分结构与生长模拟技术研究”;国家重大专项项目(E0305/1112/02):“高分湿地资源应用监测示范”;湖南省高校科技成果产业化培育项目(11CY019)
查东平(1985-),男,江西赣州人,硕士研究生,研究方向:林业遥感和地理信息系统
林 辉(1965-),女,湖北黄冈人,教授,博士,博士生导师,主要从事森林经理学、遥感技术与地理信息系统的教学和科研工作
[本文编校:谢荣秀]