GDAL源码剖析

by giser · 2019-04-05

一、GDAL简介

GDAL官方网站http://www.gdal.org/，本文章中的基本内容都是参照官网中的信息，如有错误或者与官网中的内容冲突，以官网中的为正确。

在开始文章之前，我想先提出几个问题，什么是GDAL？GDAL能做什么？GDAL怎么使用？GDAL内部结构是怎么组织的？GDAL提供的算法原理是什么？对于上面的几个问题，希望大家看完该系列文章后能对上面的几个问题少点疑惑，希望能对感兴趣的童鞋们有所帮助。

什么是GDAL？这个问题比较简单，通俗的讲，GDAL是一个读写空间数据（这里的空间数据包括栅格数据和矢量数据）的开源库（但不仅限于此，此外还提供了一些非常常用的算法和工具）。严格的讲，大家可以参考GDAL首页上的介绍。

GDAL is a translator library for raster geospatial dataformats that isreleased under an X/MIT style Open Source license by the Open Source GeospatialFoundation. As alibrary, it presents a singleabstract data model to the callingapplication for
allsupported formats. It also comes with a variety ofusefulcommandline utilities fordata translation and processing. TheNEWSpagedescribes the July 2011 GDAL/OGR 1.8.1 release.

The related OGR library(which lives within the GDAL source tree) providesa similar capability forsimple features vector data.

二、GDAL目录结构

首先对于GDAL的目录结构进行一个简单的介绍。GDAL源代码下载地址：http://trac.osgeo.org/gdal/wiki/DownloadSource，或者安装svn从源代码服务器下载，svn地址是：http://svn.osgeo.org/gdal/trunk。

如果是使用下载的压缩包，其目录结构如下图：

图1 GDAL源码压缩包文件目录结构

如果使用svn下载的源代码，目录结构如下：

图2 GDAL SVN源码文件目录结构

从上面两张图中可以看出，GDAL的目录结构不管是用什么方式获取的源代码，它的目录结构都是一样的，下面就针对目录结构中的每个文件夹和文件做一个简单的说明。（按照字母顺序来进行说明）

下面先对文件夹进行说明：

1、alg文件夹：alg文件夹中存放的是GDAL库中提供的一些算法的源代码，这些算法包括但不限于：DEM生成等高线算法；图像纠正算法（几何纠正，TPS纠正，正射RPC纠正）；栅格矢量化算法；矢量栅格化算法；格网计算算法；PCT和RGB互转算法；分类图的小碎斑块去除算法等。

2、apps文件夹：apps文件夹中存放的是GDAL库中提供的一些命令行工具集的源代码，这些工具集的介绍可以参考http://gdal.org/gdal_utilities.html，将来我会对这些工具做一个简单的说明。其中有些工具非常的有用，比如gdalinfo，可以使用该工具来查看图像的元数据信息等。

3、bridge文件夹：bridge文件夹中存放的是用来连接GDAL抽象类的定义以及GDAL自己的结构体定义和实现的源代码。具体到后面涉及到GDAL的实现原理时会对该文件夹作一个比较详细的介绍。

4、data文件夹：data文件夹中存放的是GDAL库中需要用到的一些“配置文件”（此处用配置文件可能不太准确），这些文件主要有ESRI的投影文件，ESPG的投影文件，PCI的投影和椭球体文件，autoCAD的头文件，以及其他的一些文件。在GDAL库中有很多时候会自动读取该文件夹中的文件，一般是通过环境变量来查找该文件夹，环境变量的名字叫GDAL_DATA，值就是data文件夹的路径，或者可以在你的程序中使用函数CPLSetConfigOption(“GDAL_DATA”,”C:\GDAL\data”);来进行设置该文件夹的目录，如果没有设置GDAL会自动从环境变量中查找，如果还是没有找到，那么GDAL可能会提示错误，比如如果不设置GDAL_DATA，那么在写如atuocad的dxf格式的时候就提示创建不成功，后面遇到的话会再进行说明。

5、doc文件夹：doc文件夹存放的是用来生产GDAL帮助文档的一些dox文件，dox文件是使用doxygen工具来进行生成的，后面会对doxygen工具作一个简单的介绍和说明，以及在自己的工程中怎么使用doxygen生成自己的程序的开发帮助文档等。总之一句话，这个文件夹就是用来生成GDAL库的帮助文档的一些东西。后面会告诉大家如何生成一份GDAL的帮助文档，当然你也可以把gdal.org整个网站抓下来，J

6、frmts文件夹：这个文件夹可以说是GDAL代码中东西最多的一个文件夹了，每次更新GDAL的版本后都会发现这个文件夹中会多出几个文件夹，同时GDAL也会中支持的文件格式中多出来几个新的文件格式。没错，这个文件夹存放的就是GDAL针对每种不同的特定的图像格式解析的源代码，可以举几个简单的例子，比如bmp文件夹就是解析BMP图像的，hfa文件夹就是用来解析Erdas的img图像格式，raw文件夹用来读取ENVI的hdr文件，还有pcidsk文件夹就是读取PCI的pix格式的等等。所以这个文件夹存放的是解析各个文件格式的源代码。

7、gcore文件夹：通过名字大家也应该知道这个文件夹是做什么的了，叫core的肯定都是很核心的东西了，这个文件夹就是GDAL的灵魂所在，主要存放的GDAL抽象类的数据集，波段，图像读写接口等都在这个里面实现的。如果要知道GDAL的抽象类是怎么对图像格式进行抽象的，可以看看这个里面的代码。

8、html文件夹：html文件夹如果使用压缩包的话，里面应该是空的，这个文件夹主要是用来存放GDAL的生成的帮助文档的地方，主要是使用前面介绍的doc文件夹中的dox脚本，使用doxygen工具生成的GDAL帮助文档会出现在这个文件夹中。后面会和doc文件夹一起进行一个详细介绍。

9、m4文件夹：m4文件夹存放的有好几个已m4为后缀名的文件，m4文件叫MacroProcessor Library，m4文件是编译基础中最核心的文件，这个文件主要是用autoconf来产生configure配置文件，继而自动生成Makefile文件。这个文件夹中Windows平台下貌似没什么作用，可能是我还不知道吧，在此略过。

10、man文件夹：man文件夹貌似是用来生成linux或者其他平台下的帮助文件，估计是可以使用linux下的man帮助吧。Windows平台下貌似也没什么用，略过。

11、ogr文件夹：用过GDAL的肯定知道ogr库吧，在很久很久以前，GDAL和OGR是两个库，GDAL专门负责读取栅格数据，OGR库负责读取矢量数据，然后可能是因为两个库分开有些不方便，比如GDAL的算法库中经常会用到矢量数据的读取，或者还有别的原因吧，现在将这两个库整合在了一起，目前OGR库就是GDAL库的一个子集。其实OGR库还是可以单独编译出来的。Ogr文件夹就是存放OGR库源代码的文件夹。这个文件夹里面也是有很多东西的，后面再详细进行介绍。

12、port文件夹：port文件夹中存放的是port库的东西，port库对于GDAL库来说是一个底层的支持库，port库中定义了一些字符串的操作，文件处理，网页请求，数据库连接，哈希表，字符加密文件压缩等基础的函数。比如GDAL中所有的导出函数符号CPL_DLL就是在这个port文件夹中定义的，还有frmts文件夹中，打开文件，打开数据库，打开网络路径等都是用的port库，以及字符串的处理等。

13、swig文件夹：swig文件夹主要是存放swig的脚本。Swig全称叫SimplifiedWrapper and Interface Generator，网址是www.swig.org， swig的作用就是可以将C/C++写的库封装为Python，C#，Java，Perl和 Ruby等其他语言的访问接口。网上GDAL的C#版本就是使用swig来编译出来实现的，很强大吧。后面有时间的话，会写一篇关于swig编译GDAL的文章。

14、vb6文件夹：这个文件夹中用来将GDAL编译成一个VB6的模块，对于里面具体的文件说明，以及如何编译参考文件夹中的readme.txt，由于我对VB6的使用，还是仅限于拖个按钮，写个单击事件的基础，对于模块之间的调用，实在是不太懂，所在就不说明了。想用VB6使用GDAL的可以自己研究一下，应该也不是很难。

15、wince文件夹：顾名思义，这个文件夹中的东西就是用来编译Windows CE平台下的GDAL库用的，具体编译请参考其中的说明文档。

文件夹介绍完了，下面对文件进行一个大概说明：

1、aclocal.m4：同上面的m4文件夹

2、autogen.sh：Linux平台下的shell文件，用来调用autoconf来产生configure配置文件的。

3、COMMITERS：该文件中的内容是GDAL开发人员的信息，姓名，联系邮箱以及各自负责开发的模块说明等。

4、config.guess，config.sub，configure，configure.in：这四个文件貌似都是linux平台下的配置文件，中Windows下没啥用，略过。

5、Doxyfile：Doxyfile就是前面doc文件夹中说明提到的doxygen的工程文件，用来生成帮助文档用的，后面在介绍doxygen的使用是会对该文件进行一个说明。

6、GDALmake.opt.in：这个文件是linux平台下的GDAL库编译配置文件，功能在后面的nmake.opt中介绍。

7、gdalnightlysvn.sh：Linux平台下调用svn获取GDAL源代码的一个shell脚本。

8、GNUmakefile：GNU的make文件。

9、HOWTO-RELEASE：GDAL发布版本的一些说明。

10、 install-sh：GDAL的安装shell脚本，Linux平台下。

11、 LICENSE.TXT：GDAL的许可说明文件。

12、 ltmain.sh：libtool的shell脚本，Linux平台下，Windows下貌似没用到。

13、 makefile.vc：GDAL的编译文件，用来将源代码编译成dll文件，后面会在GDAL编译中作一个介绍。

14、 makegdal_gen.bat：用来生成VS的工程文件的一个批处理文件，后面在GDAL编译中会对该文件的使用方式做一个说明。

15、 makegdal10.sln：文件夹中所有的sln文件都是VS的项目文件，文件名后的数字代表的是VS的版本号。

16、 makegdal10.vcxproj和makegdal71.vcproj：VS的工程文件，该文件可以由makegdal_gen.bat文件自动生成，后面详细介绍。

17、 mkbindist.sh，mkgdaldist.sh和mktestdist.sh：三个shell脚本文件，Windows下没用，略过。

18、 NEWS：GDAL的新增功能，以及修复的bug记录等。

19、 nmake.opt：GDAL编译选项配置文件，在编译GDAL中，可以指定GDAL绑定的其他库等都在这个里面进行设置。在后面的GDAL编译中会详细介绍说明。

20、 nmake-wince.opt：编译wince版本的编译选项配置文件。

21、 PROVENANCE.TXT：GDAL目录说明文件，如果上面说明的不够清楚，可以参考这个文件。

22、 submake.bat：一个编译的批处理文件，目前没啥用。

23、 svnkeywords.sh：又是svn的一个shell脚本。

24、 VERSION：GDAL版本信息。

一、简单的编译

1、使用VisualStudio IDE编译

首先进入GDAL的源代码目录，可以看到有几个sln为后缀的文件名，比如makegdal10.sln，makegdal80.sln，makegdal71.sln，makegdal90.sln 。这些文件是VisualStudio的工程文件，后面的数字对应的VS的版本号，71表示的VS2003，80表示VS2005，90表示VS2008，还有10表示VS2010等。根据自己电脑安装的VS版本，打开对应的文件，如下图所示（使用VS2008SP1版本，打开makegdal90.sln文件）：

图3 VS2008打开编译GDAL1.8.1

然后在左侧解决方案右键，弹出菜单中选择“生成”或者“重新生成”命令，然后GDAL就会开始编译，等待输出窗口中提示，执行完成，生成成功等信息后，就表示GDAL已经完成编译。同时会在GDAL的源代码目录中会出现gdal.lib，gdal_i.lib，gdal18.dll等文件，如果你没有修改GDAL中的nmake.opt文件的话，那么同时会在你的C盘中会出现一个叫“C:\warmerda\bld”的文件夹，里面会包含三个文件夹，分别是bin，data和html。其中bin文件夹中存放的是编译出来的GDAL的可执行程序，包括GDAL提供的十几个工具集；data文件夹就是在第一节中的介绍的data文件夹；html文件夹中存放的是各种数据格式的说明文档。

2、使用cmd命令行编译

使用cmd命令行编译，首先在“开始菜单\所有程序\Microsoft Visual Studio 2008\Visual Studio Tools\ Visual Studio 2008命令提示”，点击“Visual Studio 2008 命令提示”会弹出下面的界面：

然后使用cd命令，切换到GDAL的源代码目录，如下图所示：

切换到GDAL的源代码目录后，依次敲入下面的命令行后回车，等待编译结束即可。

nmake -f makefile.vc

nmake -f makefile.vc install

nmake -f makefile.vc devinstall

同时还有其他的命令，如：

nmake -f makefile.vc clean

nmake -f makefile.vc MSVC_VER=1400clean

nmake -f makefile.vc MSVC_VER=1400DEBUG=1

上面六行的命令含义依次是：

编译GDAL库

编译GDAL库，并安装（这里安装的意思就是将生成的dll，exe等文件拷贝到C:\warmerda\bld目录），

编译GDAL库，并安装开发者模式（安装的意思同上，开发者模式意思是将开发用的include文件夹中的头文件和lib文件一同拷贝到C:\warmerda\bld目录，此时会在C:\warmerda\bld目录中多出来两个文件夹，分别是include和lib，分别存放的是GDAL的头文件和lib文件，用于调用GDAL库使用）。

清理GDAL库，同时会删除编译GDAL库所生成的临时文件，作用相当于在VS环境中的清理命令。

作用同上，但是添加了一个MSVC_VER=1400，表示使用VS2005编译。

编译GDAL库的debug模式，可以用来调试GDAL源码。

二、自定义编译

GDAL的强大之处不单单在于可以读取栅格和矢量数据，同时它的强大之处还在于下面几个方面，第一可以进行矢量图形之间的一些常用操作，比如：求交，求并，缓冲区等等。第二可以进行投影和坐标转换。如果使用GDAL默认的编译方式，那么上述的两个非常强大的功能您将不能使用，因为GDAL这两大功能是基于另外的两个开源库GEOS（Geometry Engine, Open Source）库和PROJ4库来实现的。下面对这两大库分别做一个简单的说明，以及如何修改编译文件，让GDAL能够拥有这两大功能。

1、集成GEOS

关于GEOS库的说明，网上有很多，同时在GEOS的官网http://geos.osgeo.org有详细的说明，简单的来说，GEOS提供了OGC规范中简单几何要素对象操作的C++语言的实现。在地理信息系统领域，拓扑模型是重要的，其计算方法简单但是很难得以实现。使得GEOS不同于其他项目的也正是“空间谓词”与“空间操作”。空间谓词是比较两个空间对象并返回一个布尔变量值作为结果，它表明了存在于两个空间对象之间特殊的关系。比如典型的空间谓词有Contains(), Intersects(),
Touches(), andCrosses()函数等。GEOS项目中对该些函数的实现是异常强壮的，即使是奇异几何对象或是临时的坐标系统运算也不能使其运算不正常或计算错误。目前绝大多数的商业软件仍然在最基础的空间谓词处理上相对成熟，这正是GEOS项目的重要意义。“空间操作”则主要是对两个几何对象进行计算并且返回一个新的几何实体。比较典型的操作函数如Difference(), Union()以及Buffer()等。GEOS中的操作算法已经被广泛的经过了测试。GEOS类库被各类开源空间信息软件项目广泛应用，使用GEOS，它们可以基于最新的规范的几何实体来完成，同时也拥有了复杂空间方法的实现。

关于GEOS的说明和编译，后面会单独写一篇文件进行介绍，这里假设已经下载的是编译好的GEOS库。

首先使用记事本或者其他的文本编辑器打开GDAL源代码目录下的nmake.opt文件，找到“# Uncomment forGEOS support”这句，大概在405行左右，将下面三行代码：

#GEOS_DIR=C:/warmerda/geos

#GEOS_CFLAGS =-I$(GEOS_DIR)/capi -I$(GEOS_DIR)/source/headers -DHAVE_GEOS

#GEOS_LIB =$(GEOS_DIR)/source/geos_c_i.lib

修改为：

GEOS_DIR=F:\Work\3rdPart\geos-3.2.2

GEOS_CFLAGS =-I$(GEOS_DIR)/capi -I$(GEOS_DIR)/source/headers -DHAVE_GEOS

GEOS_LIB = $(GEOS_DIR)/source/geos_c_i.lib

其中F:\Work\3rdPart\geos-3.2.2是我本机的GEOS存放的主目录，后面两行设置的是GEOS的头文件目录和lib文件路径。设置好后保存即可。对比结果如下图如下：

保存完nmake.opt之后，按照第一步中的编译方式进行编译即可。编译后的GDAL就将会支持图形之间的操作等处理。函数主要是在OGR库中，后面会在OGR库中进行详细的介绍说明。编译后，千万别忘记将geos_c.dll文件拷贝到gdal18.dll的同级目录下，否则会提示你找不到geos_c.dll文件。

2、集成Proj4

Proj4是一套开源的坐标投影转换类库，它可以完成在两套不同制图投影系统之间的转换，同样不同的椭球体或大地基准面之间也可以成功的完成转换。GDAL中用到的坐标转换，投影转换，几何纠正，正射纠正等算法，都离不开坐标转换，也就是说要使用这些算法，必须有proj4库的支持才行。同GEOS库的配置方法，在nmake.opt文件中，找到proj4库的位置，大概在352行左右。将下面的三行代码：

#PROJ_FLAGS =-DPROJ_STATIC

#PROJ_INCLUDE =-Id:\projects\proj.4\src

#PROJ_LIBRARY =d:\projects\proj.4\src\proj_i.lib

修改为：

#PROJ_FLAGS =-DPROJ_STATIC

PROJ_INCLUDE =-IF:\Work\3rdPart\proj-4.7.0\src

PROJ_LIBRARY =F:\Work\3rdPart\proj-4.7.0\src\proj_i.lib

其中第一行表示是否使用静态链接的方式，第二行的路径表示，proj库存放的位置，第三行为proj库的lib文件所在路径。修改后保存即可，对比结果如下图如下：

同GEOS库一样，保存完nmake.opt之后，按照第一步中的编译方式进行编译即可。对于Proj库的使用后面会在有一篇文章对其做一个简单的介绍说明。编译后，同样千万别忘记将proj.dll文件拷贝到gdal18.dll的同级目录下，否则会提示你找不到proj.dll文件。

3、集成HDF数据读取

通过上面GEOS和PROJ库的介绍，相信对gdal的配置文件，nmake.opt有一个比较初步的了解了吧，那么下面对于使用GDAL支持hdf数据的读取也是同样，先下载好hdf4和hdf5的库，我用的是HDF4.2.6和HDF5-1.8.7两个库，在hdf的官方网上有编译好的库，直接下载编译好的库即可，对于hdf库的编译，我没有进行编译过，应该和其他的开源库都是差不多吧。同时官网提供了32位的库和64位的库，这里都是按照32位的库进行介绍，后面会有一个gdal的64位库的编译介绍。

将下载好的HDF4.2.6和HDF5-1.8.7两个库解压，然后修改nmake.opt文件中的278行左右，代码如下：

# Uncomment thefollowing and update to enable NCSA HDF Release 4 support.

#HDF4_PLUGIN = NO

#HDF4_DIR = D:\warmerda\HDF41r5

#HDF4_LIB = /LIBPATH:$(HDF4_DIR)\lib Ws2_32.lib

# Uncomment thefollowing and update to enable NCSA HDF Release 5 support.

#HDF5_PLUGIN = NO

#HDF5_DIR = c:\warmerda\supportlibs\hdf5\5-164-win

#HDF5_LIB = $(HDF5_DIR)\dll\hdf5dll.lib

修改为下面的代码：

# Uncomment thefollowing and update to enable NCSA HDF Release 4 support.

HDF4_PLUGIN = NO

HDF4_DIR = F:\Work\3rdPart\HDF4.2.6_win_x86

HDF4_LIB = $(HDF4_DIR)\dll\hd426m.lib$(HDF4_DIR)\dll\hm426m.lib \

$(HDF4_DIR)\lib\hd426.lib$(HDF4_DIR)\lib\hm426.lib Ws2_32.lib

# Uncomment thefollowing and update to enable NCSA HDF Release 5 support.

HDF5_PLUGIN = NO

HDF5_DIR = F:\Work\3rdPart\HDF5-1.8.7_win_x86

HDF5_LIB = $(HDF5_DIR)\dll\hdf5dll.lib

对比代码如下图：

保存，然后编译gdal即可，同时将hdf库中的dll文件夹下的dll文件拷贝到gdal18.dll的同级目录下。

三、其他方面

1、makegdal_gen.bat使用

对于makegdal_gen.bat的作用，在上一篇文章中已经进行了介绍，下面对怎么使用该文件生成VS的工程文件做一个说明。

首先打开cmd命令行窗口，使用cd命令切换到GDAL源代码目录，然后输入makegdal_gen.bat回车，会得到该工具的一个简单实用帮助，如下图所示：

该工具的使用方法是带有命令行参数的一个批处理工具，（在后面对GDAL工具集的介绍中会对带有参数的命令行程序，以及编写带有命令行的程序有一个比较详细的说明）。通过上图可以看出该工具的基本语法是：

makegdal_gen 7.10 >makegdal71.vcproj

makegdal_gen 8.00 >makegdal80.vcproj

通过上面的示例可以看出，该工具的命令行参数分别是，首先是VS的版本号，具体版本号参考本文第一小节，然后跟一个大于号“>”，最后是输出的VS的工程的名字。那么现在我要使用该命令行生成一个VS2008版本的工程文件，我可以输入下面的命令，然后回车即可：

makegdal_gen 9.00 >makegdal90.vcproj

2、编译64位系统下的GDAL

对于GDAL的64位系统的编译，基本和32位系统的编译一样，首先在VS的工程中，打开配置管理器，然后再活动解决方案平台的下拉列表中选择新建，然后弹出，新建解决方案平台对话框，选择新平台为x64（需要在安装VS的时候安装64位的编译环境），然后点击确定即可。最后在VS中选择X64进行编译即可。如下图所示：

对于使用cmd命令行编译，基本同本文开始，不一样的只有，在开始菜单选择的不是“Visual Studio 2008 命令提示”而是“Visual Studio 2008 x64 兼容工具命令提示”，剩下的编译步骤跟前面的一样。

在编译开始之前，还需要打开nmake.opt文件，找到131行处的“#WIN64=YES”，将前面的“#”去掉，保存，然后开始编译。如果就这样编译过去的话，那么恭喜你，如果不能顺利编译过去的，那么需要按照下面的步骤进行一点点设置。

1：在GDAL目录下的nmake.opt文件中，找到SYM_PREFIX的定义，应该在438行左右

将SYM_PREFIX=_ 改为SYM_PREFIX= 就是将最后的下划线去掉

2：在GDAL目录下的makefile.vc文件中，找到46行左右的代码，如下：

BASE_INCLUDE =/INCLUDE:_GDALSimpleImageWarp@36 \

/INCLUDE:_GDALReprojectImage@48 \

/INCLUDE:_GDALComputeMedianCutPCT@32 \

/INCLUDE:_GDALDitherRGB2PCT@28 \

/INCLUDE:_OCTNewCoordinateTransformation@8$(VB6_SAFEARRAYSYM)

修改为：

BASE_INCLUDE =/INCLUDE:$(SYM_PREFIX)GDALSimpleImageWarp \

/INCLUDE:$(SYM_PREFIX)GDALReprojectImage \

/INCLUDE:$(SYM_PREFIX)GDALComputeMedianCutPCT\

/INCLUDE:$(SYM_PREFIX)GDALDitherRGB2PCT \

/INCLUDE:$(SYM_PREFIX)OCTNewCoordinateTransformation $(VB6_SAFEARRAYSYM)

就是将后面的@开始，后面的数字删除。现在开始编译吧

最近写的项目需要在64位的服务器上，结果32位下编译的不能用，只好重新编译一套64位的。在编译GDAL时，出现了连接错误，如下：

LINK : error LNK2001: unresolved external symbol >_OCTNewCoordinateTransformation at 8

LINK : error LNK2001: unresolved external symbol _vbSafeArrayToPtr at16

LINK : error LNK2001: unresolved external symbol _GDALDitherRGB2PCT at 28

LINK : error LNK2001: unresolved external symbol > >_GDALComputeMedianCutPCT at 32

LINK : error LNK2001: unresolved external symbol _GDALReprojectImage at 48

LINK : error LNK2001: unresolved external symbol _GDALSimpleImageWarp at36

LINK : error LNK2001: unresolved external symbol _OGRRegisterAll

LINK : error LNK2001: unresolved external symbol _OGR_G_GetPointCount

LINK : error LNK2001: unresolved external symbol _OPTGetProjectionMethods

LINK : error LNK2001: unresolved external symbol _OSRValidate

LINK : error LNK2001: unresolved external symbol_OGRFeatureStylePuller

按照GDAL论坛中的方式（http://www.osgeo.org/pipermail/gdal-dev/2006-May/008885.html）不是很清楚。

以下是我编译的步骤，哥研究了很长时间滴：

1：在GDAL目录下的nmake.opt文件中，找到SYM_PREFIX的定义，应该在438行左右

将SYM_PREFIX=_ 改为SYM_PREFIX= 就是将最后的下划线去掉

2：在GDAL目录下的makefile.vc文件中，找到46行左右的代码，如下：

BASE_INCLUDE = /INCLUDE:_GDALSimpleImageWarp@36 \

/INCLUDE:_GDALReprojectImage@48 \

/INCLUDE:_GDALComputeMedianCutPCT@32 \

/INCLUDE:_GDALDitherRGB2PCT@28 \

/INCLUDE:_OCTNewCoordinateTransformation@8 $(VB6_SAFEARRAYSYM)

修改为：

BASE_INCLUDE = /INCLUDE:$(SYM_PREFIX)GDALSimpleImageWarp \

/INCLUDE:$(SYM_PREFIX)GDALReprojectImage \

/INCLUDE:$(SYM_PREFIX)GDALComputeMedianCutPCT \

/INCLUDE:$(SYM_PREFIX)GDALDitherRGB2PCT \

/INCLUDE:$(SYM_PREFIX)OCTNewCoordinateTransformation $(VB6_SAFEARRAYSYM)

就是将后面的@开始，后面的数字删除。

最新版的GDAL1.73版本直接可以在VS2008中选择x64编译，是直接可以编译过去的。前提是需要修改nmake.opt大概第一百行左右中的，Win64 = yes

3、GDAL1.8后打开中文路径失败问题

GDAL1.8.0发布很久了，今天将其更新到1.80发现含有中文路径的文件都不能打开了，影像和矢量文件都是。仔细对比了GDAL1.72和GDAL1.80的代码，终于发现了问题的所在之处，详细代码在GDAL_HOME\port\cpl_vsil_win32.cpp文件中的类VSIWin32FilesystemHandler中，以Stat()函数为例(435行)，其他函数类似。代码如下：

GDAL1.8.0代码（部分）：

[cpp] view plaincopyprint?

1. /************************************************************************/

2. /* Stat() */

3. /************************************************************************/

4. int VSIWin32FilesystemHandler::Stat( constchar * pszFilename,

5. VSIStatBufL * pStatBuf,

6. int nFlags )

7. {

8. (void)nFlags;

10. #if (defined(WIN32) && _MSC_VER>= 1310) || __MSVCRT_VERSION__ >= 0x0601

11. if( CSLTestBoolean(

12. CPLGetConfigOption( “GDAL_FILENAME_IS_UTF8”, “YES” ) ))

13. {

14. int nResult;

15. wchar_t *pwszFilename =

16. CPLRecodeToWChar( pszFilename, CPL_ENC_UTF8, CPL_ENC_UCS2 );

17.

18. nResult = _wstat64( pwszFilename, pStatBuf );

19. CPLFree( pwszFilename );

20.

21. return nResult;

22. }

23. else

24. #endif

25. {

26. return( VSI_STAT64( pszFilename, pStatBuf) );

27. }

28. }

GDAL1.7.2代码（部分）：

[cpp] view plaincopyprint?

1. /************************************************************************/

2. /* Stat() */

3. /************************************************************************/

4. int VSIWin32FilesystemHandler::Stat( constchar * pszFilename,

5. VSIStatBufL * pStatBuf )

6. {

7. return( VSI_STAT64( pszFilename, pStatBuf) );

8. }

通过上面的代码对比，就会看到，原来在函数中添加了一个CPLGetConfigOption( “GDAL_FILENAME_IS_UTF8”, “YES” )判断，通过判断是否是UTF8的编码，而且指定的默认值还是UTF8编码，在含有中文路径的字符串大多数的编码应该是GBK的编码，这样，系统就将GBK的编码当做UTF8的编码来进行转换，结果就是汉字全部是乱码，导致的结果就是找不到文件，所以打不开。

知道原因，那么解决的方式就知道了，大概有下面几种，各有优劣，供大家选择

1：不改变GDAL源代码，在自己调用GDALRegisterAll()和OGRAllRegiser()函数后，加上下面一句即可。

CPLSetConfigOption(“GDAL_FILENAME_IS_UTF8″,”NO”);

这样的优点是，不用改动GDAL的源代码，但是如果自己的工程中经常打开图像的话，每次都要加，比较麻烦。

2：修改GDAL源代码，将下面一句

CPLSetConfigOption(“GDAL_FILENAME_IS_UTF8″,”NO”);

分别添加到GDALAllRegister()函数【GDAL_HOME\frmts\gdalallregister.cpp73行左右】和OGRRegisterAll()函数【GDAL_HOME\ogr\ogrsf_frmts\generic\ogrregisterall.cpp38行左右】中，然后重新编译GDAL即可。这样的方式就和使用以前版本的GDAL一样了，不用改动自己的代码，推荐使用这种方式。

3：修改GDAL源代码，GDAL_HOME\port\cpl_vsil_win32.cpp文件中的全部去掉CPLGetConfigOption全部去掉，或者将后面的YES改为NO，但是该工作量巨大，而且有好多地方，这种方式不推荐。

希望对那些还在为GDAL180中文路径乱码纠结的人们有所帮助。尤其是看到好多人在外面先把中文路径转成utf8的编码，然后再调用GDAL的函数。

一、Swig编译

1、Swig介绍

SWIG全称是Simplified Wrapper and Interface Generator，官方网站：http://www.swig.org/。SWIG是个帮助使用C或者C++编写的软件能与其它各种高级编程语言进行嵌入联接的开发工具。SWIG能应用于各种不同类型的语言包括常用脚本编译语言例如Perl, PHP, Python, Tcl, Ruby and PHP。支持语言列表中也包括非脚本编译语言，例如C#, Common Lisp (CLISP, Allegro CL,
CFFI, UFFI), Java, Modula-3, OCAML以及R，甚至是编译器或者汇编的计划应用（Guile, MzScheme, Chicken）。SWIG普遍应用于创建高级语言解析或汇编程序环境，用户接口，作为一种用来测试C/C++或进行原型设计的工具。SWIG还能够导出XML或Lisp s-expressions格式的解析树。SWIG可以被自由使用，发布，修改用于商业或非商业中。[摘自SWIG官网http://www.swig.org/translations/chinese/index.html]。

下载安装Swig的时候注意下载Swigwin（我下的是swigwin-2.0.4.zip），不要下载源代码，否则不能在windwos下用。下载后解压，将swigwin-2.0.4的解压目录也添加到环境变量Path中去，否则会出现一些些该配置文件的麻烦。检验swig是否成功设置到环境变量Path中的最简单的方式就是在运行中输入swig后回车，如果提示windows找不到swig，那么说明没有设置成功；如果出现一个黑屏一闪而过，那么说明你设置成功了。

2、编译C#版本GDAL

首先，打开nmake.opt文件，找到SWIG=swig.exe这一句，假如没有将swig的目录添加到环境变量中，那么将这句后面的swig.exe修改为swig.exe的全路径，如F:\Work\3rdPart\swigwin-2.0.4\swig.exe。如果设置了环境变量，那么就不需要进行修改了。

然后按照第二篇中的使用cmd命令编译GDAL的方式来进行编译，打开“Visual Studio2008命令提示”并定位到GDAL源代码目录，然后依次执行下面三行命令：

[python] view plaincopyprint?

nmake /f makefile.vc
nmake /f makefile.vc install
nmake /f makefile.vc devinstall

执行完之后会在GDAL_HOME的目录下生成编译好的dll和exe文件以及include、lib等文件夹。接下来，使用cd命令，进入swig\csharp文件夹中：

[python] view plaincopyprint?

cd swig\csharp
#nmake /f makefile.vc interface #这句话暂时不需要
nmake /f makefile.vc
nmake /f makefile.vc install

执行完上面两句后，会在csharp文件夹下生成8个dll文件，并将这8个dll文件拷贝到GDAL输出目录下的csharp文件夹中。然后就可以在C#工程中引用这几个dll了，需要注意的是，在使用这几个dll的时候需要将gdal18.dll以及其依赖的其他dll都要拷贝到同一个目录中才能正常运行。

3、编译Java版本GDAL

首先需要JDK环境，没有的可以下载安装一个。然后打开nmake.opt文件，找到大概86行左右的位置：

[python] view plaincopyprint?

JAVA_HOME = C:\Program Files\Java\jdk1.6.0_26

将后面的路径修改为JDK的目录后保存。在打开的命令行中使用cd命令定位到swig目录后，输入下面的命令即可：

[python] view plaincopyprint?

nmake /f makefile.vc java

4、编译Python版本GDAL

首先在编译Python版本之前，确保自己的电脑中安装了Python，相信大家都装ArcGIS了吧，那么恭喜你，你不用安装Python了，因为在安装ArcGIS的时候必须要安装Python的，同样，将Python的bin目录添加到环境变量Path中去。确保这步完成之后，接下来就是编译Python版本的GDAL。

打开nmake.opt文件，找到大概76行左右的：

[python] view plaincopyprint?

PYDIR = “C:\Python26”

将后面的路径修改为Python的安装路径后保存。使用cd命令将命令行定位到swig\python目录后，依次输入：

[python] view plaincopyprint?

python setup.py build
python setup.py install

执行完上述命令后，会在python的site-packages目录看到多了gdal和ogr的文件以及一个osgeo的文件夹。我本机的目录是C:\Python26\Lib\site-packages。然后将编译出来的gdal18.dll以及它依赖的所有的dll都拷贝到site-packages文件夹中的osgeo文件夹中，最后就可以在Python中使用imort gdal来使用GDAL了。测试Python版本的gdal是否安装成功，可以使用下面命令：

[python] view plaincopyprint?

from osgeo import gdal
from osgeo.gdalconst import *
dataset=gdal.Open(“F:\Work\Data\envi.img”,GA_ReadOnly)
dataset.GetDriver().ShortName
#’ENVI’ #图像格式
dataset.GetProjectionRef()
#’PROJCS[“unnamed”,GEOGCS[“NAD83”,DATUM[“North_American_Datum_1983”,SPHEROID[“GRS 1980”,6378137,298.257222101,AUTHORITY[“EPSG”,”7019″]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[“EPSG”,”6269″]],PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,”8901″]],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,”9108″]],AUTHORITY[“EPSG”,”4269″]],PROJECTION[“Transverse_Mercator”],PARAMETER[“latitude_of_origin”,30],PARAMETER[“central_meridian”,-84.16666666666667],PARAMETER[“scale_factor”,0.9999],PARAMETER[“false_easting”,2296587.926509186],PARAMETER[“false_northing”,0],UNIT[“Foot
(International)”,0.3048]]’ #图像投影信息

注意，上面命令中开头为#号的表示的是Python输出的。如果能输出上面的信息，那么GDAL for Python就安装成功了。

二、开发帮助文档生成

1、Doxygen介绍

以下内容摘自维基百科，地址http://zh.wikipedia.org/wiki/Doxygen。

Doxygen 是一个 C++, C, Java,Objective-C、Python、IDL (CORBA 和 Microsoft flavors)、Fortran、VHDL、PHP、C#和D语言的文档生成器。可以运行在大多数类Unix系统，以及Mac OS X操作系统和Microsoft Windows 。初始版本的Doxygen借鉴了一些老版本DOC++的代码；随后，Doxygen源代码由Dimitri van Heesch重写。

Doxygen是一个编写软件参考文档的工具。该文档是直接写在代码中，因此比较容易保持更新。 Doxygen 可以交叉引用文档和代码，使文件的读者可以很容易地引用实际的代码。

KDE 使用Doxygen作为其部分文档且 KDevelop 具有内置的支持。 Doxygen的发布遵守GNU General Public License，并且是自由软件。

需要注意的是在使用doxygen的时候，会自动生成一些类图，以及函数调用关系图，如果要生成这些图，还需要另外一个很牛叉的开源库Graphviz。Graphviz （Graph Visualization Software的缩写）是一个由AT&T实验室启动的开源工具包，用于绘制DOT语言脚本描述的图形。它也提供了供其它软件使用的库。Graphviz是一个自由软件，其授权为Common Public License。其Mac版本曾经获得2004年的苹果设计奖。更多关于Graphviz的介绍参考其官方网站：http://www.graphviz.org/。

2、生成GDAL开发帮助文档

首先下载安装Doxygen和Graphviz，安装完之后最好将Doxygen和Graphviz的bin目录添加到系统环境变量Path中去。负责要设置一些参数，比较麻烦，还是放到Path中方便。

安装完Doxygen后，会在开始菜单中有个doxywizard.exe的程序，打开，然后在File->Open…菜单中选择GDAL源代码目录下的Doxyfile文件，然后，切换到Run标签，点击Run Doxygen按钮，接下来就会自动提取源代码中的注释生成一份gdal的帮助文档，默认的输出目录是GDAL目录下的html目录。等待生成结束后，点击左下角的Show Html Output后会打开生成的帮助文档。如下图：

再次，基本上GDAL的帮助文档生成完成，打开后会发现和GDAL的官方网站一模一样。但是还有个问题，这里只是生成的是GDAL的帮助文档，没有OGR的帮助文档，同样按照上面的步骤，打开GDAL目录下的OGR文件夹下的Doxyfile。然后点击生成，生成的目录默认为ogr文件夹下有个html文件夹，将该文件夹重命名为ogr，然后整个拷贝到上一层的html中，同时将GDAL目录中的doc文件夹中的文件除dox文件以外的文件全部拷贝到html文件夹中。

最后，还记得编译GDAL后生成的html文件夹吗？将这两个文件夹中的内容进行合并，然后你就得到了一份完完整整的GDAL的帮助文档。首页是html文件夹下的index.html。如下图：

3、Doxygen脚本配置以及代码注释书写规范

对于Doxygen的脚本配置，可以打开doxywizard后，在Wizard和Expert 标签中在每一项中，鼠标移动到上面，在左下角的试图区域会显示很详细的帮助信息（英文的，不是很难），这里有篇说明文档《 Doxygen配置使用指南》，http://read.pudn.com/downloads151/doc/653112/Doxygen_Using_Manual.pdf。文档中包含了所有的东西，很不错的东西

一、GDAL工具通用命令

下面的工具主要参考的GDAL官方网站中提供的帮助文档说明，此外还有我的一些经验，GDAL官方具体地址为：http://gdal.org/gdal_utilities.html。

在所有的GDAL工具集中都会支持下面的通用命令行参数，其形式一般是以两个减号（–）开始，下面详细介绍：

1. –version

输出GDAL的版本信息，即版本号。

2. –formats

输出GDAL支持的所有图像格式说明。包括只读和读写。格式支持描述如下：“ro”是只读驱动；“rw”是读写驱动（比如支持CreateCopy方法）；“rw+”是读写和更新驱动（比如支持Create方法），支持所有的读写更新操作。

3. –format format

输出GDAL单个格式驱动的细节信息。格式名需要是在–formats 后列出所要输出的格式名。比如GTiff，HFA,PCISK等。

4. –optfile file

读取指定名称的文件并把其中的内容当成参数传入命令行列表。如果行首以#开头的行将被忽略。多字组成的参数（即中间有空格隔开的参数）需要用双引号来保正其为单一的参数。

5. –config key value

设置配置，把指定键设置为某个值，从而不必把他们设置为环境变量。一些命令参数键是GDAL_CACHEMAX(用于缓存的内存有多少M)以及GDAL_DATA（gdal的数据路径）。比如在GDAL1.8之后，经常会发现打不开中文路径的文件，那么可以用这个来设置，具体为“—configGDAL_FILENAME_IS_UTF8 NO”。同时对于每一种驱动都会有各自的配置，具体参考各个驱动的说明。更多的配置选项参考该网址：http://trac.osgeo.org/gdal/wiki/ConfigOptions。

6. –debug value

控制调试信息的打印输出。ON值表示允许调试信息输出，OFF值表示不要输出调试信息。

7. –help-general

输出各个工具的命令行参数帮助信息。不同的命令输出的内容不同。

下面还有一些通用的命令，是用来创建文件来使用的。创建不同的格式需要的参数都是不相同的，尤其是在特殊的情况下，比如创建的Erdas的img格式需不需要使用压缩等特殊的需求。这些参数一般使用一个减号（-）开始。下面对这些参数进行一个简单的说明。

8. -of format

选择要创建新的文件的格式。这个格式被指定为类似GTiff（GeoTIFF格式）或者HFA（ERDAS格式）。所有的支持格式列表可以用–formats 参数列出来。但是只有格式列表“(rw)”可以被写入和创建。许多工具如果没有指定，默认是创建GeoTIFF格式的文件。文件扩展名不会自动添加，如果没有指定文件名的后缀名，gdal一般不会添加任何扩展名。各个工具的命令行参数帮助信息。不同的命令输出的内容不同。

9. -co NAME=VALUE

创建文件选项，许多格式会有一个或者更多的创建参数来控制文件创建的细节。比如GeoTIFF或者Erdas的img格式可以用创建参数控制压缩，或者控制是否用分片还是分带来进行存储。

可以使用的创建参数根据格式驱动不同而不同。而一些简单的格式根本就没有创建参数。虽然某个格式可以用”–format<format>”参数列出所有可用的参数列表，但是更详细的信息可以在格式介绍网页中查到。对于不同的文件格式，请参考对应文件格式说明网页。

10. -a_srs SRS

指定输出文件的投影信息（坐标系统）。输出各个工具的命令行参数帮助信息。不同的命令输出的内容不同。有几个工具(如gdal_translate、gdalwarp)可以在命令行中通过类似-a_srs(分配输出SRS)、-s_srs(源SRS)、-t_srs(目标SRS)来指定坐标系统。这些工具允许以一系列格式定义坐标系统（SRS就是空间参考系统spatial referencesystem）。SRS通常可以使用下面几种方式来指定：

NAD27/NAD83/WGS84/WGS72:这些常见的地理坐标系统可以通过名字来直接使用。
EPSG:n：坐标系统（投影或者地理坐标）可以通过EPSG码来选择。例如EPSG 27700是英国国家网格。更多的EPSG坐标系统可以在GDAL数据文件gcs.csv和pcs.csv中找到（位于GDAL目录中的data文件夹中）。
PROJ.4定义：一个PROJ4定义字符串可以用作坐标系统定义。例如“+proj=utm +zone=11 +datum=WGS84”。注意在命令行中要保持Proj4字符串在一起作为一个单独的参数（一般用双引号引起来）。
OpenGIS WKT字符串: OpenGIS标准定义了一个文本格式来描述坐标系统作为简单要素规范的一个部分。这个格式是gdal中使用的坐标系统的内部工作格式。包含wkt坐标系统描述的文件的文件名可以被用来作为坐标系统参数，或者坐标系统元素本身也可以被用来作为命令行参数。
ESRI WKT字符串：ESRI 在他们的ArcGIS产品（ArcGIS中的.prj文件）中使用了一种经过精简OGC WKT的格式，而且这个格式被用在一个和wkt相似的风格的文件中。但是文件名要被加以ESRI::前缀。比如”ESRI::NAD 1927 StatePlane Wyoming West FIPS 4904.prj”。
空间参考网址URLs：可以使用一个空间参考的网址来指定，如：http://spatialreference.org/ref/user/north-pacific-albers-conic-equal-area/。
文件名：可以使用一个包含WKT、Proj.4的字符串，或者XML/GML格式的坐标系统定义的文件。

二、GDAL工具说明

1. gdalinfo 输出文件信息

用法：

gdalinfo[–help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]

[-noct] [-nofl][-checksum] [-proj4][-mdd domain]*

[-sd subdataset] datasetname

参数说明：

gdalinfo程序输出gdal支持的栅格格式的一系列信息。

-mm

强制计算栅格每个波段的最大最小值。

-stats

读取和现实图像统计信息，如果指定该参数，将强制计算图像的统计信息，如各个波段的最大值、最小值、均值、标准差等。

-hist

输出所有波段的直方图信息。

-nogcp

禁止地面控制点（GCP）列表打印。这可能对大量的GCP的数据集来说是十分有用的。比如L1B AVHRR或者hdf4MODIS数据，这些数据包含了成千上万的地面控制点。

-nomd

禁止元数据打印，一些数据集可能包含极多的元数据字符串。

-noct

禁止输出颜色表。

-checksum

强制计算数据集中所有波段的checksum。

-mdd domain

输出指定区域的元数据信息。

-nofl

仅显示文件列表中的第一个文件信息。GDAL1.9.0开始支持该参数。

-sd subdataset

如果输入的数据集包含几个子数据集，那么将使用指定的数字来替代（从1开始）完整的子数据集名称。GDAL1.9.0开始支持该参数。

-proj4

输出文件的坐标系统按照PROJ.4类型的字符串输出。GDAL1.9.0开始支持该参数。

gdalinfo同时会输出如下的信息（如果有的话）：

当前文件的格式驱动信息
栅格数据大小(行列数)
文件的坐标系统(OGC WKT形式)
图像关联到地理的转换参数(当前不包含旋转系数)
地理上的边界坐标，如果可能的话还有基于经纬度的完整的地理转换参数（如果是GCPs就没有）
地面控制点(GCPs)
所有的（包括子栅格的元数据）文件元数据
波段数据类型
波段颜色信息（RGB,Gray等）
波段颜色表信息
波段瓦片大小（文件块大小）
波段描述
波段最大最小值（已经经过计算的）
波段CheckSum值（已经经过计算的）
波段无意义值（NODATA值）
波段可获得的略缩图分辨率
波段单位类型(如：波段的高程是米制还是英制)
波段的假颜色列表

举例：

gdalinfoF:/Work/Data/utm.tif

Driver:GTiff/GeoTIFF

Sizeis 512, 512

CoordinateSystem is:

PROJCS[“NAD27/ UTM zone 11N”,

GEOGCS[“NAD27”,

DATUM[“North_American_Datum_1927”,

SPHEROID[“Clarke1866”,6378206.4,294.978698213901]],

PRIMEM[“Greenwich”,0],

UNIT[“degree”,0.0174532925199433]],

PROJECTION[“Transverse_Mercator”],

PARAMETER[“latitude_of_origin”,0],

PARAMETER[“central_meridian”,-117],

PARAMETER[“scale_factor”,0.9996],

PARAMETER[“false_easting”,500000],

PARAMETER[“false_northing”,0],

UNIT[“metre”,1]]

Origin= (440720.000000,3751320.000000)

PixelSize = (60.000000,-60.000000)

CornerCoordinates:

UpperLeft ( 440720.000, 3751320.000) (117d38’28.21″W,33d54’8.47″N)

LowerLeft ( 440720.000, 3720600.000) (117d38’20.79″W,33d37’31.04″N)

UpperRight ( 471440.000, 3751320.000)(117d18’32.07″W,33d54’13.08″N)

LowerRight ( 471440.000, 3720600.000)(117d18’28.50″W,33d37’35.61″N)

Center ( 456080.000, 3735960.000)(117d28’27.39″W, 33d45’52.46″N)

Band1 Block=512×16 Type=Byte, ColorInterp=Gray

2. gdal_translate 格式转换

用法：

gdal_translate[–help-general]

[-ot{Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}][-strict]

[-of format] [-b band] [-mask band][-expand{gray|rgb|rgba}]

[-outsize xsize[%] ysize[%]]

[-unscale] [-scale [src_minsrc_max[dst_min dst_max]]]

[-srcwin xoff yoff xsizeysize][-projwin ulx uly lrx lry]

[-a_srs srs_def] [-a_ullr ulx ulylrxlry] [-a_nodata value]

[-gcp pixel line eastingnorthing[elevation]]*

[-mo “META-TAG=VALUE”]*[-q][-sds]

[-co “NAME=VALUE”]*[-stats]

src_dataset dst_dataset

参数说明：

gdal_translate工具可以用来在不同格式间转换栅格数据。同时还可以做一些诸如提取子栅格、重采样和数据类型转换等操作。

-ot: type

指定输出波段的数据类型。

-strict:

对于这转换过程中出现丢失数据等错误直接报错，不进行忽略处理。之前的参数叫-not_strict，表示对其进行宽大处理，不报错。

-of format:

选择输出格式。默认是GeoTiff（GTiff）。注意，指定的时候用格式简称。

-b band:

选择一个波段来输出。波段号从1开始编号。多个 -b 参数可以用于选择输出某几个波段或者重新对波段进行排序。从GDAL1.8.0开始，波段可以设置为“mask,1”（或者直接mask）来讲输入数据集的第一个波段作为一个mask波段来使用。

-mask band:

（从GDAL1.8.0开始）(GDAL >= 1.8.0)选择一个输入波段来作为创建输出数据的掩码波段。波段数是从1开始，band可以设置为none，用来避免复制整个输入数据集作为掩码，否则在默认情况下（“auto”），除非掩码是一个alpha通道，或者使用参数-b mask来进行输出，参数band也可以设置为“mask,1”（或者直接mask）来讲输入数据集的第一个波段作为一个mask波段来使用。

-expand gray|rgb|rgba:

（从GDAL1.6.0开始）将带有颜色表信息的单波段文件展开为RGB三波段文件或者RGBA四波段文件。对于输出格式为JPEG，JPEG2000，MrSID，ECW等不支持颜色表的数据来说很有用。灰度值（从GDAL1.7.0开始）也可以使用颜色表展开为一个数据集，但输出文件中仅仅包含灰度级别的一个索引。

-outsize xsize[%] ysize[%]:

设置输出文件的大小。大小以象元为单位，除非用“%”来标记，这时，表示的是输出为输入图层大小的百分比。

-scale [src_min src_max[dst_mindst_max]]:

重新组织输入象元的值。将它们从src_min~src_max范围内缩放到dst_min ~ dst_max范围内。如果省略，输出范围将为0~255。输入范围将由源数据自动计算。

-unscale:

大概意思是，不对波段中的数据进行缩放转换，在使用-ot设置输出文件类型时这个参数往往是很有用的。（这个参数的英文有点绕口，不好翻译，有知道的同学麻烦告诉我一下）。

-srcwin xoff yoff xsize ysize:

选择一个取值窗口，通过该窗口在原图像中的行列位置来拷贝数值。

-projwin ulx uly lrx lry:

选择一个取值窗口，通过该窗口在原图像中地理坐标范围来拷贝数据（类似srcwin）。参数中的四个值，使用的投影坐标。

-a_srs srs_def:

给输出文件投影强制指定坐标系。srs_def可以是任何常用的GDAL/OGR格式的投影信息，如：WKT、Proj4、EPSG:n 或者一个包含着wkt的文件的文件名。

-a_ullr ulx uly lrx lry:

强制指定输出文件的空间转换边界范围（图像的四至范围）。而将原图像的四至范围忽略掉。

-a_nodata value:

指定一个无意义值到输出波段。从GDAl1.8.0开始，可以设置为none来使用原文件中的nodata值作为输出文件的nodata值。

-mo “META-TAG=VALUE”:

如果可以，给输出数据设置一个元数据的键和其对应的值。

-co “NAME=VALUE”:

通过一个创建参数来指定输出格式特殊创建要求。多个-co 参数可以组合起来使用。创建参数可以参考各个数据格式本身说明。

-gcp pixel line eastingnorthingelevation:

添加指定地面控制点到输出数据集。这个选项可以被多次使用，以提供一系列的地面控制点GCPs 。GCP格式为：列号行号横坐标纵坐标高程值。

-q:

安静模式，不输出进度信息以及其他非错误信息。

-sds:

拷贝文件中所有子数据集到各自的输出文件中。通常这个参数用在HDF或者OGDI这样有子数据集的格式中。

-stats:

强制计算（重新计算）数据的统计信息。自GDAL1.8.0开始支持该参数。

src_dataset:

输入数据集名称，可以是文件名，或者是一个多数据集文件中的一个子数据集的URL地址（比如HDF数据集中的一个子数据集）。

dst_dataset:

输出文件名。

举例：

将utm.tif转换为一个以分块存储的GeoTiff文件。

gdal_translate -of GTiff-co “TILED=YES” utm.tif utm_tiled.tif

创建一个JPEG压缩的Tiff图像，同时使用内部掩码从一个RGBA数据集中。

gdal_translate rgba.tifwithmask.tif -b 1 -b 2 -b 3 -mask 4 -coCOMPRESS=JPEG -co PHOTOMETRIC=YCBCR–config GDAL_TIFF_INTERNAL_MASK YES

创建一个RGBA图像从一个RGB数据中使用一个掩码。

gdal_translatewithmask.tif rgba.tif -b 1 -b 2 -b 3 -b mask

3. gdaladdo 建立金字塔

用法：

gdaladdo[-r {nearest,average,gauss,cubic,average_mp,average_magphase,mode}]

[-ro] [-clean][–help-general]filename levels

参数说明：

gdaladdo工具可以用于为大多数支持的格式建立或者重建金字塔。可以使用下面几种重采样算法中的一种来进行缩小重采样操作。

-r {nearest(default),average,gauss,cubic,average_mp,average_magphase,mode}:

指定重采样方法，分别是最邻近（默认），均值，高斯，立方卷积等。

-ro:

（从GDAL1.6.0开始有效），以只读的方式打开数据，来创建外部金字塔（特别对于GeoTIFF来说）。

-clean:

删除所有的金字塔。（从GDAL1.7.0开始有效）。

filename:

要建立金字塔的文件名。

levels:

要建立略缩图的层号的列表。如果使用-clean选项是将被忽略。

模式，（从GDAL1.6.0开始有效）选择最常用的重采样方式。average_mp 是不适合使用的， Average_magphase用于复数数据空间的图像。Nearest 和average 用于普通的图像。Nearest 使用最邻近采样（简单采样），它计算所有的有效值的均值来进行计算。Cubic 采样（从GDAL1.7.0开始有效）使用一个4×4的近似立方卷积内核。 Gauss 采样（从GDAL1.6.0开始有效）使用高斯内核计算。这种对于高对比度和图案边界比较明显的图像效果比较好。一般建议的采样比值是2,4,8,…，使用3×3重采样作为高斯采样的计算窗口。

gdaladdo将遵守正确NODATA_VALUES元组（特殊的数据集元数据），因此，只有一个给定的的RGB三元组（在一个RGB图像的情况下）作为NODATA值，而不是每个波段有独立的NODATA值。

选择一个缩放级别值如2表示略缩图缩放程度是源图像每个维上分辨率的1/2。如果文件在所选缩放级别上已经存在略缩图，那么这个缩放级别上的缩略图将被重新计算并覆盖写入。

一些格式根本不支持金字塔。所以许多格式在文件以外以扩展名.ovr存储金字塔，TIFF就是如此。GeoTIFF格式直接把金字塔存储到原有的文件中。除非使用-ro标记来指定。在TIFF中创建金字塔可以通过用COMPRESS_OVERVIEW配置参数进行压缩。所有GeoTIFF支持的压缩方法，可以在这里获得（如：–config COMPRESS_OVERVIEW DEFLATE）。

大多数驱动也支持一个备用的略缩图格式（使用的是Erdas图像格式）。要使用这个备用格式使用USE_RRD=YES 来设置参数。这样做会把GDAL程序创建的金字塔放到一个辅助的.aux文件中使得可以该金字塔可以直接在Erdas中使用或者也可以在ArcGIS中使用。关于如何使用GDAL创建Erdas格式的金字塔，请参考我的博文：http://blog.csdn.net/liminlu0314/article/details/6127755。

举例：

在所指定的TIFF文件内部创建金字塔：

gdaladdo -r average abc.tif 2 4 8 16

从一个ERDAS.IMG文件中创建一个外部的压缩的金字塔文件：

gdaladdo –config COMPRESS_OVERVIEWDEFLATE erdas.img 2 4 8 16

为给定JPEG文件创建一个Erdas Imagine 格式金字塔：

gdaladdo –config USE_RRD YES airphoto.jpg3 9 27 81

其中有个nearblack,gdalbuildvrt工具，没有做说明，以后加上，里面可能有很多不足，望大家批评指正，谢谢。

1. gdalwarp 图像纠正

用法：

gdalwarp [–help-general] [–formats]

[-s_srs srs_def] [-t_srs srs_def] [-to”NAME=VALUE”]

[-order n | -tps | -rpc | -geoloc] [-et err_threshold]

[-refine_gcps tolerance [minimum_gcps]]

[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-tswidth height]

[-wo “NAME=VALUE”] [-ot Byte/Int16/…] [-wtByte/Int16]

[-srcnodata “value [value…]”] [-dstnodata”value [value…]”] -dstalpha

[-r resampling_method] [-wm memory_in_mb] [-multi] [-q]

[-cutline datasource] [-cl layer] [-cwhere expression]

[-csql statement] [-cblend dist_in_pixels][-crop_to_cutline]

[-of format] [-co “NAME=VALUE”]* [-overwrite]

srcfile* dstfile

参数说明：

gdalwarp工具是一个图像镶嵌、重投影、和纠正的工具。程序可以重投影到任何支持的投影。

-s_srs srs def:

设置原始空间参考，坐标系统是可以使用函数OGRSpatialReference.SetFromUserInput()调用的就行，包括EPSG PCS，PROJ4或者后缀名为.prf的wkt文本文件。

-t_srs srs_def:

设置目标空间参考，坐标系统是可以使用函数OGRSpatialReference.SetFromUserInput()调用的就行，包括EPSG PCS，PROJ4或者后缀名为.prf的wkt文本文件。

-to NAME=VALUE:

设置转换参数选项，具体选项参考函数GDALCreateGenImgProjTransformer2()支持的选项。

-order n:

多项式纠正次数（1到3），默认的多项式次数根据输入的GCP点个数自动计算。

-tps:

强制使用TPS（thin plate spline）纠正方法来纠正图像。

-rpc:

强制使用RPC参数纠正。

-geoloc:

强制使用Geolocation数组。（这个没用过，不太清楚）

-et err_threshold:

指定变换的近似误差阈值，默认为0.125个像元大小（使用像元为单位）。

-refine_gcps tolerance minimum_gcps:

(GDAL >= 1.9.0) refines the GCPs by automatically eliminating outliers.Outliers will be eliminated until minimum_gcps are left or when no outliers canbe detected. The tolerance is passed to adjust when a GCP will be eliminated.Not that GCP refinement
only works with polynomial interpolation. The toleranceis in pixel units if no projection is available, otherwise it is in SRS units.If minimum_gcps is not provided, the minimum GCPs according to the polynomialmodel is used. （这个参数不太清楚，没用过）

-te xmin ymin xmax ymax:

设置输出文件的地理范围（在目标空间参考中）。

-tr xres yres:

设置输出图像的分辨率。

-tap:

(GDAL >= 1.8.0) (target aligned pixels) align the coordinates of theextent of the output file to the values of the -tr, such that the alignedextent includes the minimum extent.

-ts width height:

设置输出文件的宽高。如果宽或者高有一个为0，那么将自动计算一个值，注意-ts和-tr不能同时使用。

-wo “NAME=VALUE”:

设置纠正选项。具体参考GDALWarpOptions::papszWarpOptions的帮助文档。

-ot type:

指定输出波段的数据类型。

-wt type:

计算的数据类型。

-r resampling_method:

重采样方式，主要有下面几种方式：

near:

最邻近采样方法（默认值，算法较快，但是质量较差）。

bilinear:

双线性内插采样。

cubic:

立方卷积采样。

cubicspline:

立方样条采样。

lanczos:

Lanczos 窗口辛克采样。

-srcnodata value [value…]:

设置输入波段的Nodata值，可以为不同的波段指定不同的值。。如果有多个值，就需要把他们用双引号括起来，以保持在命令参数中作为单一参数输入。掩膜值不会在内插中处理。

-dstnodata value [value…]:

设置输出波段的Nodata值。

-dstalpha:

创建一个Alpha波段在输出文件中。

-wm memory_in_mb:

设置纠正API使用的内存大小，以MB为单位。

-multi:

是否使用多线程纠正图像，多线程用来分块处理，同时在读取和写入图像均使用多线程技术。

-q:

不在控制台输出提示信息。

-of format:

输出文件格式，默认为GeoTiff。

-co “NAME=VALUE”:

指定创建图像选项，具体参考不同的格式说明。

-cutline datasource:

使用使用OGR支持的矢量数据进行裁切图像。

-cl layername:

指定裁切矢量的图层名称。

-cwhere expression:

从裁切矢量中根据属性表查询指定的要素来裁切图像。

-csql query:

使用SQL语句来从裁切矢量的属性表中查询要素来裁切图像。

-cblend distance:

Set a blend distance to use to blend over cutlines (in pixels).（这个参数不太清楚，没用过）

-crop_to_cutline:

(GDAL >= 1.8.0) 使用矢量边界的外接矩形大小作为输出影像的范围。

-overwrite:

(GDAL >= 1.8.0) 如果结果数据存在，那么覆盖结果数据。

srcfile:

输入数据文件名（可以为多个，使用空格隔开）。

dstfile:

输出数据文件名。

如果输出文件已经存在，那么镶嵌到这个文件是可以的。但是数据的空间范围等信息不会被修改，如果要修改为新的数据的空间信息，那么需要使用-overwrite选项来覆盖原文件。

举例：

在此不再列举例子，后面会对gdalwarp做一个比较详细的介绍。

2. gdaltindex 创建一个栅格文件块索引Shp文件

用法：

gdaltindex [-tileindex field_name] [-write_absolute_path][-skip_different_projection] index_file [gdal_file]*

参数说明：

这个程序创建一个shape文件，里面记录每一个输入的栅格文件的文件路径，以及栅格数据的落图文件。这个输出文件主要是为MapServer 来使用的。应该是创建栅格图像的四至范围矢量文件。

· 如果shape文件不存在会自己创建，如果存在，将被追加到存在的文件中。

· 默认的索引字段为location。

· 默认填写的是栅格文件名，如果指定了-write_absolute_path选项，那么写入的是栅格文件的绝对路径。

· 如果指定-skip_different_projection参数，只用相同的投影的栅格文件会被插入到tileindex文件中。

· 生成一个简单的矩形多边形，其坐标系统和栅格的坐标系统相同。

举例：

gdaltindex doq_index.shp doq/*.tif

3. gdalbuildvrt 创建VRT数据

用法：

gdalbuildvrt [-tileindex field_name] [-resolution{highest|lowest|average|user}]

[-tr xres yres] [-tap] [-separate] [-allow_projection_difference] [-q]

[-te xmin ymin xmax ymax] [-addalpha] [-hidenodata]

[-srcnodata “value [value…]”] [-vrtnodata “value[value…]”]

[-input_file_list my_liste.txt] [-overwrite] output.vrt [gdalfile]*

参数说明：

从指定的输入栅格文件中创建一个VRT数据。从gdal1.6开始支持。

该工具没有进行说明。

4. gdal_contour DEM生成等高线

用法：

gdal_contour [-b <band>] [-a <attribute_name>] [-3d][-inodata]

[-snodata n] [-f <formatname>] [-i <interval>]

[-off <offset>] [-fl <level> <level>…]

[-nln <outlayername>]

<src_filename> <dst_filename>

参数说明：

使用DEM生成等高线矢量数据。

-b band:

指定要处理的波段序号，默认为第一个波段。

-a name:

指定一个属性表的字段名称用于保存高程值，如果不指定则不会在属性表中填写。

-3d:

强制使用3D定点而不是2D，即高程值将被写入等高线矢量中。

-inodata:

忽略任何nodata值。

-snodata value:

设置输入文件的nodata值。

-f format:

创建输出矢量的格式，默认为ESRI的shape文件。

-i interval:

等高线间距。

-off offset:

设置高程相对系数。

-fl level:

Name one or more “fixed levels” to extract.

-nln outlayername:

指定输出矢量文件的图层名称，默认为“contour”。

举例：

使用dem.tif生成一个等高距为10米的等高线矢量文件contour.shp。

gdal_contour -a elev dem.tif contour.shp -i 10.0

5. gdaldem DEM算法

用法：

– To generate a shaded relief map from any GDAL-supported elevation raster:

gdaldem hillshade input_dem output_hillshade

[-z ZFactor (default=1)] [-s scale* (default=1)]”

[-az Azimuth (default=315)] [-alt Altitude (default=45)]

[-alg ZevenbergenThorne]

[-compute_edges] [-b Band (default=1)] [-of format] [-co”NAME=VALUE”]* [-q]

– To generate a slope map from any GDAL-supported elevation raster :

gdaldem slope input_dem output_slope_map”

[-p use percent slope (default=degrees)] [-s scale* (default=1)]

[-alg ZevenbergenThorne]

[-compute_edges] [-b Band (default=1)] [-of format] [-co”NAME=VALUE”]* [-q]

– To generate an aspect map from any GDAL-supported elevation raster

Outputs a 32-bit float raster with pixel values from 0-360indicating azimuth :

gdaldem aspect input_dem output_aspect_map”

[-trigonometric] [-zero_for_flat]

[-alg ZevenbergenThorne]

[-compute_edges] [-b Band (default=1)] [-of format] [-co”NAME=VALUE”]* [-q]

– To generate a color relief map from any GDAL-supported elevation raster

gdaldem color-relief input_dem color_text_fileoutput_color_relief_map

[-alpha] [-exact_color_entry | -nearest_color_entry]

[-b Band (default=1)] [-of format] [-co “NAME=VALUE”]* [-q]

where color_text_file contains lines of the format”elevation_value red green blue”

– To generate a Terrain Ruggedness Index (TRI) map from any GDAL-supportedelevation raster:

gdaldem TRI input_dem output_TRI_map

[-compute_edges] [-b Band (default=1)] [-of format] [-q]

– To generate a Topographic Position Index (TPI) map from any GDAL-supportedelevation raster:

gdaldem TPI input_dem output_TPI_map

[-compute_edges] [-b Band (default=1)] [-of format] [-q]

– To generate a roughness map from any GDAL-supported elevation raster:

gdaldem roughness input_dem output_roughness_map

[-compute_edges] [-b Band (default=1)] [-of format] [-q]

参数说明：

该工具提供了7个模块。分别如下：

hillshade（山地阴影）

用来生成山地阴影图像。

slope

用来生成坡度图像。

aspect

用来生成坡向图像。

color-relief

用来生成一个颜色渲染地图。

TRI

用来生成地表耐用指数（TRI）图像。

TPI

用来生成地形位置指数（TPI）图像。

roughness

用来生成粗糙度图像。

下面这些选项是通用选项：

input_dem:

输入DEM数据。

output_xxx_map:

输出处理的结果数据。

-of format:

输出文件格式，默认为GeoTiff。

-compute_edges:

(GDAL >= 1.8.0) 处理栅格边界数据和nodata附近的数据。

-alg ZevenbergenThorne:

(GDAL >= 1.8.0) 使用Zevenbergen & Thorne 规则代替Horn规则计算坡度坡向。研究表明Zevenbergen &Thorne用于平缓的地形，而Horn在处理比较粗糙的地形比较好。

-b band:

选择要处理的波段序号，默认为1.

-co “NAME=VALUE”:

创建输出文件选项，根据创建的文件格式不同而不同。

-q:

安静模式，不输出进度信息。

在所有的算法中，除了颜色渲染，其他的算法全部是基于一个3×3的窗口运算，如果在这个窗口中有一个值为nodata，那么结果值就是nodata，在GDAL1.8之后，如果指定了-compute_edges那么会在遇到nodata值时，会使用插值方式来确定nodata值，然后再进行处理。

算法模块

山地阴影

这个模块输出是一个8位栅格图，他在地形可视化中非常有用。您可以指定太阳高度角和太阳方位角海拔高度，垂直夸张因子和比例因子来改变结果。0作为输出的nodata值，具体的其他参数如下：

-z zFactor:

垂直夸张系数。

-s scale:

垂直系数的单位，如果水平单位是度，比如WGS84投影下的经纬度，可以将该参数设置为111120那么垂直单位是米，如果设置为370400，那么垂直单位为英尺。

-az azimuth:

光线方位角，度单位，0为正上方，90度为东边（右边）以此类推，默认值为315度。

-alt altitude:

光线的高度角，度单位，90度为DEM的正上方，0度为水平。

坡度

此模块将DEM栅格计算输出坡度值是32位浮数的栅格图像。你必须指定你想要的坡度值类型的选项：度或百分比坡度。同时也可以提供一个比例因子。-9999的值是用来作为输出NODATA值。以下具体的选项可供选择：

-p :

如果指定该选项，将使用坡度百分比计算，其他情况下使用度单位计算。

-s scale:

垂直系数的单位，如果水平单位是度，比如WGS84投影下的经纬度，可以将该参数设置为111120那么垂直单位是米，如果设置为370400，那么垂直单位为英尺。

坡向

该模块输出为一个32位浮点数的栅格数据，其值范围为0度到360度之间的坡向图， 0°表示坡向超正北方， 90°坡向朝向正东方，180°坡向朝向正南方，270°坡向朝向正西方。-9999表示无效值包括平地即坡度等于0的地方，详细的参数说明如下：

-trigonometric:

返回三角角度而不是方位，0°是指北，90°东，180°西，270°南。

-zero_for_flat:

返回0表示平地，而不是-9999。

如果使用这两个选项，使用gdaldem工具生成的坡度坡向结果同GRASS中的r.slope.aspect工具生成的结果一致。否则结果和Matthew Perry的aspect.cpp工具计算的结果一致。

color-relief

该功能是使用DEM数据和一个文本存储的颜色配置文件计算输出一个三波段（RGB）或者四波段（RGBA）的栅格图像。颜色配置文件中存储的是高程值和其对应的颜色值。默认情况下，给定的两个高程值之间的颜色使用该高程值对应的颜色插值来进行计算。参数-exact_color_entry或者参数-nearest_color_entry可以避免使用线性插值来确定从配置文件中没有找到对应的颜色值。

同时还有下面的选项可以使用：

color_text_file:

基于文本的颜色配置文件

-alpha :

在输出文件中添加Alpha通道

-exact_color_entry :

使用严格的匹配，在配置文件中搜索对应的高程值，如果没有找到匹配的颜色，使用“0,0,0,0”来进行填充。

-nearest_color_entry :

在颜色配置文件中使用与RGBA颜色相近的颜色进行填充。

color-relief是支持输出VRT格式的，在这种情形下，会将颜色配置文件转换为一个LUT（Look up table颜色查找表）。值得注意的是，使用百分比指定的高程值会被转换为绝对高程值。

基于文本格式的颜色配置文件一般每行都有四列，分别是：高程值，红，绿，蓝（0-255之间）。高程值可以是任何浮点数，或者使用nv来表示无效值。高程值统一可以使用百分比来表示，0%表示DEM中的最小值，100%表示DEM数据中的最大值。

此外还有一列是可选的就是Alpha部分，如果没有指定，默认值为255，表示不透明。同时每列之间的分隔符可以是：逗号（，），制表符（tab键），空格（），冒号（:），以上符合均为英文输入法下的符号。通常使用的颜色单词也可以用来替换RGB值，包括： white白色, black黑色, red红色, green绿色, blue蓝色, yellow黄色, magenta品红, cyan青色, aqua浅绿色, grey/gray灰色, orange橙色, brown棕色, purple/violet
紫色和indigo靛蓝色。

从GDAL1.8.0开始，GMT的.cpt调色板文件页可以支持（仅在COLOR_MODEL = RGB时可用）。

注意：颜色配置文件的句法同GRASS中的r.colors工具一致，ESRI的HDR颜色表文件（.clr）也是这种格式，Alpha部分的支持是GDAL的一个特殊扩展。

颜色配置文件示例：

3500 white

2500 235:220:175

50% 190 185 135

700 240 250 150

0 50 180 50

nv 0 0 0 0

TRI

该工具生成一个单波段的栅格图像，TRI是地形粗糙指数（Terrain Ruggedness Index）的缩写，代表的是中心点象元与周围8个象元点的高程差的绝对值的平均值。(参考Wilson et al 2007, Marine Geodesy 30:3-35)。

输出文件中的-9999代表无效值。该功能没有其他选项。

TPI

该工具生成一个单波段的栅格图像，TPI是地形方位指数（Topographic Position Index）的缩写，代表的是中心点象元高程与周围8个象元点的高程平均值的差值。(参考Wilson et al 2007, Marine Geodesy 30:3-35)。

输出文件中的-9999代表无效值。该功能没有其他选项。

roughness

该工具使用DEM数据计算一个单波段的栅格图像，粗糙度是指当前窗口中的最大值和最小值的差值。(参考Wilson et al 2007, Marine Geodesy 30:3-35)。

输出文件中的-9999代表无效值。该功能没有其他选项。

6. gdal_rasterize 矢量栅格化

用法：

gdal_rasterize [-b band]* [-i] [-at]

[-burn value]* | [-a attribute_name][-3d]

[-l layername]* [-where expression][-sql select_statement]

[-of format] [-a_srs srs_def] [-co”NAME=VALUE”]*

[-a_nodata value] [-init value]*

[-te xmin ymin xmax ymax] [-tr xresyres] [-tap] [-ts width height]

[-ot{Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}] [-q]

<src_datasource><dst_filename>

参数说明：

该程序将矢量几何图形（点，线，面）栅格化到栅格波段中。矢量必须是OGR库支持的矢量格式。注意的是矢量数据和栅格数据必须具有相同的坐标系统，不支持动态重投影功能。从GDAL1.8.0开始，支持创建结果图像，必须使用-tr或者-ts选项才行。

-b band:

要写入的栅格波段，可以使用-b参数来进行多选，默认是写入第一个波段中。

-i:

栅格反转，. Burn the fixed burn value, or the burn value associatedwith the first feature into all parts of the image not inside the provided apolygon.

-at:

Enables the ALL_TOUCHED rasterization option so that all pixels touched bylines or polygons will be updated not just those one the line render path, orwhose center point is within the polygon. Defaults to disabled for normalrendering rules.

-burn value:

栅格象元值，将所有的对象都用该值写入波段中，可以使用一个列表，来指定每个波段写入的象元值。

-a attribute_name:

指定属性字段中的字段值作为栅格值写入栅格文件中。

-3d:

将坐标点的Z值写入栅格中，假如参数-burn和-a attribute_name 被指定，那么该参数将会无效，目前仅仅支持3D的点和线数据。

-l layername:

指定矢量文件中的图层名作为输入要素，可以指定多次，但是至少有一个图层名或者-sql必须被指定。

-where expression:

SQL WHERE格式的查询语句，从指定的图层中选择特殊的要素。

-sql select_statement:

使用SQL语句来指定矢量数据中那些要素将要被栅格化。

-of format:

(GDAL >= 1.8.0)输出文件格式，默认为GeoTiff格式，使用短格式名称。

-a_nodata value:

(GDAL >= 1.8.0) 指定输出文件的无效值。

-init value:

(GDAL >= 1.8.0) 指定初始化栅格波段中的象元值，但是该值不能作为输出文件的无效值，如果只指定了一个值，那么其他的波段都会使用该值来填充。

-a_srs srs_def:

(GDAL >= 1.8.0) 覆盖输出文件的投影，如果没有指定，将使用输入的矢量文件中的投影作为输出投影，如果输入和输出投影不一致，将尝试重投影矢量。该参数的格式可以是GDAL/OGR支持的各种投影格式，如完整的WKT，PROJ.4，EPSG:n或者存储WKT的文本文件。

-co “NAME=VALUE”:

(GDAL >= 1.8.0) 创建图像选项，具体参考不同的栅格格式。

-te xmin ymin xmax ymax :

(GDAL >= 1.8.0) 设置地理范围，该值必须使用地理单位，如果没有指定，该范围同输入的矢量图层范围一致。

-tr xres yres :

(GDAL >= 1.8.0)设置输出分辨率，该值必须是有地理单位，两个必须都是正数。

-tap:

(GDAL >= 1.8.0) (target aligned pixels) align the coordinates of theextent of the output file to the values of the -tr, such that the aligned extentincludes the minimum extent.（不清楚该参数，同gdalwarp中的-tap）

-ts width height:

(GDAL >= 1.8.0) 设置输出文件大小，注意使用-ts是不能同时使用-tr。

-ot type:

(GDAL >= 1.8.0) 输出波段的数据类型，默认为64位浮点数。

-q:

(GDAL >= 1.8.0) 不输出进度信息和错误信息。

src_datasource:

任何OGR支持的矢量数据。

dst_filename:

GDAL支持的输出文件，必须支持更新模式访问，在GDAL1.8.0之前，该工具不能创建结果图像。

举例：

下面的例子，将mask.shp中的所有多边形栅格化为一个RGB的Tif图像work.tif，使用红色表示（RGB = 255,0,0）。

gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 0 -burn 0 -l mask mask.shpwork.tif

下面的例子将所有的class属性值为A的全部要素输出到一个高程文件中，高程文件的高程值为属性字段ROOF_H 的属性值。

gdal_rasterize -a ROOF_H -where ‘class=”A”‘ -l footprintsfootprints.shp city_dem.tif

7. gdaltransform 坐标系转换工具

用法：

gdaltransform [–help-general]

[-i] [-s_srs srs_def] [-t_srs srs_def] [-to”NAME=VALUE”]

[-order n] [-tps] [-rpc] [-geoloc]

[-gcp pixel line easting northing [elevation]]*

[srcfile [dstfile]]

参数说明：

gdaltransform工具用来转换坐标，从支持的投影，包括GCP点的变换。

-s_srs srs def:

原始空间参考，该值将被OGRSpatialReference.SetFromUserInput()函数调用，支持包括EPSG PCS 和GCSes (如EPSG:4296), PROJ.4声明或者一个后缀名为.prf的文件存储的WKT串。

-t_srs srs_def:

目标空间参考，该值将被OGRSpatialReference.SetFromUserInput()函数调用，支持包括EPSG PCS 和GCSes (如EPSG:4296), PROJ.4声明或者一个后缀名为.prf的文件存储的WKT串。

-to NAME=VALUE:

设置转换参数，将会使函数GDALCreateGenImgProjTransformer2()来调用。

-order n:

几何多项式变换次数（1到3），默认情况下会根据输入的GCP点的个数自动确定。

-tps:

强制将GCP点使用TPS转换方式转换。

-rpc:

强制使用RPC参数转换。

-geoloc:

强制使用Geolocation数组转换。

-i

逆转换，从目标到原始投影转换。

-gcppixel line easting northing [elevation]:

指定GCP点，通常至少需要三个。

srcfile:

原始投影下定义的GCP点文件，如果没有指定，需要从命令行-s_srs或者-gcp参数输入。

dstfile:

目标投影定义文件。

读取的坐标必须是成对出现（或者三个一组），输出也以同样的方式输出，所有的变换共gdalwarp相同，包括使用GCP的变换。注意输入输出必须使用十进制小数，当前不支持DMS（度分秒）的输入输出。如果指定了输入图像，那么输入的坐标是基于图像的行列号，如果输出图像指定，输出的坐标也是图像的行列号。

举例：

重投影示例

简单的投影转换，从一个投影转换到另一个投影：

gdaltransform -s_srs EPSG:28992 -t_srs EPSG:31370

177502 311865

Produces the following output in meters in the “Belge 1972 / BelgianLambert 72” projection:

244510.77404604 166154.532871342 -1046.79270555763

图像RPC变换示例：

下面的命令行使用基于RPC变换，并且使用-i选项来将经纬度坐标转为图像的行列号：

gdaltransform -i -rpc 06OCT20025052-P2AS-005553965230_01_P001.TIF

125.67206 39.8530750

Produces this output measured in pixels and lines on the image:

3499.49282422381 2910.83892848414 50

8. nearblack

用法：

gdalinfo [–help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]

[-noct] [-nofl][-checksum] [-proj4] [-mdd domain]*

[-sd subdataset] datasetname

参数说明：

该工具没有进行说明。

9. gdal_grid 插值工具

用法：

gdal_grid [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}]

[-of format] [-co “NAME=VALUE”]

[-zfield field_name]

[-a_srs srs_def] [-spat xmin ymin xmax ymax]

[-clipsrc <xminymin xmax ymax>|WKT|datasource|spat_extent]

[-clipsrcsqlsql_statement] [-clipsrclayer layer]

[-clipsrcwhereexpression]

[-l layername]* [-where expression] [-sql select_statement]

[-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]

[-a algorithm[:parameter1=value1]*] [-q]

<src_datasource> <dst_filename>

参数说明：

该程序从OGR支持的矢量数据中创建一个规则的栅格数据。输入数据必须有用来插值的数据，可以选择下面的插值方法。

-ot type:

指定输出文件个数据类型。

-of format:

输出文件格式，默认为GeoTiff格式。

-txe xmin xmax:

设置输出文件的地理范围X。

-tye ymin ymax:

设置输出文件的地理范围Y。

-outsize xsize ysize:

设置输出文件的大小。

-a_srs srs_def:

设置输出文件的投影，可以是WKT串，PROJ.4，EPSG:n或者含有wkt的文件。

-zfield field_name:

指定矢量要素中获取Z值的字段。

-a [algorithm[:parameter1=value1][:parameter2=value2]…]:

设置插值算法或者数据单位参数，具体参考 INTERPOLATION ALGORITHMS 和DATA METRICS 。

-spat xmin ymin xmax ymax:

指定空间范围进行过滤(xmin, ymin) – (xmax, ymax)。

-clipsrc [xmin ymin xmax ymax]|WKT|datasource|spat_extent:

添加一个空间过滤来选择处理的要素，包括空间范围（SRS表示），WKT几何形状（POLYGON 或MULTIPOLYGON），也可以使用下面的参数来进行过滤：-clipsrclayer，-clipsrcwhere或-clipsrcsql。

-clipsrcsql sql_statement:

使用SQL语句对原始数据进行查询过滤。

-clipsrclayer layername:

指定源数据中图层名称。

-clipsrcwhere expression:

使用属性表来限制几何体。

-l layername:

指定图层名称。

-where expression:

使用SQL WHERE方式来过滤图层中的数据。

-sql select_statement:

使用SQL语句来过滤图层中的数据。

-co “NAME=VALUE”:

创建文件参数选项。

-q:

不输出进度信息和错误信息。

src_datasource:

输入矢量，任何OGR库支持的矢量数据。

dst_filename:

输出栅格，GDAL支持的栅格数据。

INTERPOLATION ALGORITHMS（插值算法）

插值算法可以选择下面几种：

invdist

距离次方反比法，这个是默认算法，该算法还有下列参数：

power:

权重指数，默认为2.0。

smoothing:

平滑参数，默认为0.0。

radius1:

搜索椭圆第一半径（X轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

radius2:

搜索椭圆第二半径（Y轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

angle:

搜索椭圆的旋转角度，使用度单位。默认为0.0，逆时针方向。

max_points:

使用的数据点的最大个数，不会搜索多余这个数字的点，如果搜索椭圆设置将会使用该值（两个半径均不为0时）。默认为0，为0表示使用所有的点。

min_points:

使用的最小数据点个数，如果数据点小于该值，并且找到的网络节点值为空，将使用Nodata值来标记。该值仅仅当搜索椭圆设置（即两个半径均不为0），默认该值为0。

nodata:

设置NODATA值，默认为0.0。

average

移动均值算法，该算法参数如下：

radius1:

搜索椭圆第一半径（X轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

radius2:

搜索椭圆第二半径（Y轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

angle:

搜索椭圆的旋转角度，使用度单位。默认为0.0，逆时针方向。

min_points:

使用数据点的最小个数，如果数据点小于该值，并且找到的网络节点值为空，将使用Nodata值来标记，默认该值为0。

nodata:

设置NODATA值。默认值为0.0。

注意：使用移动均值算法必须设置搜索椭圆，当计算网格节点值时使用窗口的均值。

nearest

最邻近算法。该算法参数如下：

radius1:

搜索椭圆第一半径（X轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

radius2:

搜索椭圆第二半径（Y轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

angle:

搜索椭圆的旋转角度，使用度单位。默认为0.0，逆时针方向。

nodata:

设置NODATA值。默认值为0.0。

DATA METRICS（数据单位）

Besides the interpolation functionality gdal_grid 除了插值功能外，还有使用指定窗口和输出网格来计算一些数据指标，这些指标有：

minimum:

搜索椭圆中网格节点的最小值。

maximum:

搜索椭圆中网格节点的最大值。

range:

在搜索椭圆中查找到的最大值和最小值的差值。

count:

在搜索椭圆中找到的点个数。

average_distance:

网格节点（搜索椭圆中心）和所有在网格节点搜索椭圆的数据点之间的平均距离。

average_distance_pts:

在网格节点搜索椭圆中发现的数据点之间的平均距离。椭圆内每对点之间距离的计算和所有的距离平均是作为一个网格节点的值设置。

所有的指标都具有下面相同的选项：

radius1:

搜索椭圆第一半径（X轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

radius2:

搜索椭圆第二半径（Y轴的旋转角度为0），设置该参数为0时将使用所有的点数组，默认为0.0。

angle:

搜索椭圆的旋转角度，使用度单位。默认为0.0，逆时针方向。

min_points:

使用的数据点的最小个数。如果找到的数据点小于将设置为NODATA值，该值仅仅用于设置搜索椭圆时可用（两个半径均不为0时）。默认值为0。

nodata:

设置NODATA值。默认值为0.0。

READING COMMA SEPARATED VALUES（读取逗号分割符文件CSV）

通常使用一个文本文件，使用逗号分隔XYZ值，该文件通常叫做CSV文件。可以在gdal_grid中使用该文件。你需要做的就是要对你的CSV文件创建一个VRT文件头，并将该文件作为gdal_grid的输入数据。关于VRT格式详情参考 Virtual Format 页面。

下面有个小例子。假如我们有个CSV文件叫dem.csv ，里面内容为：

Easting,Northing,Elevation

86943.4,891957,139.13

87124.3,892075,135.01

86962.4,892321,182.04

87077.6,891995,135.01

…

对于上面的数据，创建一个dem.vrt 的头文件，内容如下：

<GeometryType>wkbPoint</GeometryType>

</OGRVRTLayer>

</OGRVRTDataSource>

该说明是一个叫做2.5D的几何数据，含有三个坐标XYZ，Z值将被用来插值。现在可以使用dem.vrt来进行测试，可以先使用ogrinfo来测试该文件是否可用。数据中包括一个叫dem的点图层，数据包括在CSV文件中。使用这种方式，可以处理超过三列，交换列的CSV文件。

如果你的CSV文件不包括列头，那么可以在VRT头文件中加入下面的内容来进行处理：<GeometryFieldencoding=”PointFromColumns” x=”field_1″y=”field_2″ z=”field_3″/>

Comma Separated Value 页面描述了OGR库对CSV文件的支持描述信息。

举例：

下面将创建一个TIFF的图像数据，使用VRT数据源，从一个CSV文件中，使用反距离次方算法来进行插值。

gdal_grid -a invdist:power=2.0:smoothing=1.0 -txe 85000 89000 -tye 894000890000 -outsize 400 400 -of GTiff -ot Float64 -l dem dem.vrt dem.tiff

下面的命令行同上面的类似，但是读取的数据是使用-zfield 选项来指定的。

gdal_grid -zfield “Elevation” -a invdist:power=2.0:smoothing=1.0-txe 85000 89000 -tye 894000 890000 -outsize 400 400 -of GTiff -ot Float64 -ldem dem.vrt dem.tiff

10. gdallocationinfo 栅格查询工具

用法：

gdallocationinfo [–help-general] [-xml] [-lifonly] [-valonly]

[-b band]* [-l_srs srs_def] [-geoloc] [-wgs84]

srcfile [x y]

参数说明：

Gdallocationinfo程序通过指定一种坐标系中的坐标来进行对其位置的像元值进行输出显示。

-xml:

输出XML格式，方便后期处理。

-lifonly:

输出LocationInfo 中的信息。

-valonly:

仅仅输出指定位置的每个波段的像元值。

-b band:

选择一个波段进行处理，多波段文件将被列出显示。默认是查询所有的波段。

-l_srs srs def:

指定输入的x,y坐标的坐标系。

-geoloc:

表示输入的x,y坐标是地理参考坐标系。

-wgs84:

表示输入的x,y是WGS84坐标系下的经纬度坐标。

srcfile:

输入栅格图像的名称。

查询的X坐标。默认为图像行列号，如果使用-l_srs,-wgs84或者-geoloc时按照指定的坐标来处理。

查询的Y坐标。默认为图像行列号，如果使用-l_srs,-wgs84或者-geoloc时按照指定的坐标来处理

该工具的目的是输出一个像素的各种信息。目前支持下面四项：

· 像素的行列号。

· 输出元数据中的LocationInfo 信息，目前只支持VRT文件。

· 全部波段或子文件中的所有波段中的像元值。

· 未缩放的像元值，如果对波段进行了缩放和偏移操作。

输入x和y坐标是，是通过命令行来进行输入的，一般是先x后y，即先列号后行号；如果使用-geoloc, -wgs84, 或-l_srs 选项后，会自动根据输入进行交换。默认的输出信息是纯文本文件，也可以使用-xml选项将其使用xml格式进行输出。将来会添加其他的信息。

11. gdalsrsinfo 将给定的CRS按照不同的格式显示

用法：

gdalsrsinfo [options] srs_def

列出给定的CRS的各种表示方法。如ESRI的wkt格式，OGC的wkt格式和Proj4的格式等。

srs_def 可以是一个GDAL或者OGR支持的文件，或者是任何GDAL和OGR支持的CRS格式（包括WKT,PROJ4,EPSG:n等）。

参数说明：

Options:

[–help-general] [-h] 显示帮助信息并退出

[-p] 输出信息（如WKT字串格式）

[-V] 验证SRS

[-o out_type] 输出类型，默认为all，所有类型见下：

{all,proj4,wkt,wkt_simple,wkt_old,wkt_esri,mapinfo,xml}

举例：

gdalsrsinfo -o all “EPSG:4326”

PROJ.4 : ‘+proj=longlat +datum=WGS84 +no_defs ‘

OGC WKT :

GEOGCS[“WGS 84”,

DATUM[“WGS_1984”,

SPHEROID[“WGS84”,6378137,298.257223563,

AUTHORITY[“EPSG”,”7030″]],

AUTHORITY[“EPSG”,”6326″]],

PRIMEM[“Greenwich”,0,

AUTHORITY[“EPSG”,”8901″]],

UNIT[“degree”,0.0174532925199433,

AUTHORITY[“EPSG”,”9122″]],

AUTHORITY[“EPSG”,”4326″]]

OGC WKT (simple) :

GEOGCS[“WGS 84”,

DATUM[“WGS_1984”,

SPHEROID[“WGS84”,6378137,298.257223563]],

PRIMEM[“Greenwich”,0],

UNIT[“degree”,0.0174532925199433]]

ESRI WKT :

GEOGCS[“GCS_WGS_1984”,

DATUM[“D_WGS_1984”,

SPHEROID[“WGS_1984”,6378137,298.257223563]],

PRIMEM[“Greenwich”,0],

UNIT[“Degree”,0.017453292519943295]]

12. gdal-config GDAL配置参数

用法：

gdal-config [OPTIONS]

Options:

[–prefix[=DIR]]

[–libs]

[–cflags]

[–version]

[–ogr-enabled]

[–formats]

参数说明：

该脚本工具可以用来查看关于GDAL的一些安装信息等内容。

–prefix:

显示GDAL的安装目录。

–libs:

显示GDAL所依赖的库。

–cflags:

编译GDAL时包括的宏定义等信息。

–version:

显示GDAL版本信息。

–ogr-enabled:

显示是否包含OGR库，如果输出”yes”说明OGR可以使用，如果是”no”则不能使用OGR库中的东西。

–formats:

输出GDAL支持的格式。

一、 GDAL Python工具

本文主要介绍的是GDAL工具集中的Python脚本命令，需要的环境必须是有Python环境和GDAL的Python版本。这是必须的，否则这些工具都不能用。对于已经安装ArcGIS的同学来说，Python都已经安装好了，可以直接下载GDAL的Python版本，然后就可以使用下面这些工具。

1. rgb2pct.py 转换24位RGB图为8位图

用法：

rgb2pct.py [-n colors | -pct palette_file] [-of format] source_file dest_file

参数说明：

该工具会自动根据指定的RGB图像计算最合适的假彩色颜色表。然后对结果影像使用该颜色表。简而言之，就是将RGB彩色图像转换为一个单波段的图像，使用颜色表来表示颜色。

-n colors:

指定生成颜色表的颜色数目，默认是256，其值必须是2到256之间的整数值。

-pct palette_file:

从调色板文件中提取颜色表而不从图像中计算。调色板文件必须是GDAL支持的调色板格式。

-of format:

输出文件格式，默认为GeoTiff格式，而且输出格式必须支持颜色表。

source_file:

输入的RGB图像。

dest_file:

输出的图像路径。如果图像不存在会创建一个。

举例：

如果希望指定调色板信息，比如简单的文本格式，如GDAL VRT格式，在下面的例子中将指定一个使用文本编辑器创建的VRT格式的调色板文件，一共有四个颜色，RGBA值分别是：238/238/238/255，237/237/237/255，36/236/236/255和229/229/229/255。

% rgb2pct.py -pct palette.vrt rgb.tif pseudo-colored.tif

% more < palette.vrt

<ColorInterp>Palette</ColorInterp>

</ColorTable>

</VRTRasterBand>

</VRTDataset>

2. pct2rgb.py 转换8位图为24位RGB图

用法：

pct2rgb.py [-of format] [-b band] [-rgba] source_file dest_file

参数说明：

该工具是将一个带有颜色表的图像转为RGB图像。

-of format:

输出文件格式，默认为GeoTiff格式。

-b band:

指定要转换的波段序号，默认是第一个波段。

-rgba:

生成RGBA文件（默认是生成RGB文件）。

source_file:

输入文件。

dest_file:

输出文件。

3. gdal_merge.py 镶嵌图像

用法：

gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*

[-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-v] [-pct]

[-ul_lr ulx uly lrx lry] [-n nodata_value] [-init “value [value…]”]

[-ot datatype] [-createonly] input_files

参数说明：

该工具会自动镶嵌指定的图像。所有的图像必须具有相同的坐标系统和相同的波段数目；输入图像可能是有重叠区或者不同的分辨率。在重叠区部分最后的图像将会覆盖之前的图像值。

-o out_filename:

输出文件，如果不存在将会创建一个新图像。如果该值不指定，将会创建一个叫out.tif的图像。

-of format:

输出文件格式，默认为GeoTiff格式。

-co NAME=VALUE:

创建图像选项，具体参考具体图像格式说明。

-ot datatype:

指定输出数据类型，使用类型名称，如Byte，Int16等。

-ps pixelsize_x pixelsize_y:

输出图像的象元大小，如果不指定，将以第一个图像的分辨率为基准。

-tap:

(GDAL >= 1.8.0) (target aligned pixels) align the coordinates of theextent of the output file to the values of the -tr, such that the alignedextent includes the minimum extent.

-ul_lr ulx uly lrx lry:

输出文件范围，如果不指定，将会以所有的输入图像的范围并集为输出文件范围。

-v:

生成详细的镶嵌操作信息，当执行结束的时候。

-separate:

Place each input file into a separate stacked band.

-pct:

提取颜色表从第一个输入图像，并在结果图像中使用，按照这种方式镶嵌图像的话，后面的图像的颜色表都会被忽略，以第一个为准。

-n nodata_value:

指定nodata值，在镶嵌操作的时候会忽略该值。

-init “value(s)”:

指定输出文件的初始值，不会在输出文件中指定nodata值，如果只指定了一个值，那么所有的波段都会使用该值来作为初始值。

-createonly:

输出文件已经创建好了。但是输入文件数据没有写入。

举例：

创建一个图像，将所有的波段都指定为255。

% gdal_merge.py -init 255 -o out.tif in1.tif in2.tif

创建一个RGB图像，并将前两个波段初始化为0，第三个波段为255。

% gdal_merge.py -init “0 0 255” -o out.tif in1.tif in2.tif

4. gdal2tiles.py 生成TMS切片

用法：

gdal2tiles.py [-title “Title”] [-publishurlhttp://yourserver/dir/]

[-nogooglemaps] [-noopenlayers] [-nokml]

[-googlemapskey KEY] [-forcekml] [-v]

input_file [output_dir]

参数说明：

该工具按照OSGeo的切片地图服务说明书来生成小的切片数据和元数据。基于Google地图和OpenLayers简单的网页地图浏览。

GDAL2Tiles 工具创建Google Earth (KMLSuperOverlay)必需的元数据文件，按照EPSG:4326投影。同时会创建世界文件（World file），但是页可以创建其他的投影的文件。

-title “Title”:

生成元数据，网页视图和KML文件。

-publishurl http://yourserver/dir/:

添加发布的网页地址，会将结果上传到该地址。

-nogooglemaps:

不要生成基于Google地图的html页面。

-noopenlayers:

不要生成基于OpenLayers的html页面。

-nokml:

不要生成Google Earth的kml文件。

Do not generate KML files for Google Earth.

-googlemapskey KEY:

指定生成区域的范围，使用Google Maps的API。参考(http://www.google.com/apis/maps/signup.html).

-forcekml

生成kml文件，输入图像必需是EPSG:4326的投影！

-v

处理结束后输出详细信息。

5. gdal_retile.py

用法：

gdal_retile.py [-v] [-co NAME=VALUE]* [-of out_format] [-ps pixelWidthpixelHeight]

[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}]’

[ -tileIndex tileIndexName [-tileIndexField tileIndexFieldName]]

[ -csv fileName [-csvDelim delimiter]]

[-s_srs srs_def] [-pyramidOnly]

[-r {near/bilinear/cubic/cubicspline/lanczos}]

-levels numberoflevels

[-useDirForEachRow]

-targetDir TileDirectory input_files

参数说明：

该工具对输入图像创建分块。所有的输入图像必需具有相同的投影和相同的波段数。可以选择生成金字塔。会生成输出瓦片的矢量文件。

-targetDir directory:

输出瓦片存放的目录。金字塔文件会存放在子文件夹中，子文件夹的命名从1开始，数字代表金字塔级别。创建的输出名称由源文件名称和一个序号组成。

-of format:

输出文件格式，默认为GeoTiff格式。

-co NAME=VALUE:

创建选项，具体参考具体的格式说明。

-ot datatype:

指定输出图像的数据类型，使用类型名称，如Byte，Int16，UInt16等。

-ps pixelsize_x pixelsize_y:

输出文件大小，如果不指定，默认是256*256。

-levels numberOfLevels:

创建的金字塔级数。

-v:

输出相信的信息。

-pyramidOnly:

不建立瓦片文件，只建立金字塔。

-r algorithm:

重采样方式，默认为最邻近采样。

-s_srs srs_def:

源文件的空间参考。坐标系统可以是任何OGRSpatialReference.SetFro函数支持的坐标系统，包括EPSG PCS和GCSes，PROJ.4声明或者其他的包含WKT字符串的prf文件。如果不知道，将会从指定输入图像中读取。同时该空间参考用来创建生成的瓦片的矢量文件。

-tileIndex tileIndexName:

包含结果瓦片索引的shp名称。

-tileIndexField tileIndexFieldName:

包含瓦片名称的属性表字段名称。

-csv csvFileName:

指定包含瓦片地理信息的csv文件名称。该文件包含五列，分别是：（瓦片名称，最小x，最大x，最小y，最大y）tilename,minx,maxx,miny,maxy。

-csvDelim column delimiter:

CSV文件的列之间分隔符，默认使用分号风格。

-useDirForEachRow:

通常情况下，输出的瓦片文件存放在-targetDir指定的文件夹中。对于大图像，在有些文件系统中对于文件夹中的文件数过多可能会有一些问题，从而导致gdal_retile程序不能正常执行完毕。使用该参数可以创建不同的目录结构。原始的瓦片文件存放在名字为0的子文件夹中，金字塔依次存放在子文件1，2，3…中，数字代表金字塔的级别。

6. gdal_proximity.py

用法：

gdal_proximity.py srcfile dstfile [-srcband n] [-dstband n]

[-of format] [-co name=value]*

[-ot Byte/Int16/Int32/Float32/etc]

[-values n,n,n] [-distunits PIXEL/GEO]

[-maxdist n] [-nodata n] [-fixed-buf-val n]

参数说明：

gdal_proximity.py脚本用来生成栅格图像的距离图，即计算中心点的象元与临近的象元之间的距离作为结果象元值。

srcfile

输入的栅格图像。=

dstfile

输出影像，可以是一个已经存在的文件，大小和原始图像一致，如果不存在，将会创建一个新图像。

-srcband n

指定用来计算的波段序号，默认为第一个波段。

-dstband n

指定输出的波段序号，默认为第一个波段。

-of format:

输出图像格式，默认为GeoTiff格式，使用短格式名称。

-co “NAME=VALUE”:

创建图像选项，具体参考对应格式的说明。

-ot datatype:

输出文件的数据类型，使用短名称，如Byte，Int16。

-values n,n,n:

A list of target pixel values in the source image to be considered targetpixels. If not specified, all non-zero pixels will be considered target pixels.指定结果象元值在原始影像如果没有指定，所有的非0象元将被用来计算结果图。

-distunits PIXEL/GEO:

指定距离的单位，默认为像素单位，可以为像素单位或者地理坐标单位。

-maxdist n:

最大距离值，所有超过该值的象元具体将被设置为nodata值或者65535。距离使用像素单位时，如果使用GEO参数的话，不会执行该操作。

-nodata n:

指定nodata值。

-fixed-buf-val n:

Specify a value to be applied to all pixels that are within the -maxdistof target pixels (including the target pixels) instead of a distance value.

7. gdal_polygonize.py 栅格矢量化工具

用法：

gdalinfo [–help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]

[-noct] [-nofl][-checksum] [-proj4] [-mdd domain]*

[-sd subdataset] datasetname

参数说明：

gdalinfo程序输出gdal支持的栅格格式的一系列信息。

-mm

强制计算栅格每个波段的最大最小值。

-stats

读取和现实图像统计信息，如果指定该参数，将强制计算图像的统计信息，如各个波段的最大值、最小值、均值、标准差等。

举例：

Generate polygons from raster.

8. gdal_sieve.py 小斑去除滤波工具

用法：

gdalinfo [–help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]

[-noct] [-nofl][-checksum] [-proj4] [-mdd domain]*

[-sd subdataset] datasetname

参数说明：

gdalinfo程序输出gdal支持的栅格格式的一系列信息。

-mm

强制计算栅格每个波段的最大最小值。

-stats

读取和现实图像统计信息，如果指定该参数，将强制计算图像的统计信息，如各个波段的最大值、最小值、均值、标准差等。

举例：

Raster Sieve filter.

9. gdal_fillnodata.py 填充NoData区域

用法：

gdalinfo [–help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]

[-noct] [-nofl][-checksum] [-proj4] [-mdd domain]*

[-sd subdataset] datasetname

参数说明：

gdalinfo程序输出gdal支持的栅格格式的一系列信息。

-mm

强制计算栅格每个波段的最大最小值。

-stats

读取和现实图像统计信息，如果指定该参数，将强制计算图像的统计信息，如各个波段的最大值、最小值、均值、标准差等。

二、其他说明

关于GDAL的Python编译、安装和使用方法，网上有很多的博文。由于Python是不用编译的脚本语言，所以不用编译就可以运行，语法简单，受到很多人的喜爱，我也是Python的爱好者。尤其是对于一些简单的功能，完全没有必要写一个C或者C++以及其他的语言程序，还需要编译后才能使用，总觉得不是那么的方便。如果用Python，写完后直接可以运行。

关于Python版本的编译和安装，参考此系列的博文三：http://blog.csdn.net/liminlu0314/article/details/6945452；此外，对于Python版本GDAL的使用，主要参考李林（lilin）写的《GDAL库学习笔记》系列文章（想当年我就是看这七篇文章来学习GDAL的）。链接地址为：http://wiki.woodpecker.org.cn/moin/lilin。

也可以直接下载编译好的Python版本GDAL，具体地址为：http://pypi.python.org/pypi/GDAL/

一、简单的调用

关于GDAL的使用，网上的资料都很多，主要还是要熟悉GDAL的组织结构，类以及类的函数等，熟悉了，使用GDAL就不在话下了。最常用的就是动态库的GDAL，当然你也可以使用静态库，这里只是简单的介绍使用动态GDAL库来做开发。

首先打开VS，新建一个工程，控制台的就成。然后在工程的属性对话框中，找到【配置属性】-【C/C++】-【常规】，右侧的【附加包含目录】中，将GDAL的include文件夹路径填写到这里，如下图：

第二、继续在属性对话框中，找到【配置属性】-【链接器】-【常规】，右侧的【附加库目录】中，将GDAL的lib文件夹路径填写到这里，如下图：

第三、在【配置属性】-【链接器】-【输入】，右侧的【附加依赖项】中，将gdal_i.lib填写到此处。然后点击确定即可。至此，使用GDAL的环境全部搭建完成，剩下的就是在您的代码中使用GDAL了。

将下面的代码（代码摘自GDAL官方指南：http://gdal.org/gdal_tutorial.html）贴到刚才新建的工程中的cpp文件中，保存后编译，正常情况下会提示生成成功，然后运行，会在控制台上将图像的信息输出。

[cpp] view plaincopyprint?

#include “gdal_priv.h”
#include “cpl_conv.h” //for CPLMalloc()
int main()
{
//注册文件格式
GDALAllRegister();
const char* pszFile = “C:\\Test.img”;
GDALDataset *poDataset;
//使用只读方式打开图像
poDataset = (GDALDataset*) GDALOpen( pszFile,GA_ReadOnly );
if( poDataset == NULL )
{
printf( “File: %s不能打开！\n”,pszFile);
return 0;
}
//输出图像的格式信息
printf( “Driver:%s/%s\n”,
poDataset->GetDriver()->GetDescription(),
poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME) );
//输出图像的大小和波段个数
printf( “Size is%dx%dx%d\n”,
poDataset->GetRasterXSize(),poDataset->GetRasterYSize(),
poDataset->GetRasterCount());
//输出图像的投影信息
if( poDataset->GetProjectionRef() != NULL )
printf( “Projectionis `%s’\n”, poDataset->GetProjectionRef() );
//输出图像的坐标和分辨率信息
double adfGeoTransform[6];
if( poDataset->GetGeoTransform( adfGeoTransform) == CE_None )
{
printf( “Origin =(%.6f,%.6f)\n”,
adfGeoTransform[0], adfGeoTransform[3]);
printf( “PixelSize = (%.6f,%.6f)\n”,
adfGeoTransform[1], adfGeoTransform[5]);
}
GDALRasterBand *poBand;
int nBlockXSize, nBlockYSize;
int bGotMin, bGotMax;
double adfMinMax[2];
//读取第一个波段
poBand = poDataset->GetRasterBand( 1 );
//获取图像的块大小并输出
poBand->GetBlockSize(&nBlockXSize, &nBlockYSize );
printf( “Block=%dx%dType=%s, ColorInterp=%s\n”,
nBlockXSize, nBlockYSize,
GDALGetDataTypeName(poBand->GetRasterDataType()),
GDALGetColorInterpretationName(
poBand->GetColorInterpretation()));
//获取该波段的最大值最小值，如果获取失败，则进行统计
adfMinMax[0] = poBand->GetMinimum( &bGotMin);
adfMinMax[1] = poBand->GetMaximum( &bGotMax);
if( ! (bGotMin&& bGotMax) )
GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
printf( “Min=%.3fd,Max=%.3f\n”, adfMinMax[0], adfMinMax[1] );
//输出图像的金字塔信息
if( poBand->GetOverviewCount() > 0 )
printf( “Band has%d overviews.\n”, poBand->GetOverviewCount() );
//输出图像的颜色表信息
if( poBand->GetColorTable() != NULL)
printf( “Band hasa color table with %d entries.\n”,
poBand->GetColorTable()->GetColorEntryCount() );
float *pafScanline;
int nXSize = poBand->GetXSize();
//读取图像的第一行数据
pafScanline = (float*) CPLMalloc(sizeof(float)*nXSize);
poBand->RasterIO(GF_Read, 0, 0, nXSize,1,
pafScanline, nXSize,1, GDT_Float32, 0, 0 );
CPLFree(pafScanline);
//关闭文件
GDALClose((GDALDatasetH)poDataset);
}

二、GDAL源代码调试

在很多时候我们需要看看GDAL的内部实现，当然可以直接查看GDAL的源代码，但是直接看源代码，不能很好的理解，这时候就需要调试查看源代码中变量的内容。调试GDAL的源代码，需要GDAL的debug版本，以及编译GDAL的时候的pdb等调试文件。当然也可以把GDAL的源代码加入到你的工程中，但是这样太费时费力。

下面就有一个很简单的方法，可以直接调试进GDAL的源代码中，首先编译一下GDAL的debug版本，将编译生成的文件，主要有gdal18.dll，gdal_i.exp，gdal_i.lib，gdal.lib，gdal18.pdb，gdal18.ilk，gdal18.exp等文件，将gdal18开头的文件拷贝到自己工程的生成目录中，然后调试自己的程序，在执行到GDALOpen函数（或者其他GDAL的函数）时按F11键，就会进入到GDAL的源代码中进行调试GDAL代码。

三、GDAL使用示例

1、使用GDAL进行图像裁切，参考http://blog.csdn.net/liminlu0314/article/details/6136512

2、使用GDAL进行图像重采样，参考http://blog.csdn.net/liminlu0314/article/details/6130064

3、使用GDAL创建金字塔，参考http://blog.csdn.net/liminlu0314/article/details/6127755

一、关于RasterIO

在GDAL中读写图像是最基本的操作，那么RasterIO也就是最基本的函数了，关于RasterIO有很多方式，这个函数的功能相当强大，下面慢慢说明。RasterIO一共有两个，一个是GDALRasterBand::RasterIO，另一个是GDALDataset::RasterIO，这两个RasterIO都可以对图像数据来进行读写，大多数情况下是一样的，但是还是有一些区别的。

二、RasterIO参数说明

下面对两个RasterIO的参数进行一个简单的说明：首先是GDALRasterBand::RasterIO ，该函数的声明如下，具体可以参考下面网址：http://gdal.org/classGDALRasterBand.html#5497e8d29e743ee9177202cb3f61c3c7

[cpp] view plaincopyprint?

CPLErr GDALRasterBand::RasterIO ( GDALRWFlag eRWFlag,
int nXOff,
int nYOff,
int nXSize,
int nYSize,
void * pData,
int nBufXSize,
int nBufYSize,
GDALDataType eBufType,
int nPixelSpace,
int nLineSpace
)

接下来是GDALDataset::RasterIO ，该函数的什么如下，具体形式可以参考下面的网址：http://gdal.org/classGDALDataset.html#e077c53268d2272eebed10b891a05743

[cpp] view plaincopyprint?

CPLErr GDALDataset::RasterIO ( GDALRWFlag eRWFlag,
int nXOff,
int nYOff,
int nXSize,
int nYSize,
void * pData,
int nBufXSize,
int nBufYSize,
GDALDataType eBufType,
int nBandCount,
int * panBandMap,
int nPixelSpace,
int nLineSpace,
int nBandSpace
)

由于这两个函数的参数基本一致，大多数参数都是一样的，下面就一起进行说明：

第一个参数RWFlag来指定是读数据还是写入数据，其值只能有两个即：GF_Read 和GF_Write，分别表示读取数据和写入数据。

第二个和第三个参数nXOff, nYOff表示读取或者写入图像数据的起始坐标图像的左上角坐标为（0，0）。

第四个和第五个参数 nXSize, nYSize表示读取或者写入图像数据的窗口大小，nXSize表示宽度，nYSize表示高度，均使用像素为单位，该宽度和高度是从第二个和第三个参数处开始计算。这两个参数和第二第三个参数一起表示就是，读取和写入图像的窗口位置和大小。

第六个参数pData是指向存储数据的一个指针。如果是写入数据，那么会将pData中的数据写入到栅格图像中去；如果是读取数据，那么会将栅格数据中的数据读入到pData中。pData的真实数据类型是通过后面的eBufType参数来指定的。如GDT_Byte就是代表的是一个8U的数据类型，如果是GDT_Float32就表示的是一个32F（float）的数据类型。RasterIO会自动将读入的数据按照参数eBufType指定的数据类型进行转换，需要注意的是，将浮点数转换为整数时，将对数据进行四舍五入处理；而且从一个大的存储单位转换到一个较小的存储单位是所进行的操作是截断操作而不是按照比例缩小操作，比如原来的实际数据中float的存取数据范围超过了缓冲数据指定的数据byte类型的最大存储范围，操作将把超过byte存储的范围外的数据进行丢弃处理，而不是将float缩小到0～255。

第七个和第八个参数nBufXSize和nBufYSize参数指定缓冲区的大小。注意pData的大小应当是nBufXSize×nBufYSize。当读取的数据是完整分辨率的数据（原始数据，没有进行缩放操作），他们应该设置和取值窗口的大小相同，也就是与第四个和第五个参数相等，但是在读取时使用了缩小或者放大系数，那么他们需要根据这个缩放系数进行调整。在这种情况下，RasterIO将会使用缩略图组overviews（金字塔）中某个合适的缩略图来进行读取数据。

第九个和第十个参数（对于GDALDataset来说还有第十一个参数nBandSpace）nPixelSpace和nLineSpace（以及nBandSpace）参数一般情况下是将0作为缺省值。但是，他们可以用于控制存取的内存数据的排列顺序，可以使用这两个参数将图像数据按照另一种组织形式读取内存缓冲区中。也就是说这两个（三个）参数可以读取或者写入非常规组织的缓冲数据。这个首先可以用于在一个缓冲区中包含多个波段数据，并且各个数据之间是交叉排列的，比如一个图像中的数据组织是RGBRGBRGB…，而普通的数据可能是RRR…GGG…BBB…，我们一般读取到的数据就是RGBRGBRGB…这种排列，现在需要使用RRR…GGG…BBB…这样的排列，一般想法就是自己写个for循环之类的，重新组织一次，其实完全没有必要，只要设置这两个（三个）参数就可以达到这个目的。

nPixelSpace表示的是在一个扫描行中一个像元的字节偏移起始点到下一个像元字节偏移起始点之间的字节间隔，如果默认使用0，那么将使用eBufType作为实际的两个像元之间的字节间隔。

nLineSpace表示在一行数据和下一行数据之间的起始字节见的间隔，如果使用0，那么将会使用eBufType*nBufXSize来表示实际间隔。

此外如果使用的GDALDataset::RasterIO函数,最后还有一个参数叫nBandSpace,同上,这个参数的意思就是一个波段与下一个波段之间的起始字节间的间隔,如果使用0,实际将使用eBufType*nBufXSize*nBufYSize来表示。

对于GDALDataset::RasterIO函数还有两个参数nBandCount和panBandMap,分别表示要读取的波段个数和波段序号，尤其是后一个参数波段序号，可以自定义先读取那一个波段，后读取那一个。具体使用方法见下。

三、RasterIO使用方法示例

使用示例一,在Windsow位图数据颜色排列是BGR，但是图像存储的可能是按照RGB来存储的，一般的做法是将数据按照每个波段读出来，然后再认为的按照BGR来进行组织，其实完全可以使用后面三个参数来将读出来的数据自动按照BGR的方式组织好。只要将参数设置为：nBandCount=3;panBandMap=new int[]{3,2,1}即可。

[cpp] view plaincopyprint?

int panBandMap [3]= {3,2,1}; //按照BGR BGR BGR … 来读取数据组织
DT_8U *pData = new DT_8U[iWidth*iHeight*3];
poDataset ->RasterIO(GF_Read, 0, 0, iWidth,iHeight, pData,iWidth, iHeight, (GDALDataType)iDataType, 3, panBandMap,iDataType*3, iDataType*iMthWidth*3, iDataType);

使用示例二，实现的是将7波段图像中的第2 3 4波段按照3 42的顺序读入内存中，逐像素存储：

[cpp] view plaincopyprint?

intbandmap[7];
bandmap[0]=3;
bandmap[1]=4;
bandmap[2]=2;
Scanline=(nImagSizex*8+31)/32*4;
BYTE pafScan=( BYTE )CPLMalloc(sizeof(byte)*nImgSizeX*nImgSizeY*3)
poDataset->RasterIO(GF_Read, 0, 0,nImgSizeX,nImgSizeY, pafScan,
ImgSizeX,nImgSizeY,GDT_Byte,3,bandmap,3,Scanline*3,1 );

以这种方式读取之后，直接可构建位图进行显示。这里可以按照自己的需要进行其他方式读取。以上读取方式仅仅为了显示方便，如进行图像处理相关运算，则按波段全部读出会比较方便，即按照常规的方式读取处理：

[cpp] view plaincopyprint?

poDataset->RasterIO( GF_Read, 0,0,nImgSizeX,nImgSizeY, pafScan,
nImgSizeX,nImgSizeY,GDT_Byte,bandcount,0,0,0,0);

之前申请的内存改为：

[cpp] view plaincopyprint?

BYTE pafScan=new byte[nImgSizeX*nImgSizeY*bandcount];

将图像数据读入内存后，即可通过指针pafScan对图像进行你想要进行的操作了。

以下内容更新于2012年5月25日

通过上面的介绍，发现很多人还是不太理解，下面举六个例子，基本上涵盖了RasterIO的所有用法：

1、第一种方式

常规的读图方式，这里只读取图像的第二波段的第二行数据，代码如下：

[cpp] view plaincopyprint?

//获取第二波段的波段指针，参数就是表示第几波段的意思
GDALRasterBand *pBand =poDataset->GetRasterBand(2);
//获取该波段的数据类型，如8U，16U等
GDALDataType dataType =pBand->GetRasterDataType();
//分配存储一行数据的空间，按照8U的数据类型
DT_8U *pBuf = newDT_8U[iWidth];
pBand->RasterIO(GF_Read,0, 1, iWidth, 1, pBuf, iWidth, 1, GDT_Byte, 0, 0);
delete []pBuf;
pBuf = NULL;

2、第二种方式

使用数据集读取，按照8bit的数据来进行处理，读取整幅图像，按照BGRBGR…BGR的顺序来组织，代码如下：

[cpp] view plaincopyprint?

DT_8U *pBuf = newDT_8U[iWidth*iWidth*iBandCount]; //分配存储空间
int panBandMap [3]={3,2,1}; //如果想读取为RGB，那么将数据换成1，2，3
poDataset->RasterIO(GF_Read, 0, 0, iWidth,iHeight, pBuf,iWidth, iHeight,
GDT_Byte, 3, panBandMap, 3, iWidth*3, 1);
delete []pBuf;
pBuf = NULL;

3、第三种方式

读取一个矩形区域的像元，这里读取第三波段的数据，读取范围为从100行，200列开始，到400行，400列结束的一个矩形区域，代码如下：

[cpp] view plaincopyprint?

//获取第三波段的波段指针，参数就是表示第几波段的意思
GDALRasterBand *pBand =poDataset->GetRasterBand(3);
//获取该波段的数据类型，如8U，16U等
GDALDataType dataType =pBand->GetRasterDataType();
//分配存储数据的空间，按照8U的数据类型，大小按照矩形区域来确定
DT_8U *pBuf = newDT_8U[200*300];
pBand->RasterIO(GF_Read,200, 100, 200, 300, pBuf, 200, 300, GDT_Byte, 0, 0);
delete []pBuf;
pBuf = NULL;

4、第四种方式

读取缩放区域，即读取一个区域的数据，然后把这个区域进行缩放，放大或者缩小，这个非常有用，比如在显示大图像的时候，就可以用这个方式来进行读取数据，这里依旧读取第三波段的数据，读取范围为从100行，200列开始，到400行，400列结束的一个矩形区域，按照0.5的缩放比例对区域进行缩放，所以原始图像的区域大小为200*300，缩放后就是100*150。代码如下：

[cpp] view plaincopyprint?

//为了举例，这里读取第三波段的数据，读取范围为从100行，200列开始，到400行，400列结束的一个矩形区域
GDALRasterBand *pBand =poDataset->GetRasterBand(3);
//获取该波段的数据类型，如8U，16U等
GDALDataType dataType =pBand->GetRasterDataType();
//按照0.5的区域进行缩放，所以原始图像的区域大小为200*300，缩放后就是100*150
//分配存储数据的空间，按照8U的数据类型，大小按照矩形区域来确定
DT_8U *pBuf = newDT_8U[100*150];
pBand->RasterIO(GF_Read,200, 100, 200, 300, pBuf, 100, 150, GDT_Byte, 0, 0);
delete []pBuf;
pBuf = NULL;

5、第五种方式

将8bit读取为16bit，为了方便说明，这里读取第一波段的第一行数据，这里假设打开的数据是一个8bit的数据。代码如下：

[cpp] view plaincopyprint?

GDALRasterBand *pBand =poDataset->GetRasterBand(1);
//获取该波段的数据类型，如8U，16U等
GDALDataType dataType =pBand->GetRasterDataType();
//分配存储一行数据的空间，按照16U的数据类型
DT_16U *pBuf = newDT_16U[iWidth];
pBand->RasterIO(GF_Read,0, 0, iWidth, 1, pBuf, iWidth, 1, GDT_UInt16, 0, 0);
delete []pBuf;
pBuf = NULL;

6、第六种方式

将16bit读取为8bit，为了方便说明，这里依旧读取第一波段的第一行数据，这里假设打开的数据是一个16bit的数据。在这里需要特别的注意，虽然GDAL提供这种可以把16bit的数据读取到8bit的存储空间中，GDAL会自动进行转换，但是这种转换是有损的，和把int值转为unsigned char类似，超过存储范围的就会被截取掉，这里也是一样，比如，16bit中的数据有个超过了255，那么使用8bit的数据存储来只剩下255了，所以一般不建议从这里转换/最好自己写个算法来进行转换，比如线性拉伸等，代码如下：

[cpp] view plaincopyprint?

//！！！！！！！！！！！！！！！！！！
//强烈注意，虽然GDAL提供这种可以把16bit的数据读取到8bit的存储空间中，但是会自动进行转换这种转换是有损的，和把int值转为unsigned char类似，超过存储范围的就会被截取掉，这里也是一样。比如，16bit中的数据有个超过了255，那么使用8bit的数据存储来只剩下255了，所以一般不建议从这里转换。最好自己写个算法来进行转换，比如线性拉伸等。
//！！！！！！！！！！！！！！！！！！
GDALRasterBand *pBand =poDataset->GetRasterBand(1);
GDALDataType dataType =pBand->GetRasterDataType();
DT_8U *pBuf = new DT_8U[iWidth];
pBand->RasterIO(GF_Read,0, 0, iWidth, 1, pBuf, iWidth, 1, GDT_Byte, 0, 0);
…
delete []pBuf;
pBuf = NULL;

使用GDAL的MEM内存文件保存临时文件

在使用GDAL编写算法的时候，经常会将计算的中间结果存在一个临时的图像文件中，然后使用完再将其删除，如果临时文件就一个的话，创建一个也无所谓，但是当一个复杂的算法中可能会出现很多个临时文件的时候（我在编写Hariss角点自动匹配算法的时候有4个临时文件），这种情况下总觉得临时文件很不爽，此外第一个不爽的地方；第二个图像太大的时候，临时文件也会占用很大的空间，假如空间不足或者给定的临时文件路径不可写等问题会让人头疼；第三，在创建临时文件的读写上会耗用比较多的时间，尤其在磁盘的IO时，耗时比较多。

基于上面的几个原因，总想不用临时文件，发现GDAL里面提供了一个VRT的虚拟文件（参考http://www.gdal.org/gdal_vrttut.html），查看GDAL的App中的gdalenhance.cpp文件中，发现VRT使用起来比较复杂，或者我不太熟悉而已，过段时间仔细研究一下这个东西。除了VRT，GDAL中还有一个叫MEM的内存文件，简介参考这个网页：http://gdal.org/frmt_mem.html。关于MEM文件的使用，和普通的文件使用一样，参考下面的代码：

[cpp] view plaincopyprint?

#include “gdal_priv.h”
int main()
{
GDALAllRegister();
//打开原始数据
GDALDataset *poSrcDS;
poSrcDS = (GDALDataset *) GDALOpen( “C:\\test.img”, GA_ReadOnly );
if( poSrcDS == NULL )
return -1;
//获取MEM的文件驱动
GDALDriver *poDriver;
poDriver = GetGDALDriverManager()->GetDriverByName(“MEM”);
if( poDriver == NULL )
return -1;
//这里或者使用Create函数创建一个MEM文件，文件路径给个空串就可以了
GDALDataset *poDstDS;
poDstDS = poDriver->CreateCopy( “”, poSrcDS, FALSE, NULL, NULL, NULL );
{
//do something
}
//使用完直接Close掉就好了，不用删除临时文件了
GDALClose( (GDALDatasetH) poDstDS );
GDALClose( (GDALDatasetH) poSrcDS );
return 0;
}

在调用GDAL算法的时候，希望能够显示其处理进度信息，其实在GDAL的算法API中，一般最后两个参数就是进度信息的指针。下面分别实现两种进度条信息，一种是在控制台中的进度条，一种是基于QT界面的进度条（你可以参考写一个MFC的）。

对于GDAL来说，本身就实现了一个基于控制台的进度条函数，名字叫GDALTermProgress，其函数说明参考这里，调用这个进度函数后，会在控制台中显示一个进度信息，形状大概就是下面的样子：

0…10…20…30…40…50…60…70…80…90…100- done.

每个2.5%就会输出一个点，然后每到10%就会输出数字，在完成结束后会输出done。

GDAL的算法中，一般最后两个参数是用来返回进度信息的，下面以栅格矢量化的接口举例说明，GDAL栅格矢量化的函数接口定义为：

CPLErr GDALPolygonize

(

GDALRasterBandH

hSrcBand,

GDALRasterBandH

hMaskBand,

OGRLayerH

hOutLayer,

int

iPixValField,

char **

papszOptions,

GDALProgressFunc

pfnProgress,

void *

pProgressArg

)

算法说明文档参考地址：http://www.gdal.org/gdal__alg_8h.html#a3f522a9035d3512b5d414fb4752671b1。从文档中可以看到，最后两个参数说明是：

pfnProgress

callback for reporting algorithm progress matching the GDALProgressFunc()semantics. May be NULL.

pProgressArg

callback argument passed to pfnProgress.

上面的意思大概就是说，参数pfnProgress 是一个回调函数指针，回调函数定义参考 GDALProgressFunc()，参数pProgressArg 是回调函数pfnProgress 的第三个参数。关于GDALProgressFunc()的定义如下：

int(*GDALProgressFunc)(double dfComplete, const char *pszMessage, void*pProgressArg);

该回调函数参数说明分别是：

dfComplete

进度值，取值范围是0.0~1.0，0.0表示开始，1.0表示完成。

pszMessage

可选项，用来显示的字符串信息。通常用NULL即可。

pProgressArg

应用程序回调函数参数，是一个自定义的类（或结构体）指针。

下面以GDAL自带的进度函数GDALTermProgress为例，进行说明。GDALTermProgress的源代码在GDALSrc\gcore\gdal_misc.cpp，762行左右（GDAL1.8.1版本），代码如下：

[cpp] view plaincopyprint?

int CPL_STDCALL GDALTermProgress( double dfComplete, const char *pszMessage,
void * pProgressArg )
{
static int nLastTick = -1;
int nThisTick = (int) (dfComplete * 40.0);
(void) pProgressArg;
nThisTick = MIN(40,MAX(0,nThisTick));
// Have we started a new progress run?
if( nThisTick < nLastTick && nLastTick >= 39 )
nLastTick = -1;
if( nThisTick <= nLastTick )
return TRUE;
while( nThisTick > nLastTick )
{
nLastTick++;
if( nLastTick % 4 == 0 )
fprintf( stdout, “%d”, (nLastTick / 4) * 10 );
else
fprintf( stdout, “.” );
}
if( nThisTick == 40 )
fprintf( stdout, ” – done.\n” );
else
fflush( stdout );
return TRUE;
}

可以看出，上面的函数中，只用了第一个参数，后面两个参数没用，就是将第一个参数使用printf函数输出在屏幕上。基于上面的分析，可以自己实现一个自己的GDAL进度条类，首先看类声明：

[cpp] view plaincopyprint?

/**
* @brief 进度条基类
*
* 提供进度条基类接口，来反映当前算法的进度值
*/
class IMGALG_API CProcessBase
{
public:
/**
* @brief 构造函数
*/
CProcessBase()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
m_bIsContinue = true;
}
/**
* @brief 析构函数
*/
virtual ~CProcessBase() {}
/**
* @brief 设置进度信息
* @param pszMsg 进度信息
*/
virtual void SetMessage(const char* pszMsg) = 0;
/**
* @brief 设置进度值
* @param dPosition 进度值
* @return 返回是否取消的状态，true为不取消，false为取消
*/
virtual bool SetPosition(double dPosition) = 0;
/**
* @brief 进度条前进一步，返回true表示继续，false表示取消
* @return 返回是否取消的状态，true为不取消，false为取消
*/
virtual bool StepIt() = 0;
/**
* @brief 设置进度个数
* @param iStepCount 进度个数
*/
virtual void SetStepCount(int iStepCount)
{
ReSetProcess();
m_iStepCount = iStepCount;
}
/**
* @brief 获取进度信息
* @return 返回当前进度信息
*/
string GetMessage()
{
return m_strMessage;
}
/**
* @brief 获取进度值
* @return 返回当前进度值
*/
double GetPosition()
{
return m_dPosition;
}
/**
* @brief 重置进度条
*/
void ReSetProcess()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
m_bIsContinue = true;
}
/*! 进度信息 */
string m_strMessage;
/*! 进度值 */
double m_dPosition;
/*! 进度个数 */
int m_iStepCount;
/*! 进度当前个数 */
int m_iCurStep;
/*! 是否取消，值为false时表示计算取消 */
bool m_bIsContinue;
};

该类是一个进度条基类，将一些通用的函数进行实现，比如获取进度值，获取进度信息等之类的，然后将一些需要根据使用方式不同的设置为虚函数，比如设置进度值等。下面给出一个控制台进度条类，用于显示基于控制台程序的进度值显示，参考GDALTermProgress函数写的：

[cpp] view plaincopyprint?

/**
* @brief 控制台进度条类
*
* 提供控制台程序的进度条类接口，来反映当前算法的进度值
*/
class CConsoleProcess : public CProcessBase
{
public:
/**
* @brief 构造函数
*/
CConsoleProcess()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
};
/**
* @brief 析构函数
*/
~CConsoleProcess()
{
//remove(m_pszFile);
};
/**
* @brief 设置进度信息
* @param pszMsg 进度信息
*/
void SetMessage(const char* pszMsg)
{
m_strMessage = pszMsg;
printf(“%s\n”, pszMsg);
}
/**
* @brief 设置进度值
* @param dPosition 进度值
* @return 返回是否取消的状态，true为不取消，false为取消
*/
bool SetPosition(double dPosition)
{
m_dPosition = dPosition;
TermProgress(m_dPosition);
m_bIsContinue = true;
return true;
}
/**
* @brief 进度条前进一步
* @return 返回是否取消的状态，true为不取消，false为取消
*/
bool StepIt()
{
m_iCurStep ++;
m_dPosition = m_iCurStep*1.0 / m_iStepCount;
TermProgress(m_dPosition);
m_bIsContinue = true;
return true;
}
private:
void TermProgress(double dfComplete)
{
static int nLastTick = -1;
int nThisTick = (int) (dfComplete * 40.0);
nThisTick = MIN(40,MAX(0,nThisTick));
// Have we started a new progress run?
if( nThisTick < nLastTick && nLastTick >= 39 )
nLastTick = -1;
if( nThisTick <= nLastTick )
return ;
while( nThisTick > nLastTick )
{
nLastTick++;
if( nLastTick % 4 == 0 )
fprintf( stdout, “%d”, (nLastTick / 4) * 10 );
else
fprintf( stdout, “.” );
}
if( nThisTick == 40 )
fprintf( stdout, ” – done.\n” );
else
fflush( stdout );
}
};

至此，一个控制台的进度条就完成了，使用方式如下（注意：之前的博客中出现的LT_ConsoleProgress就是上面这个类）：

[cpp] view plaincopyprint?

void main()
{
CConsoleProcess *pProgress = new CConsoleProcess();
//参考建立Erdas金字塔的那篇博客：http://blog.csdn.net/liminlu0314/article/details/6127755
int f = CreatePyramids(“C://Work//Data//ttttt.img”, pProgress);
if (f == RE_SUCCESS)
printf(“计算成功/n”);
else
printf(“计算失败/n”);
delete pProgress;
}

效果图如下：

下面开始编写一个基于QT界面的进度条类，首先是类定义，基于QT界面的，对于MFC界面的可以参考这个写一个：

[cpp] view plaincopyprint?

/**
* @brief 进度条对话框类
*
* 提供GUI程序的进度条类接口，来反映当前算法的进度值
*/
class CProcessDlg :
public QProgressDialog,
public CProcessBase
{
Q_OBJECT
public:
/**
* @brief 构造函数
*/
CProcessDlg(QWidget *parent = 0);
/**
* @brief 析构函数
*/
~CProcessDlg();
/**
* @brief 设置进度信息
* @param pszMsg 进度信息
*/
void SetMessage(const char* pszMsg);
/**
* @brief 设置进度值
* @param dPosition 进度值
*/
bool SetPosition(double dPosition);
/**
* @brief 进度条前进一步
*/
bool StepIt();
public slots:
void updateProgress(int);
};

接下来是实现部分代码，具体和上面的控制台类一样，只不过是设置进度值都是设置在界面上的相应的控件中：

[cpp] view plaincopyprint?

/**
* @brief 构造函数
*/
CProcessDlg::CProcessDlg(QWidget *parent)
:QProgressDialog(parent)
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
setModal(true);
setLabelText(tr(“处理中…”));
setAutoClose(false);
setAutoReset(false);
setCancelButtonText(tr(“取消”));
setWindowTitle(tr(“进度”));
//禁用关闭按钮
setWindowFlags(Qt::WindowTitleHint | Qt::FramelessWindowHint | Qt::WindowStaysOnTopHint);
};
/**
* @brief 析构函数
*/
CProcessDlg::~CProcessDlg()
{
};
/**
* @brief 设置进度信息
* @param pszMsg 进度信息
*/
void CProcessDlg::SetMessage(const char* pszMsg)
{
if (pszMsg != NULL)
{
m_strMessage = pszMsg;
setLabelText(QString(pszMsg));
}
}
/**
* @brief 设置进度值
* @param dPosition 进度值
*/
bool CProcessDlg::SetPosition(double dPosition)
{
m_dPosition = dPosition;
setValue( std::min( 100u, ( uint )( m_dPosition*100.0 ) ) );
QCoreApplication::instance()->processEvents();
if(this->wasCanceled())
return false;
return true;
}
/**
* @brief 进度条前进一步,返回false表示终止操作
*/
bool CProcessDlg::StepIt()
{
m_iCurStep ++;
m_dPosition = m_iCurStep*1.0 / m_iStepCount;
setValue( std::min( 100u, ( uint )( m_dPosition*100.0 ) ) );
QCoreApplication::instance()->processEvents();
if(this->wasCanceled())
return false;
return true;
}
void CProcessDlg::updateProgress(int step)
{
this->setValue(step);
QCoreApplication::instance()->processEvents();
}

至此，QT界面的进度条类已经完成，调用方式和上面控制台的基本类似，示例代码：

[cpp] view plaincopyprint?

//省略部分代码
CProcessDlg *pPro = new CProcessDlg();
pPro->setWindowTitle(tr(“正在执行栅格矢量化”));
pPro->show();
int iRev = ImagePolygonize(m_strInputName.c_str(), m_strOutputName.c_str(), pszFormat, pPro);
delete pPro;

这样，我们基本上就实现了两个进度条，一个控制台，一个QT界面的。执行的效果图如下：

那么在自己的算法中如何使用这个进度条类呢，下面给出一个简单的算法，图像反色算法：

[cpp] view plaincopyprint?

/**
* @brief 图像反色
* @param pszSrcFile 输入文件路径
* @param pszDstFile 输出文件路径
* @param pszFormat 输出文件格式，详细参考GDAL支持数据类型
* @param pProcess 进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageAnticolor(const char* pszSrcFile, const char* pszDstFile, const char* pszFormat, CProcessBase* pProcess = NULL)
{
if(pProcess != NULL)
{
pProcess->ReSetProcess();
pProcess->SetMessage(“执行图像反色…”);
}
if(pszSrcFile == NULL || pszDstFile == NULL)
return RE_PARAMERROR;
GDALAllRegister();
GDALDataset *poSrcDS = (GDALDataset *) GDALOpen( pszSrcFile, GA_ReadOnly );
if( poSrcDS == NULL )
{
if(pProcess != NULL)
pProcess->SetMessage(“输入文件不能打开，请检查文件是否存在！”);
return RE_FILENOTEXIST;
}
GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if( poDriver == NULL )
{
if(pProcess != NULL)
pProcess->SetMessage(“不能创建制定类型的文件，请检查该文件类型GDAL是否支持创建！”);
GDALClose( (GDALDatasetH) poSrcDS );
return RE_FILENOTSUPPORT;
}
//获取图像宽高和波段数
int iXSize = poSrcDS->GetRasterXSize();
int iYSize = poSrcDS->GetRasterYSize();
int iBandCount = poSrcDS->GetRasterCount();
//假设输入图像也是8U的数据
GDALDataset *poDstDS = poDriver->Create(pszDstFile, iXSize, iYSize, iBandCount, GDT_Byte, NULL);
double dGeoTrans[6] = {0};
//设置仿射变换参数
poSrcDS->GetGeoTransform(dGeoTrans);
poDstDS->SetGeoTransform(dGeoTrans);
//设置图像投影信息
poDstDS->SetProjection(poSrcDS->GetProjectionRef());
DT_8U *pSrcData = new DT_8U[iXSize];
DT_8U *pDstData = new DT_8U[iXSize];
if(m_pProcess != NULL)
{
m_pProcess->SetMessage(“计算图像反色…”);
m_pProcess->SetStepCount(iYSize*iBandCount); //设置进度条的总步长
}
//循环波段
for(int iBand=1; iBand<=iBandCount; iBand++)
{
GDALRasterBand *pSrcBand = poSrcDS->GetRasterBand(iBand);
GDALRasterBand *pDstBand = poDstDS->GetRasterBand(iBand);
for(int i=0; i<iYSize; i++) //循环图像高
{
pSrcBand->RasterIO(GF_Read, 0, i, iXSize, 1, pSrcData, iXSize, 1, GDT_Byte, 0, 0);
for(int j=0; j<iXSize; j++) //循环图像宽
pDstData[j] = 255 – pSrcData[j];
pDstBand->RasterIO(GF_Write, 0, i, iXSize, 1, pDstData, iXSize, 1, GDT_Byte, 0, 0);
if(m_pProcess != NULL)
{
bool bIsCancel = m_pProcess->StepIt();
if(!bIsCancel)
{
RELEASE(pSrcData);
RELEASE(pDstData);
//关闭原始图像和结果图像
GDALClose( (GDALDatasetH) poDstDS );
GDALClose( (GDALDatasetH) poSrcDS );
remove(pszDstFile); //删除结果图像
if(pProcess != NULL)
pProcess->SetMessage(“图像反色取消！”);
return RE_CANCEL;
}
}
}
}
RELEASE(pSrcData);
RELEASE(pDstData);
//关闭原始图像和结果图像
GDALClose( (GDALDatasetH) poDstDS );
GDALClose( (GDALDatasetH) poSrcDS );
if(pProcess != NULL)
pProcess->SetMessage(“图像反色完成！”);
return RE_SUCCESS;
}

如果我要使用GDAL提供的算法，该怎么使用上面的进度条呢，很简单，具体是要实现一个GDALProgressFunc()函数相同的一个使用__stdcall导出的函数即可，函数定义如下：

[cpp] view plaincopyprint?

/**
* @brief 导出符号定义
*/
#ifndef STD_API
#define STD_API __stdcall
#endif
/**
* \brief 调用GDAL进度条接口
*
* 该函数用于将GDAL算法中的进度信息导出到CProcessBase基类中，供给界面显示
*
* @param dfComplete 完成进度值，其取值为 0.0 到 1.0 之间
* @param pszMessage 进度信息
* @param pProgressArg CProcessBase的指针
*
* @return 返回TRUE表示继续计算，否则为取消
*/
int STD_API ALGTermProgress( double dfComplete, const char *pszMessage, void * pProgressArg );

下面为该函数的实现代码，请重点看第三个参数的使用方式：

[cpp] view plaincopyprint?

/**
* \brief 调用GDAL进度条接口
*
* 该函数用于将GDAL算法中的进度信息导出到CProcessBase基类中，供给界面显示
*
* @param dfComplete 完成进度值，其取值为 0.0 到 1.0 之间
* @param pszMessage 进度信息
* @param pProgressArg CProcessBase的指针
*
* @return 返回TRUE表示继续计算，否则为取消
*/
int STD_API ALGTermProgress( double dfComplete, const char *pszMessage, void * pProgressArg )
{
if(pProgressArg != NULL)
{
CProcessBase * pProcess = (CProcessBase*) pProgressArg;
pProcess->m_bIsContinue = pProcess->SetPosition(dfComplete);
if(pProcess->m_bIsContinue)
return TRUE;
else
return FALSE;
}
else
return TRUE;
}

稍微对上面的函数体进行说明一下，关键是第三个参数，第三个参数其实就是我们上面定义的控制台类获取QT界面类的对象指针，只不过是专为一个void指针，在函数体中，我们需要将其转为基类CProcessBase的指针，然后调用该基类的设置进度值函数即可。对于该函数的使用，下面还是以栅格矢量化为例来进行说明，废话不多说，看代码：

[cpp] view plaincopyprint?

//前面省略很多代码
GDALRasterBandH hSrcBand = (GDALRasterBandH)poSrcDS->GetRasterBand(1);
if(GDALPolygonize(hSrcBand, NULL, (OGRLayerH)poLayer, 0, NULL, ALGTermProgress, pProcess)!= CE_None)
{
GDALClose( (GDALDatasetH) poSrcDS );
if(pProcess != NULL)
{
if(!pProcess->m_bIsContinue)
{
pProcess->SetMessage(“计算取消！”);
return RE_USERCANCEL;
}
else
{
pProcess->SetMessage(“计算失败！”);
return RE_PARAMERROR;
}
}
return RE_PARAMERROR;
}
//后面省略部分代码

在用到GDAL时，经常会用到Proj4和GEOS，关于这两个库的作用，可以到其官网看看。下面编译是在Windows环境下，编译器使用MS的VS2008。

一、编译PROJ4

PROJ4的最新版本是4.8，官网地址为：http://trac.osgeo.org/proj/。从官网下载PROJ4的源代码，解压到文件夹中，如F:\Work\3rdPart\proj-4.8.0。

1、正常编译RELEASE版本

打开VS2008的命令行工具，然后将其工作目录切换到F:\Work\3rdPart\proj-4.8.0，如下图所示：

如果不进行输出目录设置的话，就直接在命令行中依次输入下面的命令回车即可：

nmake /f makefile.vc install-all

等待编译完成后，会默认值C盘的根目录下，创建PROJ文件夹，里面有四个文件夹，分别是bin，lib，include以及share四个文件夹，其中include和lib是用来做二次开发使用，bin存放的是dll和exe文件，share里面存储的是PROJ4所定义的一些投影文件等，在发布程序的时候，share文件夹需要一同进行发布，否则在做投影转换的时候可能因为找不到其中的文件而导致转换失败。

2、编译DEBUG版本

在有的时候需要调试PROJ4的源代码，那么需要编译DEBUG版本，编译DEBUG版本和RELEASE版本一样，只不过在是最后输入命令的时候，在后面加上DEBUG=1即可，完整命令如下：

nmake /f makefile.vc clean

nmake /f makefile.vc install-all DEBUG=1

等编译结束后，将src目录下的pdb等调试文件拷贝到你自己的工程输出目录中即可。nmake /f makefile.vc clean，这句的目的是为了清理之前编译生成的临时文件，如果之前没有编译过，可以不用。

3、编译X64版本

有时候需要在64位系统上运行，为了高效，需要编译X64的版本，编译X64的版本和上面的基本一样，只不过是在打开VS2008的命令行的时候，要使用X64兼容工具命令提示，如下图所示：

二、编译GEOS

目前GEOS的最新版本是3.3.2，官网地址是：http://trac.osgeo.org/geos/。从官网下载GEOS的源代码，解压到文件夹中，如F:\Work\3rdPart\geos-3.3.2。

1、正常编译RELEASE X86版本

编译过程和上面的PROJ4基本一致，只不过在执行nmake之前，有点小区别。打开VS2008的命令行工具，然后将其工作目录切换到F:\Work\3rdPart\geos-3.3.2，首先要执行一下，autogen,bat，然后再进行编译，如下图所示：

命令如下：

atuogen.bat

nmake /f makefile.vc src_dir

注意，上面截图的时候，后面的命令行敲错了，之前3.2.2版本是source_dir，在3.3.2版本，改为了src_dir。在编译3.3.2版本时，提示一个错误，log函数参数不匹配，用记事本打开F:\Work\3rdPart\geos-3.3.2\src\operation\buffer\BufferOp.cpp，找到第97行处的log(10)，将log(10)改为log(10.0)，保存，然后重新nmake即可。修改后的代码如下，修改的部分在第17行标注：

[cpp] view plaincopyprint?

<span style=”color:#000099;”>/*private*/
double
BufferOp::precisionScaleFactor(const Geometry *g,
double distance,
int maxPrecisionDigits)
{
const Envelope *env=g->getEnvelopeInternal();
double envMax = std::max(
std::max(fabs(env->getMaxX()), fabs(env->getMinX())),
std::max(fabs(env->getMaxY()), fabs(env->getMinY()))
);
double expandByDistance = distance > 0.0 ? distance : 0.0;
double bufEnvMax = envMax + 2 * expandByDistance;
// the smallest power of 10 greater than the buffer envelope
int bufEnvPrecisionDigits = (int) (std::log(bufEnvMax) / std::log(10.0) + 1.0);
int minUnitLog10 = maxPrecisionDigits – bufEnvPrecisionDigits;
double scaleFactor = std::pow(10.0, minUnitLog10);
return scaleFactor;
}</span>

2、编译DEBUG版本

在编译DEBUG版本和RELEASE版本一样，只不过在是最后输入命令的时候，在后面加上DEBUG=1即可，完整命令如下：

atuogen.bat

nmake /f makefile.vc clean

nmake /f makefile.vc src_dir DEBUG=1

等编译结束后，将src目录下的pdb等调试文件拷贝到你自己的工程输出目录中即可。同样，nmake /f makefile.vc clean，这句的目的是为了清理之前编译生成的临时文件，如果之前没有编译过，可以不用。

3、编译X64版本

编译X64的版本和编译PROJ4的X64一样，只不过是输入的命令不同而已。

GDAL源码剖析（九）之GDAL体系架构

在GDAL库中包含栅格数据的读写，矢量数据的读写，以及栅格和矢量数据的相关算法。下面主要对GDAL中栅格数据和矢量数据的体系架构做一个简单的说明。本人英文很烂，有些部分写出来的东西自己都看不懂，如果不懂，可以看英文。

一、GDAL体系架构

参考GDAL官方文档：http://www.gdal.org/gdal_datamodel.html。

GDAL使用抽象数据模型(abstract data model)来解析它所支持的数据格式，抽象数据模型包括数据集(dataset)，坐标系统（Coordinate System），仿射地理坐标转换(Affine GeoTransform)，大地控制点(GCPs)，元数据(Metadata)，子数据集域(Subdatasets Domain)，图像结构域(Image_Structure Domain)，RPC域（RPC Domain），XML域(XML:Domains)，栅格波段(Raster
Band)，颜色表(Color Table)和快视图（Overviews）。

1、数据集（Dataset）

数据集是由栅格波段和一些相关的信息共同组成的（具体参考GDALDataset类），在一个指定的数据集中所有的波段都具有相同的大小（图像行数和列数）；数据集同时负责对所有的波段定义地理坐标和投影等信息；数据集中还有一些元数据信息，这些信息是由字符串组成的一系列（名称：值）的字符串列表。

GDAL数据集是基于OpenGIS网格数据的说明来实现的。

2、坐标系统（CoordinateSystem）

数据集中的坐标系统是基于OpenGIS的WKT（WellKnown Text）字符串格式的。里面包含下面内容：

坐标系统名称
地理坐标系统名称
大地水准面
托球体名称，长半轴和扁率
中央经线名称和与本初子午线的偏移量
投影方式（如：横轴墨卡托）
投影参数列表（如：中央经线）
单位名称和与米或者弧度之间的转换系数
名称和坐标轴顺序
以上项对应的EPSG代码

更多关于OpenGIS的WKT坐标定义和相关信息，可以参考文档osr_tutorial。获取坐标系统可以使用函数 GDALDataset::GetProjectionRef()，如果该函数的返回值为””，说明该图像没有地理坐标系统。

1、仿射地理坐标转换（Affine GeoTransform）

GDAL数据集有两种方式来表示栅格行列号坐标和地理坐标之间的关系，第一种，也是最常用的，就是使用仿射变换来表示（另外一种是GCP点）。

仿射变换包含六个参数，可以使用函数 GDALDataset::GetGeoTransform()获取，图像行列号和地理空间之间的变换关系如下：

Xgeo = GT(0) +Xpixel*GT(1) + Yline*GT(2)

Ygeo = GT(3) +Xpixel*GT(4) + Yline*GT(5)

在北方向上的图像中，参数GT(2)和GT(4)的值是0，GT(1)表示象元宽度，GT(5)表示象元高度。点 (GT(0),GT(3))表示图像左上角的横纵坐标值。

值得注意的是：行列号坐标系统中，图像的（0.0，0.0）坐标在图像的左上角，坐标（象元宽个数,象元高个数）表示图像的右下角坐标。图像的左上角象元的中心点坐标是（0.5,0.5）。

2、大地控制点（GCPs）

数据集中可能会包括一系列的控制点，用来确定图像的地理参考坐标。所有的控制点都是在同一个坐标系统下的（可以使用函数GDALDataset::GetGCPProjection()来获取该坐标系），每一个GCP（使用GDAL_GCP类来表示）包括下面的内容：

[cpp] view plaincopyprint?

1. typedef struct

2. {

3. char *pszId;

4. char *pszInfo;

5. double dfGCPPixel;

6. double dfGCPLine;

7. double dfGCPX;

8. double dfGCPY;

9. double dfGCPZ;

10.} GDAL_GCP;

pszId字符串表示一个唯一的字符串，用来藐视在当前数据集中的GCP点中的具体的GCP点（通常情况下，可能是由数字构成，但不是所有的情况下）。pszInfo通常是一个空字符串，但是可以用来存储用户自定义用来描述GCP的文本信息。很可能包含GCP的一些状态信息，如获取时间等。

坐标（Pixel,Line）表示GCP在图像上的位置，坐标（X,Y,Z）表示在地理参考坐标系中的位置，通常Z值是0。

GDAL数据模型中不会都包含GCP点构成的转换方程，这个工作常常在应用程序中处理。通常使用1次到5次多项式进行转换。

通常数据集包括一个仿射地理坐标变换，GCP点可能会有，有时候这两个都没有。

3、元数据（Metadata）

GDAL元数据是辅助格式和应用程序指定的数据使用一个名称/值的点对列表。名称需要能够很好的表示其作用（不含有空格或者其他特殊字符）。值可以有任意长度，可以包括任何东西除了空值（ASCII的0值）。

元数据在系统中不能处理太大的数据，如果元数据超过100kb就会导致性能下降。

一些常用的的格式的元数据都可以支持，比如在TIFF驱动返回一些了标签信息作为元数据，比如日期/时间的形式如下：

TIFFTAG_DATETIME=1999:05:11 11:29:56

元数据可以按照名称组来拆分为不同的域，默认域是没有名称的，如（NULL或者””）。一些特殊的域用于特殊的目的。

下面的元数据项目定义在默认的域中：

AREA_OR_POINT: 可以是“Area”或者“Point”，默认为Area。表明一个像素值表示中心像元的像素值或者像元中心点的采样值。但这并不影响整个区域的地理参考信息。
NODATA_VALUES:NODATA值。这个值是使用空格隔开的一系列的像素值，个数和波段数相同。使用这种方式是为了使所有波段都有一个NODATA值相匹配。这个元数据没有广泛的用于数据驱动，算法或者工具集中。
MATRIX_REPRESENTATION:这个值用来表示极化SAR数据集，包含矩阵表示的数据。矩阵可以用下面的值来进行表示：
- SCATTERING 散射
- SYMMETRIZED_SCATTERING 对称散射
- COVARIANCE 协方差
- SYMMETRIZED_COVARIANCE 对称协方差
- COHERENCY 相干
- SYMMETRIZED_COHERENCY 对称相干
- KENNAUGH
- SYMMETRIZED_KENNAUGH
POLARMETRIC_INTERP:这个元数据项用来定义多波段极化SAR数据的栅格波段。指定的矩阵表示这个波段的数据。数据集就是一个散射矩阵，比如，这个元数据可接受的值有，HH、HV、VH、VV四个。当数据集是一个协方差矩阵时，这个元数据将分别是Covariance_11、Covariance_22、Covariance_33、Covariance_12、Covariance_13、Covariance_23中的一个(因为这个矩阵本身就是一个Hermitian矩阵，即所有的数据都用来表示这个矩阵)。

4、子数据集域（Subdatasets Domain）

子数据集域中含有一个子数据集的列表。通常用来存储在一个文件的多个文件（如HDF和NITF格式），例如，一个NITF格式的文件中获取子数据集的信息如下：

SUBDATASET_1_NAME=NITF_IM:0:multi_1b.ntf

SUBDATASET_1_DESC=Image 1 of multi_1b.ntf

SUBDATASET_2_NAME=NITF_IM:1:multi_1b.ntf

SUBDATASET_2_DESC=Image 2 of multi_1b.ntf

SUBDATASET_3_NAME=NITF_IM:2:multi_1b.ntf

SUBDATASET_3_DESC=Image 3 of multi_1b.ntf

SUBDATASET_4_NAME=NITF_IM:3:multi_1b.ntf

SUBDATASET_4_DESC=Image 4 of multi_1b.ntf

SUBDATASET_5_NAME=NITF_IM:4:multi_1b.ntf

SUBDATASET_5_DESC=Image 5 of multi_1b.ntf

含有_NAME的字符串，可以使用GDALOpen()函数来进行访问子文件。含有_DESC字符串是用来方便的展示给用户选择哪个文件。

5、图像结构域（Image_StructureDomain）

在默认域的元数据用来描述图像的相关信息，这些元数据并没有以特殊的方式保存在磁盘中。意思就是说，当这些元数据信息被复制到一个新的图像中的时候，会以合适的方式存在新的图像中。此外还有一些信息是和图像的格式紧密关联的。为了防止这些信息在写入新图像的时候被复制，将这些信息存储在一个叫IMAGE_STRUCTURE的特殊域中，通常这里存储的东西不会写入到新的图像中去。

目标，在IMAGE_STRUCTURE域中定义的项目主要有下面几个。

COMPRESSION:压缩方式用来指定数据集或者波段的压缩方式，这个是没有固定的压缩类型名称，但是一个指定的格式，如果指定压缩方式，那么在创建的时候会使用这种压缩方式进行压缩。
NBITS:当前波段或者当前数据集的波段中真实的数据存储bit数。通常只有非标准的数据类型才会使用该值，比如TIFF图像中的1bit数据，在GDAL会使用GDT_Byte来表示。
INTERLEAVE:这个只能用在数据集中，用来表示一个像元，一行和一个波段之间的间隔，可以用来作为读取数据的一个提示。
PIXELTYPE:这个值会出现在一个GDT_Byte类型的波段（或者相应的数据集）中，并且会使用SIGNEDBYTE值来表示无符号byte值介于128和255之间值应该转换为-128到-1之间。

6、RPC域（RPCDomain）

RPC元数据域存储的是有理函数模型的系数，该模型表示从图像行列号与空间参考位置间的变换关系，该模型具体定义如下：

· ERR_BIAS: 偏移误差，图像上所有的点在水平轴上的偏移的中误差，-1.0表示未知。

· ERR_RAND: 随机误差，图像上所有的点在水平轴上的随机中误差，-1.0表示未知。

· LINE_OFF: 行偏移量

· SAMP_OFF: 列偏移量

· LAT_OFF: 纬度偏移量

· LONG_OFF: 经度偏移量

· HEIGHT_OFF: 高程偏移量

· LINE_SCALE: 行缩放比例

· SAMP_SCALE: 列缩放比例

· LAT_SCALE: 纬度缩放比例

· LONG_SCALE: 经度缩放比例

· HEIGHT_SCALE: 高程缩放比例

· LINE_NUM_COEFF(1-20): 行分子系数，一共20个（使用空格隔开）

· LINE_DEN_COEFF(1-20): 行分母系数，一共20个（使用空格隔开）

· SAMP_NUM_COEFF(1-20): 列分子系数，一共20个（使用空格隔开）

· SAMP_DEN_COEFF(1-20): 列分母系数，一共20个（使用空格隔开）

这些内容都是来自GeoTIFF的RPC文档(http://geotiff.maptools.org/rpc_prop.html)。

7、XML域（XML:Domains）

任何前缀是“xml:”的一个字符串，但不是名称/值这种类型的，这是一个简单的XML格式的长字符串。

8、栅格波段（Raster Band）

一个栅格波段在GDAL中使用GDALRasterBand类来进行表示。他表示一个栅格波段、通道或者图层。波段不能完全用来表示整个图像，比如一个24位的RGB图像中就含有三个波段，分别是红波段，绿波段和蓝波段。

栅格波段含有下面属性：

· 图像的宽和高，这个和数据集里面的定义一样，如果这个波段是全分辨率波段的话。（这里有个说明，GDALRasterBand还可以表示金字塔的波段，如果是金字塔的波段的话，里面的宽高就和图像的宽高不一样）。

· 数据类型(GDALDataType)。应该是Byte、UInt16、Int16、UInt32、Int32、Float32、Float64以及复数类型CInt16、CInt32、CFloat32和CFloat64中的一个。

· 块大小。通过块是读取数据最高效的方式，对于分块数据，就是一个分块大小，对于大多数图像来说，一块就是一行。

· 名称/值的元数据对，格式和数据集中的一员，但是包含的信息可能是波段特有的。

· 一个可选的波段描述字符串。

· 一个可选的用来描述NODATA值的像元值。

· 一个可选的NODATA值表示的掩码波段或者在某些时候作为透明通道。

· 可选的类别名称列表（用于分类图）

· 可选的最大值和最小值。

· 可选的偏移量和缩放比例，用来对图像的像素值进行变换，比如变换高度到米等。

· 图像单位名称，可选。比如可以用来表示高程数据的海拔。

· 波段的颜色信息，是下面值中的某一个：

o GCI_Undefined:默认值，未知

o GCI_GrayIndex: 灰度图像

o GCI_PaletteIndex:颜色表图像

o GCI_RedBand: RGBA图像的R部分

o GCI_GreenBand: RGBA图像的G部分

o GCI_BlueBand: RGBA图像的B部分

o GCI_AlphaBand: RGBA图像的Alpha部分

o GCI_HueBand: HLS图像的色调部分

o GCI_SaturationBand: HLS图像的饱和度部分

o GCI_LightnessBand:HLS图像的亮度部分

o GCI_CyanBand: CMYK图像的青色部分

o GCI_MagentaBand: CMYK图像的品红部分

o GCI_YellowBand: CMYK图像的黄色部分

o GCI_BlackBand: CMYK图像的黑色部。

· 颜色表，下面有更详细的说明。

· 如果金字塔可用，含有一些关于金字塔的信息。

9、颜色表（Color Table）

颜色表的定义如下，使用C语言的风格定义：

[cpp] view plaincopyprint?

typedef struct
{
/- gray, red, cyan or hue-/
short c1;
/- green, magenta, orlightness -/
short c2;
/- blue, yellow, orsaturation -/
short c3;
/- alpha or blackband -/
short c4;
} GDALColorEntry;

颜色表通常是颜色调色板的一个值，对于c1/c2/c3/c4 四个值对应不同的调色板，其表示的含义不同，具体表示见下：

l GPI_Gray: 使用c1表示灰度值。

l GPI_RGB: c1表示红色，c2表示绿色，c3表示蓝色，c4表示alpha通道。

l GPI_CMYK: c1表示青色，c2表示洋红，c3表示黄色，c4表示黑色。

l GPI_HLS: c1表示色调，c2表示亮度，c3表示饱和度。

通过颜色表，将象元值用颜色表中的颜色来进行表示，颜色表中的值是从0开始递增。

10、快视图（Overviews）

一个波段中可能没有或者有很多个快视图。每个快视图都是一个GDALRasterBand，略缩图大小将和其下层的栅格大小不一样，但是快视图表示的区域与整个图像的区域是一致的。

快视图是用来快速显示图像用的，使用全分辨率图像进行降采样得到。

波段可以通过函数HasArbitraryOverviews，返回TRUE表示可以快速读取任何分辨率的快视图，这个主要应用在一些如FFT编码的图像（这句有点别扭，自己看英文原文吧）。

11、GDAL核心类结构设计

GDAL的核心类结构设计如图所示：

其中的类说明如下：

GDALMajorObject类：带有元数据的对象。

GDALDdataset类：通常是从一个栅格文件中提取的相关联的栅格波段集合和这些波段的元数据；GDALDdataset也负责所有栅格波段的地理坐标转换(georeferencingtransform)和坐标系定义。

GDALDriver类：文件格式驱动类，GDAL会为每一个所支持的文件格式创建一个该类的实体，来管理该文件格式。

GDALDriverManager类：文件格式驱动管理类，用来管理GDALDriver类。

二、OGR体系架构

1、OGR体系结构

OGR包括如下几部分：

Geometry：类Geometry (包括OGRGeometry等类)封装了OpenGIS的矢量数据模型，并提供了一些几何操作，WKB(Well Knows Binary)和WKT(Well Known Text)格式之间的相互转换，以及空间参考系统(投影)。

Spatial Reference：类OGRSpatialReference封装了投影和基准面的定义。

Feature：类OGRFeature封装了一个完整feature的定义，一个完整的feature包括一个geometry和geometry的一系列属性。

Feature Definition：类OGRFeatureDefn里面封装了feature的属性，类型、名称及其默认的空间参考系统等。一个OGRFeatureDefn对象通常与一个层(layer)对应。

Layer：类OGRLayer是一个抽象基类，表示数据源类OGRDataSource里面的一层要素(feature)。

Data Source：类OGRDataSource是一个抽象基类，表示含有OGRLayer对象的一个文件或一个数据库。

Drivers：类OGRSFDriver对应于每一个所支持的矢量文件格式。类OGRSFDriver由类OGRSFDriverRegistrar来注册和管理。

2、OGR的Geometry模型

OGR的Geometry模型是建立在OpenGIS的简单要素数据模型之上的。如下图所示：

图-OGR的Geometry模型关系图

OpenGIS的简单要素数据模型，其关系图如下所示：

图-OpenGIS的简单要素数据模型

由上面两图的对比，可以清楚的看到，OGR的Geometry模型是严格遵循OpenGIS的简单要素数据规范的。OGR的Geometry模型不仅在继承体系上与OpenGIS的简单要素数据模型一致，在函数接口上也向其靠拢，从基本的获取Geometry对象信息的方法如Dimension ()、GeometryType ()、SRID ()、Envelope()、AsText()、Boundary()等到判定空间未知关系的方法如Equals(anotherGeometry:Geometry)、Disjoint(anotherGeometry:Geometry)、

Intersects(anotherGeometry:Geometry)、Touches(anotherGeometry:Geometry)等都是符合其标准的。

对于OGR库后面会专门进行说明，这里简单说明一下就好。

GDAL源码剖析（十）之编写自己的扩展格式

一、简介

该节内容参考GDAL的英文原文：http://www.gdal.org/gdal_drivertut.html。

通常，可以通过从GDALDataset和GDALRasterBand继承来实现GDAL对新的数据格式支持。同时，还需要为这种格式创建一个GDALDriver的实例，让后通过GDALDriverManager将该新的驱动注册给GDAL系统。

该教程将为JDEM数据格式实现一个简单的只读驱动开始，进而使用RawRasterBand类，实现一个可创建和修改图像数据的格式以及更高级的问题。

强烈建议在看这篇文章之前，仔细理解GDAL数据模型的相关内容。地址为http://www.gdal.org /gdal_datamodel.html，或者参考上篇博文，GDAL体系架构。

该文章主要包含下面9部分内容：

从Dataset继承
从RasterBand继承
栅格驱动（Driver）
添加驱动到GDAL中
添加地理参考信息
金字塔（快视图）
创建文件
RawDataset和RawRasterBand类
元数据和其他外部扩展

二、从Dataset继承

从这里开始，我们将演示如何创建一个日本DEM数据的只读驱动。代码地址为：http:// www.gdal.org/jdemdataset.cpp.html。首先，从GDALDataset继承一个子类JDEMDataset。代码如下：

[cpp] view plaincopyprint?

class JDEMRasterBand;
class JDEMDataset : publicGDALPamDataset
{
friend class JDEMRasterBand;
VSILFILE *fp;
GByte abyHeader[1012];
public:
JDEMDataset();
~JDEMDataset();
static GDALDataset*Open( GDALOpenInfo* );
static int Identify( GDALOpenInfo* );
CPLErr GetGeoTransform( double* padfTransform );
const char *GetProjectionRef();
};

我们可以通过重载基类GDALDataset中的一些虚函数来为驱动重新实现这些特殊的功能。但是，Open()相对比较特殊，它不是基类 GDALDataset的虚函数。我们需要一个独立的函数来实现这个功能，因此我们把他声明为静态的（static）。将Open()声明为JDEMDataset的静态函数可以方便的访问JDEMDataset的私有方法去修改内容。

Open()函数具体实现如下：

[cpp] view plaincopyprint?

GDALDataset *JDEMDataset::Open(GDALOpenInfo * poOpenInfo)
{
if (!Identify(poOpenInfo))
return NULL;
/*——————————————————————– */
/* Confirm the requested access is supported. */
/* ——————————————————————–*/
if( poOpenInfo->eAccess == GA_Update)
{
CPLError( CE_Failure,CPLE_NotSupported,
“The JDEM driver does not support update access toexisting”
” datasets.\n” );
return NULL;
}
/*——————————————————————– */
/* Create a corresponding GDALDataset. */
/* ——————————————————————–*/
JDEMDataset *poDS;
poDS = new JDEMDataset();
poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, “rb”);
if (poDS->fp == NULL)
{
delete poDS;
return NULL;
}
/*——————————————————————– */
/* Read the header. */
/*——————————————————————– */
VSIFReadL( poDS->abyHeader, 1, 1012, poDS->fp );
poDS->nRasterXSize= JDEMGetField( (char*) poDS->abyHeader+ 23, 3 );
poDS->nRasterYSize= JDEMGetField( (char*) poDS->abyHeader+ 26, 3 );
if (poDS->nRasterXSize<= 0 || poDS->nRasterYSize<= 0 )
{
CPLError( CE_Failure,CPLE_AppDefined,
“Invalid dimensions : %d x %d”,
poDS->nRasterXSize,poDS->nRasterYSize);
delete poDS;
return NULL;
}
/* ——————————————————————–*/
/* Create band information objects. */
/*——————————————————————– */
poDS->SetBand(1, new JDEMRasterBand(poDS, 1 ));
/*——————————————————————– */
/* Initialize any PAM information. */
/*——————————————————————– */
poDS->SetDescription(poOpenInfo->pszFilename);
poDS->TryLoadXML();
/*——————————————————————– */
/* Check for overviews. */
/* ——————————————————————–*/
poDS->oOvManager.Initialize( poDS,poOpenInfo->pszFilename);
return( poDS );
}

第一步，打开文件，判断驱动是否支持该类型。这里我们应该知道，GDAL在打开文件的时候会调用每个注册驱动的Open函数将文件打开，直到有一个成功为止，可以参考GDAL源代码中的GDALOpen函数实现部分。如果指定的文件不是驱动所支持的类型，则它们必须返回NULL。如果格式是它们所支持的类型，但是数据已经被破坏，那么它们应该产生一个错误。

在文件打开之后文件的信息会保存在GDALOpenInfo对象中。GDALOpenInfo对象的公有成员如下：

[plain] view plaincopyprint?

char *pszFilename;
char **papszSiblingFiles;
GDALAccess eAccess;
int bStatOK;
int bIsDirectory;
FILE *fp;
int nHeaderBytes;
GByte *pabyHeader;

驱动可以通过检查这些值来判断是否支持这个格式的文件。如果pszFilename引用的是系统中的一个文件，那么bStatOK将被设置为TRUE，同时，如果文件打开成功，程序会读取前一千个字节的数据，并且将数据存储在pabyHeader中，nHeaderBytes保存实际读取的字节数。

在这个例子中，假设文件已经被打开并且已经读取到文件头部的一些信息。在这里，JDEM并没有魔术数字（magic numbers，这个没搞懂，个人理解是一个文件标识数字），因此我们只能检查不同的数据域是否合理。如果检查文件不是当前驱动所支持的格式，那么将返回NULL。检查函数的代码如下：

[cpp] view plaincopyprint?

int JDEMDataset::Identify(GDALOpenInfo * poOpenInfo)
{
/*——————————————————————– */
/* Confirm that the header has what appears to be dates in the */
/* expected locations. Sadly this isa relatively weak test. */
/*——————————————————————– */
if( poOpenInfo->nHeaderBytes < 50 )
return FALSE;
/* check if century values seem reasonable */
if( (!EQUALN((char *)poOpenInfo->pabyHeader+11,”19″,2)
&& !EQUALN((char *)poOpenInfo->pabyHeader+11,”20″,2))
|| (!EQUALN((char *)poOpenInfo->pabyHeader+15,”19″,2)
&&!EQUALN((char*)poOpenInfo->pabyHeader+15,”20″,2))
|| (!EQUALN((char *)poOpenInfo->pabyHeader+19,”19″,2)
&& !EQUALN((char *)poOpenInfo->pabyHeader+19,”20″,2)))
{
return FALSE;
}
return TRUE;
}

在真正的监测中，检查代码越全面越好。这里的监测相对比较弱。如果一个其他文件的在相应的位置具有相同的内容，那么它们很可能被错误地当作JDEM格式处理，最终导致不可预期的结果。

对于满足要求的文件格式，还有另外一个必须做的检查，就是检查文件是否可写，对于JDEM格式，需要进行判断，代码如下：

[cpp] view plaincopyprint?

if( poOpenInfo->eAccess== GA_Update )
{
CPLError( CE_Failure,CPLE_NotSupported,
“The JDEMdriver does not support update access to existing”
“datasets.\n” );
return NULL;
}

如果文件是我们支持的类型，满足上面的监测要求，我们需要创建一个数据集的实例，并在里面记录一些必要的信息。

[cpp] view plaincopyprint?

JDEMDataset *poDS;
poDS = new JDEMDataset();
poDS->fp = VSIFOpenL(poOpenInfo->pszFilename,”rb” );
if (poDS->fp == NULL)
{
delete poDS;
return NULL;
}

通常在这时我们已经将打开文件了，并且保存文件的指针。在接下来，尽可能的使用GDAL提供的接口中VSI*L的函数来访问文件，这些POSIX风格的API可以访问大文件，内存文件以及Zip压缩的文件。

接下来是要从读取到的文件头信息中提取图像的宽和高，nRasterXSize和nRasterYSize是在基类GDALDataset中定义的，这两个值必须在Open函数中进行赋值。代码如下：

[cpp] view plaincopyprint?

VSIFReadL( poDS->abyHeader,1, 1012, poDS->fp);
poDS->nRasterXSize = JDEMGetField((char *) poDS->abyHeader + 23, 3 );
poDS->nRasterYSize = JDEMGetField((char *) poDS->abyHeader + 26, 3 );
if (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0 )
{
CPLError( CE_Failure,CPLE_AppDefined,
“Invalid dimensions : %d x %d”,
poDS->nRasterXSize,poDS->nRasterYSize);
delete poDS;
return NULL;
}

通过SetBand()函数将所有的波段与当前的GDALDataset对象绑定。在下一节，我们将讨论JDEMRasterBand类的具体细节。

[cpp] view plaincopyprint?

poDS->SetBand( 1, new JDEMRasterBand( poDS,1 ));

最后，我们需要给dataset指定一个名称，并且调用GDALPamDataset类的函数TryLoadXML()重.aux.xml文件中读取一些辅助信息，更详细的内容请参阅GDALPamDataset类的相关函数。

[cpp] view plaincopyprint?

poDS->SetDescription( poOpenInfo->pszFilename );
poDS->TryLoadXML();
poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename );
return( poDS );

三、从RasterBand继承

和从GDALDataset类派生子类JDEMDataset一样，我们需要从GDALRasterBand类派生一个JDEMRasterBand类，JDEMRasterBand类的定义如下：

[cpp] view plaincopyprint?

class JDEMRasterBand : publicGDALPamRasterBand
{
friend class JDEMDataset;
int nRecordSize;
char* pszRecord;
public:
JDEMRasterBand(JDEMDataset *, int);
~JDEMRasterBand();
virtual CPLErr IReadBlock( int, int, void * );
};

构造函数可以任意定义，这里不写使用默认的就行，但是只能从Open()函数中调用。其他的虚函数必须和gdal_priv.h中定义的一致，例如IReadBlock()。构造函数实现代码如下：

[cpp] view plaincopyprint?

JDEMRasterBand::JDEMRasterBand( JDEMDataset*poDS, int nBand )
{
this->poDS = poDS;
this->nBand = nBand;
eDataType = GDT_Float32;
nBlockXSize = poDS->GetRasterXSize();
nBlockYSize = 1;
/* Cannot overflow as nBlockXSize <= 999 */
nRecordSize = nBlockXSize*5+ 9 + 2;
pszRecord = NULL;
}

下面列举的成员变量在GDALReasterBand中定义的，但是需要在子类的构造函数中进行初始化：

poDS: 指向基类GDALDataset的指针。
nBand: 对应数据集中的波段序号。
eDataType: 当前波段的数据类型。
nBlockXSize: 当前波段的块宽度。
nBlockYSize: 当前波段的块高度。

GDALDataType类型的定义在文件gdal.h中，包括GDT_Byte、GDT_UInt16、GDT_Int16和 GDT_Float32。块的大小记录数据的实际或有效的大小。对于分块数据集这个值将是一个分块的大小，而对于其他大多数数据而言这个值将是一行。

接下来，我们将实现真正读取影像数据的代码，IReadBlock()函数。

[cpp] view plaincopyprint?

CPLErr JDEMRasterBand::IReadBlock(int nBlockXOff,int nBlockYOff,
void * pImage )
{
JDEMDataset *poGDS= (JDEMDataset *) poDS;
int i;
if (pszRecord == NULL)
{
if (nRecordSize< 0)
return CE_Failure;
pszRecord = (char*) VSIMalloc(nRecordSize);
if (pszRecord == NULL)
{
CPLError(CE_Failure,CPLE_OutOfMemory,
“Cannot allocate scanline buffer”);
nRecordSize = -1;
return CE_Failure;
}
}
VSIFSeekL( poGDS->fp, 1011 + nRecordSize*nBlockYOff, SEEK_SET);
VSIFReadL( pszRecord,1, nRecordSize, poGDS->fp );
if( !EQUALN((char *) poGDS->abyHeader,pszRecord,6))
{
CPLError( CE_Failure,CPLE_AppDefined,
“JDEM Scanline corrupt. Perhaps file was not transferred\n”
“in binary mode?” );
return CE_Failure;
}
if( JDEMGetField( pszRecord + 6, 3 ) != nBlockYOff+ 1 )
{
CPLError( CE_Failure,CPLE_AppDefined,
“JDEM scanline out of order, JDEM driver doesnot\n”
“currently support partial datasets.” );
return CE_Failure;
}
for( i = 0; i < nBlockXSize;i++ )
((float *) pImage)[i] = (float)
(JDEMGetField( pszRecord+ 9 + 5 * i, 5) * 0.1);
return CE_None;
}

在这里需要注意的有：

需要将GDALRasterBand::poDS成员变量传递给子类，如果你的RasterBand类需要访问数据集的私有成员，那么需要吧RasterBand类声明为Dataset的友元类。
如果在处理中出现错误，请使用CPLError()提示，并且返回CE_Failure，否则返回 CE_None即可。
pImage缓存应该用一个块的数据来进行赋值，块大的小nBlockXSize 和nBlockYSize在波段中定义。pImage的数据类型应该和栅格波段对象中声明的数据类型eDataType一致。
nBlockXOff和nBlockYOff是块的偏移量，因此对于一个128×128大小的块，序号为1和1的块应该对应原始数据是从(128,128)到(255,255)的数据。

四、栅格驱动（Driver）

虽然JDEMDataset和JDEMRasterBand已经可以读数据，但是GDAL库仍然不清楚JDEMDataset的任何信息。这可以通过GDALDriverManager来完成。为了将我们自己实现的驱动注册到GDAL库中，我们需要重新实现一个注册函数，函数定义写在gdal_frmts.h文件中：

[cpp] view plaincopyprint?

CPL_C_START
void CPL_DLL GDALRegister_JDEM(void);
CPL_C_END

函数的实现部分如下：

[cpp] view plaincopyprint?

void GDALRegister_JDEM()
{
GDALDriver *poDriver;
if( GDALGetDriverByName(“JDEM” ) == NULL )
{
poDriver = new GDALDriver();
poDriver->SetDescription(“JDEM” );
poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
“Japanese DEM (.mem)” );
poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
“frmt_various.html#JDEM” );
poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, “mem”);
poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, “YES” );
poDriver->pfnOpen= JDEMDataset::Open;
poDriver->pfnIdentify= JDEMDataset::Identify;
GetGDALDriverManager()->RegisterDriver( poDriver );
}
}

可以使用GDAL_CHECK_VERSION宏来判断GDAL的版本（从1.5.0版本之后才有），这个宏可以用来使GDAL扩展一些额外的库，比如编写一个GDAL的插件驱动。对于这种情况，尽可能使用这个宏来判断GDAL的版本信息，对于一些老的版本，GDAL可能会不支持。

当第一次被调用的时候，注册函数将创建一个GDALDriver的对象，并且通过GDALDriverManager来注册。在用GDALDriverManager注册驱动之前，需要设置以下的成员变量：

驱动对应格式的名称，不能和其他驱动名字冲突。通常在3-5个字符长度，并且匹配驱动类名字的前缀。（手工设置）

GDAL_DMD_LONGNAME: 文件格式的详细描述，长度一般在50-60个字符。（手工设置）
GDAL_DMD_HELPTOPIC: 如果有的话，是关于这个驱动帮助文档的名称。在这个例子中，JDEM格式包含在frmt_various.html#JDEM位置描述。（可选）
GDAL_DMD_EXTENSION: 该类型文件的扩展名。如果扩展名多于一个，则最好选择最通用的那个，或者直接为空。（可选）
GDAL_DMD_MIMETYPE: 该格式数据的标准用途类型，例如“image/png”。（可选）
GDAL_DMD_CREATIONOPTIONLIST: 用于描述创建时的选项。可以参考geotiff驱动的实现代码。（可选）
GDAL_DMD_CREATIONDATATYPES: 支持创建数据集的全部类型列表。如果存在Create()方法这些类型都将被支持。如果存在CreateCopy()方法，该类型列表将是那些支持无损输出的类型。比如，一个格式使用了CreateCopy()方法，如果可以写为Float32类型，那么Byte、Int16和UInt16都应该被支持（因为它们都可以用Float32表示）。
pfnOpen: 打开这种格式文件的函数指针。（可选）
pfnCreate: 创建这中格式的函数指针。（可选）
pfnCreateCopy: 创建一个从其它数据源复制而来的新数据集，但是不需要更新的函数。（可选）
pfnDelete: 删除这种格式数据集函数指针。（可选）
pfnUnloadDriver: 这个函数只有在驱动被销毁的时候才被调用。在驱动层被用来清除数据。很少用到。（可选）

五、添加驱动到GDAL库

GDALRegister_JDEM()函数必须在上层的程序中调用，才能使支持JDEM格式的数据。通常需要以下的工作来完成：

在gdal/frmts 创建一个目录，目录的名字与驱动的短名称相同。
在目录中添加文件GNUmakefile和makefile.vc，具体格式可以参考其他目录中，如jdem目录中的文件。
添加dataset和RasterBand的实现部分（这里值得是源代码）。通常命名规则值<短名称> dataset.cpp，如jdemdataset.cpp。
在gdal/gcore/gdal_frmts.h文件中添加一个注册函数声明，如GDALRegister_JDEM())。
在frmts/gdalallregister.cpp文件中的注册函数中调用该注册函数，可以使用ifdef来进行判断，可以参考相关代码。
在文件GDALmake.opt.in (和GDALmake.opt)中的GDAL_FORMATS宏中添加格式的短名称。
在文件frmts/makefile.vc中的EXTRAFLAGS宏添加该项，具体参考相关代码。

完成上面的步骤，重新编译生成GDAL库后，那么所有的工具都会支持这种新的文件格式，可以使用gdalinfo工具来进行测试打开数据和显示数据的信息，可以使用gdal_translate工具测试图像的读取。

六、添加地理参考信息

接下来我们继续完善之前的工作，这里将使驱动支持地理参考信息，我们需要重写JDEMDataset类中的两个虚函数，这两个函数的定义在GDALRasterDataset基类中。两个函数定义如下：

[cpp] view plaincopyprint?

virtual const char *GetProjectionRef(void);
virtual CPLErr GetGeoTransform(double * );

在JDEMDataset类中添加函数声明：

[cpp] view plaincopyprint?

CPLErr GetGeoTransform(double * padfTransform);
const char *GetProjectionRef();

重新实现的GetGeoTransform()函数只是复制地理转换矩阵到缓存中。GetGeoTransform()函数可能经常使用，因此一般最好使该函数体小。在很多时候，使用Open()函数已经获取到地理转换的信息，然后用这个函数进行复制里面的信息即可。需要注意的是，转换矩阵中的起始点坐标对应像素左上角的坐标，而不是左上角象元的中心坐标。

[cpp] view plaincopyprint?

CPLErr JDEMDataset::GetGeoTransform(double * padfTransform)
{
double dfLLLat, dfLLLong,dfURLat, dfURLong;
dfLLLat = JDEMGetAngle((char *) abyHeader+ 29 );
dfLLLong = JDEMGetAngle((char *) abyHeader+ 36 );
dfURLat = JDEMGetAngle((char *) abyHeader+ 43 );
dfURLong = JDEMGetAngle((char *) abyHeader+ 50 );
padfTransform[0] = dfLLLong;
padfTransform[3] = dfURLat;
padfTransform[1] = (dfURLong- dfLLLong) / GetRasterXSize();
padfTransform[2] = 0.0;
padfTransform[4] = 0.0;
padfTransform[5] = -1 * (dfURLat- dfLLLat) / GetRasterYSize();
return CE_None;
}

GetProjectionRef ()方法将返回一个字符串，使用的是OGC WKT格式的坐标系统的定义。在这个例子中，这种格式的所有数据只有一种坐标系统。但是在比较复杂的数据格式中，我们需要使用OGRSpatialReference来针对特定情形来进行导出坐标系统。

[cpp] view plaincopyprint?

const char *JDEMDataset::GetProjectionRef()
{
return( “GEOGCS[\”Tokyo\”,DATUM[\”Tokyo\”,SPHEROID[\”Bessel1841\”,6377397.155,299.1528128,AUTHORITY[\”EPSG\”,7004]],TOWGS84[-148,507,685,0,0,0,0],AUTHORITY[\”EPSG\”,6301]],PRIMEM[\”Greenwich\”,0,AUTHORITY[\”EPSG\”,8901]],UNIT[\”DMSH\”,0.0174532925199433,AUTHORITY[\”EPSG\”,9108]],AUTHORITY[\”EPSG\”,4301]]”);
}

到这里已经完成了JDEM驱动的代码，整个源代码可以在jdemdataset.cpp（地址为http://www.gdal.org/jdemdataset.cpp.html）文件中看到。

七、金字塔（快视图）

GDAL可以使用GDALRasterBand::GetOverview()和相关的函数来对文件建立金字塔。然后，这些内容已经超出了本文档的讨论范围，具体可以参考GeoTIFF驱动(gdal/frmts/gtiff/geotiff.cpp)中的相关代码。

通过GDALRasterBand的函数HasArbitraryOverviews()可以获取文件中是否含有金字塔信息，返回TRUE表示有。在这个例子中，栅格波段对象应该重载实现RasterIO()函数，通过该函数可以有效的访问图像或者重采样等。这个函数比较复杂，需要判断很多参数是否正确等，关于这个的例子可以在OGDI和ECW格式中找到相关的代码。

然而，很多图像的金字塔都使用GDAL默认的外部金字塔格式来保存，实际上是一个TIFF文件，但是后缀名为.ovr，为了使用读取和创建金字塔，需要对GDALDataset类中的oOvManager进行初始化。通常情况下，在Open()函数的最后调用类似下面的代码（PAM TryLoadXML()函数之前）。

[cpp] view plaincopyprint?

poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename );

这样依赖，该格式可以读取和创建金字塔。一般简单的图像格式，都可以使用这种方式，除非一些特殊的格式有自定义的金字塔格式。

八、创建文件

有两种方法来创建文件。第一种是调用CreateCopy()函数，可以从源影像中获取影像信息并写入到结果图像中；第二种是调用Create动态创建文件，包含Create函数和用来设置各种信息的函数。

第一种方法的优点是在创建文件的所有的信息都从源图像中获取并写入到结果图像中。尤其是对于文件那些特别重要的信息例如颜色表、参照系等信息。关于CreateCopy方法的缺点是对于一些信息没有相应的设置方法，例如min/max、scaling、description和GCPs。

第二种方法的优点在于可以创建一个空的新文件，并且在需要的时候把结果写入其中。对于一个影像来说，动态创建不需要提前获取影像的所有信息。

对于非常重要的格式而言，两种方法最好都支持。

1、CreateCopy

GDALDriver::CreateCopy()方法可以直接调用，因此只需要知道调用的参数既可。不过需要注意以下的细节：

如果参数bStrict为FALSE，那么驱动器将对数据自行做一些合适的处理。特别是对于那些并不完全等同格式。
可以实现CreateCopy进度提示。回调函数的返回值表示是否需要继续进行，并且过程最好是在合理的时间单位中被调用（这个例子没有示范）。
特殊的创建选项必须在帮助中说明。如果存在”NAME=VALUE”的格式，可以使用CPLFetchNameValue()函数来进行设置，比如JPEGCreateCopy()函数就有QUALITY和PROGRESSIVE两个标志来说明相关的创建选项。
返回的GDALDataset指针是只读或者更新模式。在通常情况下返回更新模式，否则返回只读模式即可。

JPEG格式对应的CreateCopy的代码如下：

[cpp] view plaincopyprint?

static GDALDataset *
JPEGCreateCopy( const char * pszFilename, GDALDataset*poSrcDS,
int bStrict, char ** papszOptions,
GDALProgressFuncpfnProgress, void* pProgressData )
{
int nBands = poSrcDS->GetRasterCount();
int nXSize = poSrcDS->GetRasterXSize();
int nYSize = poSrcDS->GetRasterYSize();
int nQuality = 75;
int bProgressive = FALSE;
// ————————————————–
// Some somerudimentary checks
// ————————————————–
if( nBands != 1&& nBands != 3 )
{
CPLError( CE_Failure,CPLE_NotSupported,
“JPEG driver doesn’t support %d bands. Must be 1 (grey) ”
“or 3 (RGB) bands.\n”, nBands );
return NULL;
}
if( poSrcDS->GetRasterBand(1)->GetRasterDataType()!= GDT_Byte && bStrict )
{
CPLError( CE_Failure,CPLE_NotSupported,
“JPEG driver doesn’t support data type %s. ”
“Only eight bit byte bands supported.\n”,
GDALGetDataTypeName(
poSrcDS->GetRasterBand(1)->GetRasterDataType()) );
return NULL;
}
// ————————————————–
// What optionshas the user selected?
// ————————————————–
if( CSLFetchNameValue(papszOptions,”QUALITY”)!= NULL )
{
nQuality = atoi(CSLFetchNameValue(papszOptions,”QUALITY”));
if( nQuality <10 || nQuality > 100 )
{
CPLError( CE_Failure,CPLE_IllegalArg,
“QUALITY=%s is not a legal value in the range10-100.”,
CSLFetchNameValue(papszOptions,”QUALITY”) );
return NULL;
}
}
if( CSLFetchNameValue(papszOptions,”PROGRESSIVE”)!= NULL )
{
bProgressive = TRUE;
}
// ————————————————–
// Create thedataset.
// ————————————————–
FILE *fpImage;
fpImage = VSIFOpen(pszFilename, “wb”);
if( fpImage == NULL )
{
CPLError( CE_Failure,CPLE_OpenFailed,
“Unable to create jpeg file %s.\n”,
pszFilename );
return NULL;
}
// ————————————————–
// InitializeJPG access to the file.
// ————————————————–
struct jpeg_compress_structsCInfo;
struct jpeg_error_mgrsJErr;
sCInfo.err = jpeg_std_error( &sJErr);
jpeg_create_compress( &sCInfo );
jpeg_stdio_dest( &sCInfo,fpImage );
sCInfo.image_width= nXSize;
sCInfo.image_height= nYSize;
sCInfo.input_components= nBands;
if( nBands == 1 )
{
sCInfo.in_color_space= JCS_GRAYSCALE;
}
else
{
sCInfo.in_color_space= JCS_RGB;
}
jpeg_set_defaults( &sCInfo);
jpeg_set_quality( &sCInfo,nQuality, TRUE);
if( bProgressive )
jpeg_simple_progression( &sCInfo );
jpeg_start_compress( &sCInfo,TRUE );
// ————————————————–
// Loop overimage, copying image data.
// ————————————————–
GByte *pabyScanline;
CPLErr eErr;
pabyScanline = (GByte*) CPLMalloc( nBands* nXSize );
for( int iLine = 0; iLine< nYSize; iLine++)
{
JSAMPLE *ppSamples;
for( int iBand = 0; iBand< nBands; iBand++)
{
GDALRasterBand * poBand= poSrcDS->GetRasterBand(iBand+1 );
eErr = poBand->RasterIO( GF_Read,0, iLine, nXSize,1,
pabyScanline + iBand,nXSize, 1, GDT_Byte,
nBands, nBands* nXSize );
}
ppSamples = pabyScanline;
jpeg_write_scanlines( &sCInfo, &ppSamples, 1 );
}
CPLFree( pabyScanline);
jpeg_finish_compress( &sCInfo );
jpeg_destroy_compress( &sCInfo );
VSIFClose( fpImage);
return (GDALDataset*) GDALOpen( pszFilename,GA_ReadOnly );
}

2、动态创建

在动态创建的例子中，没有源数据集。但是提供了文件的大小、波段数和像素的数据类型。对于其他的一些信息（例如地理参照系等），将在以后通过特定的函数进行设置。

下面使用的是一个简单实现了PCI.aux文件的创建示例。使用自己写的方法创建了一个空白的影像，并且在最后调用GDALOpen，使用更新方式打开数据，返回打开的数据集指针。

[cpp] view plaincopyprint?

GDALDataset *PAuxDataset::Create(const char * pszFilename,
int nXSize, int nYSize, int nBands,
GDALDataTypeeType,
char ** // papszParmList )
{
char *pszAuxFilename;
//————————————————–
// Verify inputoptions.
//————————————————–
if( eType != GDT_Byte && eType!= GDT_Float32 && eType != GDT_UInt16
&& eType != GDT_Int16)
{
CPLError( CE_Failure,CPLE_AppDefined,
“Attempt to create PCI .Aux labelled dataset with anillegal\n”
“data type (%s).\n”,
GDALGetDataTypeName(eType));
return NULL;
}
//————————————————–
// Try to createthe file.
//————————————————–
FILE *fp;
fp = VSIFOpen( pszFilename, “w”);
if( fp == NULL )
{
CPLError( CE_Failure,CPLE_OpenFailed,
“Attempt to create file `%s’ failed.\n”,
pszFilename );
return NULL;
}
//————————————————–
// Just writeout a couple of bytes to establish the binary
// file, andthen close it.
//————————————————–
VSIFWrite( (void*) “\0\0”, 2, 1, fp );
VSIFClose( fp);
//————————————————–
// Create theaux filename.
//————————————————–
pszAuxFilename = (char*) CPLMalloc(strlen(pszFilename)+5);
strcpy( pszAuxFilename,pszFilename );;
for( int i = strlen(pszAuxFilename)-1; i> 0; i– )
{
if( pszAuxFilename[i] == ‘.’ )
{
pszAuxFilename[i]= ‘\0’;
break;
}
}
strcat( pszAuxFilename,”.aux” );
//————————————————–
// Open thefile.
//————————————————–
fp = VSIFOpen( pszAuxFilename, “wt”);
if( fp == NULL )
{
CPLError( CE_Failure,CPLE_OpenFailed,
“Attempt to create file `%s’ failed.\n”,
pszAuxFilename );
return NULL;
}
//————————————————–
// We need towrite out the original filename but without any
// pathcomponents in the AuxilaryTarget line. Do so now.
//————————————————–
int iStart;
iStart = strlen(pszFilename)-1;
while( iStart >0 && pszFilename[iStart-1] != ‘/’
&& pszFilename[iStart-1]!= ‘\\’ )
iStart–;
VSIFPrintf( fp,”AuxilaryTarget: %s\n”, pszFilename + iStart);
//————————————————–
// Write out theraw definition for the dataset as a whole.
//————————————————–
VSIFPrintf( fp,”RawDefinition: %d %d %d\n”,
nXSize, nYSize,nBands );
//————————————————–
// Write out adefinition for each band. We alwayswrite band
// sequentialfiles for now as these are pretty efficiently
// handled byGDAL.
//————————————————–
int nImgOffset = 0;
for( int iBand = 0; iBand< nBands; iBand++)
{
const char * pszTypeName;
int nPixelOffset;
int nLineOffset;
nPixelOffset = GDALGetDataTypeSize(eType)/8;
nLineOffset = nXSize* nPixelOffset;
if( eType == GDT_Float32 )
pszTypeName = “32R”;
else if( eType == GDT_Int16)
pszTypeName = “16S”;
else if( eType == GDT_UInt16)
pszTypeName = “16U”;
else
pszTypeName = “8U”;
VSIFPrintf( fp,”ChanDefinition-%d: %s %d %d %d %s\n”,
iBand+1, pszTypeName,
nImgOffset, nPixelOffset,nLineOffset,
#ifdef CPL_LSB
“Swapped”
#else
“Unswapped”
#endif
);
nImgOffset += nYSize* nLineOffset;
}
//————————————————–
// Cleanup
//————————————————–
VSIFClose( fp);
return (GDALDataset*) GDALOpen( pszFilename,GA_Update );
}

如果文件格式支持动态创建或支持更新都需要实现GDALRasterBand中的IWriteBlock()方法，它类似于IReadBlock()。并且，在GDALRasterBand的析构函数中调用FlushCache()方法是比较危险的。因此，在析构方法调用之前必须确保任何保存在缓存块的数据都已经清除。

九、RawDataset和RawRasterBand类

很多文件格式将图像数据存储在一个使用二进制文件中，使用规则的按照行存储的格式。为了读取这些类似的格式，这里使用了一个RawDataset类和RawRasterBand类来进行处理，关于这两个类的定义可以在gdal/frmts/raw中找到。

在这个例子中没有用到这个类，如果有用的，可以从RawRasterBand派生，数据集可以从RawDataset派生即可。

Open()函数调用后通过构造函数实例化栅格波段里面的所有信息。可以查看PNM的驱动就是通过这样的方式来创建他的栅格波段的。

RawRasterBand可以获取下面的参数：

poDS: 该波段的GDALDataset指针，这个数据集必须是从RawRasterDataset派生。
nBand: 波段序号，从1开始。
fpRaw: 图像数据打开后的文件句柄，类型是FILE *。
nImgOffset: 第一个象元的偏移量byte。
nPixelOffset: 第一个象元从在一行中的偏移量byte。
nLineOffset: 起始扫描线和下一个扫描线的偏移量byte。
eDataType: GDALDataType类型的数据类型，用来表示图像在磁盘上的存储类型。
bNativeOrder: （这句没看懂啥意思，原文放这里吧）FALSE if the data is not in the same endianness asthe machine GDAL is running on.个人感觉，这个的意思大概是如果GDAL读取的数据的顺序和存储的顺序不同的话，使用TRUE，然后GDAL会对这些数据自动进行交换。

一个简单的裸数据读取的例子在gdal/frmts/raw目录中，里面有大量的关于裸数据的说明。

十、元数据和其他外部扩展

在GDAL数据模型中还有很多其他项目，在GDALDataset和GDALRasterBand中都存在对应的虚函数。它们包括：

Metadata: 关于数据集或者波段的Name/value。GDALMajorObject（包括其子类）都能支持元数据，可以在Open()函数中通过调用SetMetadataItem()函数来进行设置。SAR_CEOS(frmts/ceos2/sar_ceosdataset.cpp)和GeoTIFF驱动中都有关于读取元数据的例子。
ColorTables: GDT_Byte类型的波段可以包含与它们相关联的颜色表信息。在文件frmts/png/pngdataset.cpp中就是一个关于支持颜色表的例子。
ColorInterpretation: PNG驱动中包括一个返回波段被表示为红，绿，蓝，透明度或者灰度的描述信息。
GCPs: GDALDataset含有一系列的地面控制点与他们相对应的空间参考坐标进行关联（通过GetGeoTransform来映射转换）。MFF2格式（gdal/frmts.raw.hkvdataset.cpp）就是一个简单的支持GCP的例子。
NoDataValue: GetNoDataValue()可以获取波段的NODATA值。具体可以参考frmts/raw/pauxdataset.cpp驱动的相关代码。
Category Names: GetCategoryNames()函数可以根据影象的名字将其分类。不过目前还没有哪个驱动用到这个信息。

GDAL源码剖析（十一）之OGR投影说明

一、简介

本文参考英文地址：http://www.gdal.org/ogr/osr_tutorial.html。

OGRSpatialReference类和OGRCoordinateTransformation类主要用来提供定义坐标系统（投影和水准面）和转换坐标。这两个类都基于OpenGIS的坐标转换说明，并且使用Well Known Text格式来进行表述坐标系统。

一些关于OpenGIS坐标系统的资料，以及空间参考坐标抽象模型文件可以在OGC（Open Geospatial Consortium）的网站上找到。GeoTIFF投影转换列表（GeoTIFF Projections Transform List）网页可以更好的帮助你理解WKT的规则，同时EPSG的网站也是很有用的资料。

二、定义地理坐标系统

坐标系统使用OGRSpatialReference类来进行封装。这里提供了数种初始化OGRSpatialReference类的方式。这里有两类主要的坐标系统，第一种是地理坐标系统（位置信息使用经纬度来表示的），第二种是投影坐标系统（比如UTM-通用横轴墨卡托投影，位置信息使用米或者英尺来表示）。

一个地理坐标系统包含的信息有一个大地基准（里面含有一个使用长半轴和扁率的倒数来表示的托球体），一个中央经线（通常是本初子午线，也就是0度经线Greenwich）, 此外还有一个角度的度量单位，使用度而不是弧度。如果含有这些信息，就可以构造一个有效的地理坐标系统。

[cpp] view plaincopyprint?

OGRSpatialReference oSRS;
oSRS.SetGeogCS( “Mygeographic coordinate system”,
“WGS_1984”,
“My WGS84 Spheroid”,
SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
“Greenwich”, 0.0,
“degree”, SRS_UA_DEGREE_CONV );

在上面的代码中，名称为“My geographic coordinate system”，“My WGS84 Spheroid”，“Greenwich”和“degree”的并不是关键词，这些主要是用来给用户进行说明的。然而，“WGS_1984”是一个定义大地基准的关键词，注意：这里的大地基准必须是一个有效的大地基准！（这句话的意思，前面的那些字符串就是随便指定的，用来显示的，后面的WGS_1984这个位置的字符串，必须是一个有效的，不能随便命名，具体后面会说到）。

OGRSpatialReference可以使用一些常用的字符串来进行建立一个常用的坐标系统，这些字符串包括“NAD27”、“NAD83”，“WGS72”和“WGS84”等，使用的函数是SetWellKnownGeogCS()，使用方式见下面。

[cpp] view plaincopyprint?

oSRS.SetWellKnownGeogCS( “WGS84” );

而且，还可以使用EPSG数据库定义的GCS代码来定义一个地理坐标系统，如：

[cpp] view plaincopyprint?

oSRS.SetWellKnownGeogCS( “EPSG:4326” );

为了方便的和其他库进行关联，这里可以转换为OpenGIS的WKT格式。同样OGRSpatialReference可以使用一个WKT来进行初始化，或者将里面的信息导出为WKT格式。

[cpp] view plaincopyprint?

char *pszWKT =NULL;
SRS.SetWellKnownGeogCS( “WGS84” );
oSRS.exportToWkt( &pszWKT );
printf( “%s\n”, pszWKT );

上面的语句会输出下面的内容：

[plain] view plaincopyprint?

GEOGCS[“WGS84”,DATUM[“WGS_1984”,SPHEROID[“WGS84”,6378137,298.257223563,AUTHORITY[“EPSG”,7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[“EPSG”,6326]],PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,8901]],UNIT[“DMSH”,0.0174532925199433,AUTHORITY[“EPSG”,9108]],AXIS[“Lat”,NORTH],AXIS[“Long”,EAST],AUTHORITY[“EPSG”,4326]]

将上面的字符串格式调整成更好看的样子：

[plain] view plaincopyprint?

GEOGCS[“WGS 84”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.257223563,
AUTHORITY[“EPSG”,7030]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY[“EPSG”,6326]],
PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,8901]],
UNIT[“DMSH”,0.0174532925199433,AUTHORITY[“EPSG”,9108]],
AXIS[“Lat”,NORTH],
AXIS[“Long”,EAST],
AUTHORITY[“EPSG”,4326]]

函数OGRSpatialReference::importFromWkt()可以从一个WKT定义的坐标系统来构造一个OGRSpatialReference类对象。

三、定义投影坐标系统

一个投影坐标系统（比如UTM，兰伯特等角圆锥投影等）需要建立在一个地理坐标系统之上，在投影坐标系统中，坐标点使用米或者英尺等长度单位来表示，同时也可以用经纬度的角度坐标来表示。下面将定义一个UTM的第17带的投影坐标系统，基于WGS84的大地基准椭球体。

[cpp] view plaincopyprint?

OGRSpatialReference oSRS;
oSRS.SetProjCS( “UTM 17(WGS84) in northern hemisphere.” );
oSRS.SetWellKnownGeogCS(“WGS84” );
oSRS.SetUTM( 17, TRUE );

首先调用SetProjCS()函数设置投影坐标系统的名称，然后使用函数SetWellKnownGeogCS()指定地理坐标系统，最后调用函数SetUTM()设置投影转换参数信息。完成这些工作之后就定义了一个有效的投影坐标系统。这里必须要注意定义OGRSpatialReference的顺序！

上面定义的投影使用WKT表示的形式如下，注意UTM17会使用横轴墨卡托的分带投影参数来表示。

[plain] view plaincopyprint?

PROJCS[“UTM 17 (WGS84) in northernhemisphere.”,
GEOGCS[“WGS 84”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.257223563,
AUTHORITY[“EPSG”,7030]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY[“EPSG”,6326]],
PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,8901]],
UNIT[“DMSH”,0.0174532925199433,AUTHORITY[“EPSG”,9108]],
AXIS[“Lat”,NORTH],
AXIS[“Long”,EAST],
AUTHORITY[“EPSG”,4326]],
PROJECTION[“Transverse_Mercator”],
PARAMETER[“latitude_of_origin”,0],
PARAMETER[“central_meridian”,-81],
PARAMETER[“scale_factor”,0.9996],
PARAMETER[“false_easting”,500000],
PARAMETER[“false_northing”,0]]

这里提供了很多设置投影坐标的方法，包括SetTM()(横轴墨卡托投影), SetLCC()(兰伯特等角圆锥投影)和SetMercator()。

四、获取坐标系统信息

一旦一个OGRSpatialReference对象进行创建，那么就可以获取里面的各种信息，比如可以使用OGRSpatialReference::IsProjected()和OGRSpatialReference::IsGeographic()方法来判断坐标系统是地理坐标系统还是投影坐标系统。使用函数OGRSpatialReference::GetSemiMajor()、OGRSpatialReference::GetSemiMinor()和OGRSpatialReference::GetInvFlattening()方法可以用来获取椭球体信息，分别是椭球体的长半轴，短半轴以及扁率的倒数。使用OGRSpatialReference::GetAttrValue()方法可以用来获取PROJCS、GEOGCS、DATUM、SPHEROID和PROJECTION的名称字符串。使用OGRSpatialReference::GetProjParm()方法可以获取投影参数信息。使用OGRSpatialReference::GetLinearUnits()方法可以获取长度单位类型，并将其转换为米。

下面的代码是一个获取坐标系统信息的示例，摘自ogr_srs_proj4.cpp文件中。使用GetAttrValue()获取投影名称，使用GetProjParm()获取投影参数信息。GetAttrValue()函数的第一个值节点就是WKT字符串的描述信息。投影参数的宏定义（比如SRS_PP_CENTRAL_MERIDIAN）和函数GetProjParm()一起使用，可以用来获取投影的参数。更多的使用方式可以参考文件ogrspatialreference.cpp中的相关代码。

[cpp] view plaincopyprint?

const char *pszProjection = poSRS->GetAttrValue(“PROJECTION”);
if( pszProjection == NULL )
{
if( poSRS->IsGeographic() )
sprintf(szProj4+strlen(szProj4), “+proj=longlat ” );
else
sprintf(szProj4+strlen(szProj4), “unknown ” );
}
else if(EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
sprintf(szProj4+strlen(szProj4),
“+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f “,
poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0));
}
…
const char *pszProjection = poSRS->GetAttrValue(“PROJECTION”);
if( pszProjection == NULL )
{
if( poSRS->IsGeographic() )
sprintf(szProj4+strlen(szProj4), “+proj=longlat ” );
else
sprintf(szProj4+strlen(szProj4), “unknown ” );
}
else if(EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
sprintf(szProj4+strlen(szProj4),
“+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f “,
poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0));
}
…

五、坐标转换

OGRCoordinateTransformation类可以用来在不同的坐标系统中进行坐标转换。可以使用函数OGRCreateCoordinateTransformation()创建一个新的坐标转换对象，然后使用OGRCoordinateTransformation::Transform()方法来进行坐标转换。

[cpp] view plaincopyprint?

OGRSpatialReference oSourceSRS,oTargetSRS;
OGRCoordinateTransformation *poCT;
double x, y;
oSourceSRS.importFromEPSG(atoi(papszArgv[i+1]) );
oTargetSRS.importFromEPSG(atoi(papszArgv[i+2]) );
poCT = OGRCreateCoordinateTransformation(&oSourceSRS,
&oTargetSRS );
x = atof( papszArgv[i+3] );
y = atof( papszArgv[i+4] );
if( poCT == NULL || !poCT->Transform( 1, &x,&y ) )
printf(“Transformation failed.\n” );
else
printf(“(%f,%f) -> (%f,%f)\n”,
atof(papszArgv[i+3] ),
atof(papszArgv[i+4] ),
x, y );

一对点转换失败的原因有，首先，OGRCreateCoordinateTransformation()函数执行失败，通常的原因是不知道指定的两个投影之间的转换关系。这个可能是试用了PROJ.4库不支持的投影，不同的椭球体之间的转换关系没有定义，或者是其中的一个坐标系统没有定义完全。如果函数OGRCreateCoordinateTransformation() 执行失败，那么其返回值将是NULL。

OGRCoordinateTransformation::Transform() 方法本身页可能执行失败。这个可能的原因也有上面的问题，或者是转换的点里面有一个以上没有定义的数字。函数Transform()执行成功是返回TRUE，如果有一个点转换失败都会返回FALSE。

坐标转换也可以支持三维点的坐标转换，会自动根据不同的椭球地的基准面调整高程值。这个可以用来在不同的垂直基准面进行坐标转换。如果没有Z值，系统会认为所有的点都是在水准面上。

下面的例子演示了如果从一个投影坐标系统中的点转换为该投影中的地理坐标系统中的点，将米表示的坐标转换为经纬度表示的坐标。

[cpp] view plaincopyprint?

OGRSpatialReference oUTM, *poLatLong;
OGRCoordinateTransformation *poTransform;
oUTM.SetProjCS(“UTM 17 /WGS84”);
oUTM.SetWellKnownGeogCS(“WGS84” );
oUTM.SetUTM( 17 );
poLatLong = oUTM.CloneGeogCS();
poTransform = OGRCreateCoordinateTransformation( &oUTM,poLatLong );
if( poTransform == NULL )
{
…
}
…
if( !poTransform->Transform( nPoints, x,y, z ) )
…

六、其他语言接口

OGR的空间参考提供了一个C语言的接口，定义在ogr_srs_api.h文件中，Python的接口定义在osr.py文件中。所有的接口名称和C++的接口都很相似，但是C和Python中有些方法没有进行提供。

1、C语言接口

[cpp] view plaincopyprint?

typedef void *OGRSpatialReferenceH;
typedef void *OGRCoordinateTransformationH;
OGRSpatialReferenceH OSRNewSpatialReference( const char *);
void OSRDestroySpatialReference( OGRSpatialReferenceH );
int OSRReference( OGRSpatialReferenceH );
int OSRDereference( OGRSpatialReferenceH );
OGRErr OSRImportFromEPSG( OGRSpatialReferenceH, int );
OGRErr OSRImportFromWkt( OGRSpatialReferenceH, char ** );
OGRErr OSRExportToWkt( OGRSpatialReferenceH, char ** );
OGRErr OSRSetAttrValue( OGRSpatialReferenceH hSRS, const char * pszNodePath,
const char *pszNewNodeValue );
const char *OSRGetAttrValue( OGRSpatialReferenceH hSRS,
const char * pszName, intiChild);
OGRErr OSRSetLinearUnits( OGRSpatialReferenceH, const char *, double );
double OSRGetLinearUnits( OGRSpatialReferenceH, char ** );
int OSRIsGeographic( OGRSpatialReferenceH );
int OSRIsProjected( OGRSpatialReferenceH );
int OSRIsSameGeogCS( OGRSpatialReferenceH, OGRSpatialReferenceH );
int OSRIsSame( OGRSpatialReferenceH, OGRSpatialReferenceH );
OGRErr OSRSetProjCS( OGRSpatialReferenceH hSRS, const char * pszName );
OGRErr OSRSetWellKnownGeogCS( OGRSpatialReferenceH hSRS,
const char * pszName );
OGRErr OSRSetGeogCS( OGRSpatialReferenceH hSRS,
const char * pszGeogName,
const char *pszDatumName,
const char *pszEllipsoidName,
double dfSemiMajor,double dfInvFlattening,
const char * pszPMName ,
double dfPMOffset ,
const char * pszUnits,
double dfConvertToRadians);
double OSRGetSemiMajor( OGRSpatialReferenceH, OGRErr * );
double OSRGetSemiMinor( OGRSpatialReferenceH, OGRErr * );
double OSRGetInvFlattening( OGRSpatialReferenceH, OGRErr * );
OGRErr OSRSetAuthority( OGRSpatialReferenceH hSRS,
const char *pszTargetKey,
const char *pszAuthority,
int nCode );
OGRErr OSRSetProjParm( OGRSpatialReferenceH, const char *, double );
double OSRGetProjParm( OGRSpatialReferenceH hSRS,
const char *pszParmName,
double dfDefault,
OGRErr * );
OGRErr OSRSetUTM( OGRSpatialReferenceH hSRS, int nZone, int bNorth );
int OSRGetUTMZone( OGRSpatialReferenceH hSRS, int *pbNorth );
OGRCoordinateTransformationH
OCTNewCoordinateTransformation( OGRSpatialReferenceH hSourceSRS,
OGRSpatialReferenceHhTargetSRS );
void OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH );
int OCTTransform(OGRCoordinateTransformationH hCT,
int nCount, double *x, double*y, double *z );

2、Python语言接口

[python] view plaincopyprint?

class osr.SpatialReference
def __init__(self,obj=None):
def ImportFromWkt( self, wkt ):
def ExportToWkt(self):
def ImportFromEPSG(self,code):
def IsGeographic(self):
def IsProjected(self):
def GetAttrValue(self, name, child = 0):
def SetAttrValue(self, name, value):
def SetWellKnownGeogCS(self, name):
def SetProjCS(self, name = “unnamed” ):
def IsSameGeogCS(self, other):
def IsSame(self, other):
def SetLinearUnits(self, units_name, to_meters ):
def SetUTM(self, zone, is_north = 1):
class CoordinateTransformation:
def __init__(self,source,target):
def TransformPoint(self, x, y, z = 0):
def TransformPoints(self, points):

七、内部实现

OGRCoordinateTransformation依赖于PROJ.4库。所以要使用坐标转换的内容，GDAL必须在编译的时候绑定PROJ4才可以用来进行坐标转换。如果要使用GDAL的坐标转换，重投影相关的算法，就必须要有PROJ4库的支持，否则会转换失败。关于PROJ4的编译和GDAL绑定PROJ4的内容，请参考之前的博文。

GDAL源码剖析（十二）之GDAL Warp API使用说明

一、简介

本文原文地址：http://www.gdal.org/warptut.html

GDAL Warp API（在文件gdalwarper.h中定义）是一个高效的进行图像变换的接口。主要由几何变换函数（GDALTransformerFunc），多种图像重采样方式，掩码操作选项等组成。这个接口可以对很大的图像进行处理。

下面说明示例让你如何在程序中使用变换API。首先假定你已经熟悉了GDAL的抽象数据模型，以及GDAL的API。

在程序中，首先要初始化一个GDALWarpOptions 结构体的对象，然后使用GDALWarpOptions 的对象来初始化GDALWarpOperation 的对象，最后通过调用GDALWarpKernel 类里面的GDALWarpOperation::ChunkAndWarpImage() 函数来完成图像的变换。

二、简单的影像重投影示例

首先我们以一个图像重投影的例子来入手，需要注意的是，这里假设输出结果文件已经存在，同时这里没有对错误信息进行检查，仅仅演示的最正常的处理流程。

[cpp] view plaincopyprint?

#include “gdalwarper.h”
int main()
{
GDALDatasetH hSrcDS, hDstDS;
// 打开输入输出图像
GDALAllRegister();
hSrcDS = GDALOpen(“in.tif”, GA_ReadOnly );
hDstDS = GDALOpen(“out.tif”, GA_Update );
// 建立变换选项
GDALWarpOptions*psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS =hSrcDS;
psWarpOptions->hDstDS =hDstDS;
psWarpOptions->nBandCount = 1;
psWarpOptions->panSrcBands =
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
psWarpOptions->panSrcBands[0] = 1;
psWarpOptions->panDstBands =
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->pfnProgress = GDALTermProgress;
// 创建重投影变换函数
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( hSrcDS,
GDALGetProjectionRef(hSrcDS),
hDstDS,
GDALGetProjectionRef(hDstDS),
FALSE,0.0, 1 );
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// 初始化并且执行变换操作
GDALWarpOperationoOperation;
oOperation.Initialize(psWarpOptions );
oOperation.ChunkAndWarpImage( 0, 0,
GDALGetRasterXSize( hDstDS),
GDALGetRasterYSize( hDstDS) );
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
return 0;
}

这个例子中首先打开已经存在的输入文件和输出文件（in.tif和out.tif）。使用GDALCreateWarpOptions()函数创建一个GDALWarpOptions结构体对象（结构体中的参数会指定一个比较合理的默认值），然后对这个结构体对象设置输入输出文件的句柄（就是文件指针）和要处理的波段。需要注意的是panSrcBands和panDstBands数组需要在外面动态申请，然后在调用GDALDestroyWarpOptions()函数的时候会自动释放，所以在后面就不需要我们对这两个数组进行释放了。GDALTermProgress是一个简单的控制台进度信息函数，用来显示处理进度。

GDALCreateGenImgProjTransformer()函数是用来创建一个从原始图像到结果图像的投影变换参数。我们假设这两个图像都有合适的四至范围和坐标系统，不使用GCP点。

一旦GDALWarpOptions结构体创建好了，可以使用这个GDALWarpOptions对象来初始化一个GDALWarpOperation的对象，然后再调用函数GDALWarpOperation::ChunkAndWarpImage()进行转换。在转换完成之后，转换选项，数据集等都需要进行释放。

通常应该在打开图像之后进行一系列的检查，然后在建立投影变换参数是要对返回的参数进行检查（返回值为NULL表示失败），最后还要对建立的变换操作进行检查。上面的例子为了方便，没有对这些信息进行检查，在我们自己的程序中需要对这些进行检查。

三、其他变换选项

GDALWarpOptions结构体中包含了很多参数来对变换进行设置。下面对一些比较重要的进行列举说明：

1. GDALWarpOptions::dfWarpMemoryLimit：设置GDALWarpOperation在处理图像中使用的最大内存数。单位为比特，默认值比较保守（比较小），可以根据自己的内存大小来进行调整这个值。增加这个值可以帮助提高程序的运行效率，但是需要注意内存的大小。这个大小、GDAL的缓存大小，还有你的应用程序以及系统所需要的内存加起来要小于系统的内存，否则会导致错误。一般来说，比如一个内存为256MB的系统，这个值最少设置为64MB比较合理。注意，这个值不包括GDAL读取数据使用的内存。

2. GDALWarpOpations::eResampleAlg：重采样方式，可用的值有 GRA_NearestNeighbour （最邻近采样，默认值，处理速度最快）、GRA_Bilinear（2×2双线性内插采样）和 GRA_Cubic（三次立方卷积采样）。GRA_NearestNeighbour采样方式一般用在专题图或者彩色图像中，其他的重采样方式效果比较好，尤其是在计算中改变图像分辨率的算法中。

3. GDALWarpOptions::padfSrcNoDataReal：这个数组（每个波段一个值），可以用来指定输入图像波段的NODATA值，比如图像的背景值，对于这样的值，算法不会参与运算，直接将这个值复制到结果图像中。

4. GDALWarpOptions::papszWarpOptions：这个是一个字符串列表，用来设置图像变换过程中的一些选项，样子为NAME=VALUE。更多详细的内容可以参考GDALWarpOptions::papszWarpOptions的文档，里面含有全部的选项，支持的值包括：

INIT_DEST=[value]或者INIT_DEST=NO_DATA：这个选项用来强制设置结果图像的初始值（所有的波段），初始值为指定的value，或者NODATA值。NODATA值从参数padfDstNoDataReal或者参数padfDstNoDataImag中获取。如果这个值没有设置，那么将会使用原始图像的NODATA值来覆盖。
WRITE_FLUSH=YES/NO：这个选项用来强制设置在处理完每一块后将数据写入磁盘中。在某些时候，这个选项可以更加安全的写入结果数据，但是同时会增加更多的磁盘操作。目前这个默认值为NO。

四、创建输出文件

在前面的例子中，结果图像是已经存在的。选择我们将通过预测输出文件的范围和坐标系统来创建新的图像。这个操作不是图像变换的特殊API，这个仅仅是变换的一个API。

[cpp] view plaincopyprint?

#include “gdalwarper.h”
#include”ogr_spatialref.h”
…
GDALDriverH hDriver;
GDALDataType eDT;
GDALDatasetH hDstDS;
GDALDatasetH hSrcDS;
// 打开源文件
hSrcDS = GDALOpen(“in.tif”, GA_ReadOnly );
CPLAssert( hSrcDS != NULL );
// 创建输出图像的数据类型和输入图像第一个波段类型一致
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
// 获取输出图像的驱动（GeoTIFF格式）
hDriver = GDALGetDriverByName(“GTiff” );
CPLAssert( hDriver != NULL );
// 获取源图像的坐标系统
const char *pszSrcWKT, *pszDstWKT = NULL;
pszSrcWKT = GDALGetProjectionRef( hSrcDS);
CPLAssert( pszSrcWKT != NULL &&strlen(pszSrcWKT) > 0 );
// 设置输出图像的坐标系统为UTM 11 WGS84
OGRSpatialReference oSRS;
oSRS.SetUTM( 11, TRUE );
oSRS.SetWellKnownGeogCS( “WGS84”);
oSRS.exportToWkt( &pszDstWKT );
// 创建一个变换参数，从源图像的行列号坐标到结果图像的地理坐标（没有
//结果行列）的句柄，初始值设置为NULL。
void *hTransformArg;
hTransformArg =
GDALCreateGenImgProjTransformer( hSrcDS,pszSrcWKT, NULL, pszDstWKT,
FALSE,0, 1 );
CPLAssert( hTransformArg != NULL );
// 获取输出文件的近似地理范围和分辨率
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
CPLErr eErr;
eErr = GDALSuggestedWarpOutput( hSrcDS,
GDALGenImgProjTransform, hTransformArg,
adfDstGeoTransform, &nPixels, &nLines );
CPLAssert( eErr == CE_None );
GDALDestroyGenImgProjTransformer(hTransformArg );
// 创建输出文件
hDstDS = GDALCreate(hDriver, “out.tif”, nPixels, nLines,
GDALGetRasterCount(hSrcDS),eDT, NULL );
CPLAssert( hDstDS != NULL );
// 写入投影
GDALSetProjection( hDstDS,pszDstWKT );
GDALSetGeoTransform( hDstDS,adfDstGeoTransform );
// 复制颜色表，如果有的话
GDALColorTableH hCT;
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1));
if( hCT != NULL )
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1),hCT );
… 变换处理之前做的工作…

这里需要注意的一些逻辑关系：

我们需要创建的输出坐标是使用地理坐标表示的，而不是行列号坐标。同时输出的坐标是在结果坐标系统之下的坐标，不是原始坐标系统下的。
GDALSuggestedWarpOutput()函数返回值adfDstGeoTransform、nPixels 和nLines这三个值是描述输出图像的大小和四至范围，这个四至范围包含了源图像的所有的像素点。里面的分辨率是根据原始图像估算出来的，输出图像的分辨率一致是横向和纵向保持一致，不受输入图像限制。
变化要求输出文件格式是可以“随机”写入的。一般使用Create()（不能用CreateCopy()函数）函数创建非压缩的图像格式。如果要创建压缩的图像格式，或者使用CreateCopy()函数创建的话，系统内部会创建一个临时文件，然后使用CreateCopy()函数写入到结果图像中。
变换API仅仅处理图像的象元值。所有的颜色表，地理参考以及其他元数据信息必须使用程序自行设置到结果图像中去，变换是不会设置这些信息的。

五、性能优化

下面几点可以在使用变换API的时候提高处理效率。

增加变换API使用的内存，可以使程序执行的更快。这个参数叫GDALWarpOptions::dfWarpMemoryLimit。理论上，越大的块对于数据I/O效率越高，并且变换的效率也会越高。然而，变换所需的内存和GDAL缓存应该小于机器的内存，最多是内存的2/3左右。
增加GDAL的缓存。这个尤其对于在处理非常大的输入图像很有用。如果所有的输入输出图像的数据频繁的读取会极大的降低运行效率。使用函数GDALSetCacheMax()来设置GDAL使用的最大缓存。
使用近似的变换代替精确的变换（精确的变换是每个象元都会计算一次）。下面的代码用来说明近似变换的使用方式，近似变换要求指定一个变换的误差dfErrorThreshold，这个误差单位是输出图像的象元个数。
[plain] view plaincopyprint?
1. <span style=”font-size:16px;color:#000099;”>hTransformArg =
2. GDALCreateApproxTransformer(GDALGenImgProjTransform,
3. hGenImgProjArg, dfErrorThreshold );
4. pfnTransformer = GDALApproxTransform;
5. </span>
当写入一个空的输出文件，在GDALWarpOptions::papszWarpOptions 对象中使用INIT_DEST选项来进行设置。这样可以很大的减少磁盘的IO操作。
输入和输出文件使用分块存储的文件格式。分块文件格式可以快速的访问一块数据，相比来说，普通的大数据量的顺序存储文件格式在IO操作中要花费更多的时间。
在一次调用中处理所有的波段。一次处理所有的波段比每次处理一个波段要保险。
使用GDALWarpOperation::ChunkAndWarpMulti()方法来代替GDALWarpOperation::ChunkAndWarpImage()方法。这个使用多个独立的线程来进行IO操作和变换影像处理，能够提高CPU和IO的效率。执行这个操作，GDAL需要多线程的支持（从GDAL1.8.0开始，默认在Win32平台、Unix平台是支持的，对于之前的版本，需要在配置中进行修改才行）。
重采样方式，最邻近象元执行速度最快，双线性内插次之，三次立方卷次最慢。不要使用根据复杂的重采样函数。
避免使用复杂的掩码选项。比如，最邻近采样在处理没有掩码的8bit数据要比一般的效率高很多。

六、关于掩码选项

GDALWarpOptions包含了一些处理复杂的掩码的能力，比如掩码的有效性，对输入和输出数据的掩码。这些功能还没有做足够的测试。其他每个波段的有效的掩码在使用的时候要小心。

GDAL源码剖析（十三）之GDAL网格插值说明

一、简介

英文网址：http://www.gdal.org/grid_tutorial.html。

网格插值的意思就是从离散的数据点创建一个栅格图像的过程。通常情况下，你有一系列研究区域的离散点，如果你想将这些点转换为规则的网格数据来进行进一步的处理，或者和其他网格数据进行合并等处理。下图是网格插值的一个示意图：

网格插值示意图

使用数据插值和逼近算法可以用来解决这个问题，但是插值的作用并不仅仅用来处理这个问题。有时候，你并不需要对数据进行插值处理，你需要做的就是计算一下数据覆盖区域的统计值和其他指标。统计值本身是十分有用的，或者你可以根据这些来选择更好的插值算法和参数。

这就是GDAL网格插值API所要理解的东西。它可以帮助你更好的对数据进行插值（参考下面的离散数据插值）或者计算数据的指标信息（参加下面的数据指标计算）。

一共有两种方式来使用这个接口。一种是使用GDALGridCreate的C语言接口编程实现，另一种是使用gdal_grid工具来使用。本文后面的内容将详细介绍GDAL网格插值API的算法原理及其参数表示的含义。

二、离散数据插值

1、反距离权重插值（Inverse Distanceto a Power）

反距离权重插值方法是一种加权平均插值方式。你应该提供离散的数据值和每个点的坐标信息以及输出的格网，然后函数会插值计算输出格网节点的数据值。每一个格网节点的计算方式如下：

上式中字母的含义分别是：

Zi是已知点i的值；
r是格网节点到点i的距离；
p是权重指数；
n是搜索椭圆中的点个数。

在这个算法中权重系数ω的计算方式是：

可以参考函数GDALGridCreate中结构体GDALGridInverseDistanceToAPowerOptions的参数列表和gdal_grid工具的invdist选项列表。

2、移动平均值（Moving Average）

移动平均值是一个简单的数据平均算法，具体是用一个椭圆形的移动窗口，然后搜索在移动窗口中的所有的离散点，然后计算平均值。搜索椭圆可以指定旋转角度，椭圆的中心点位于网格节点。数据点的最少个数的平均值可以进行设置，如果没有足够的点在移动窗口中，这个网格节点就是空的，然后使用NODATA值来进行填充。

算法的数学公式可用下面的公式来表示：

上式中字母的含义分别是：

Z是插值结果，
Zi是已知点i的值；
n是搜索椭圆中的点个数。

可以参考函数GDALGridCreate中结构体GDALGridMovingAverageOptions的参数列表和gdal_grid 工具的average选项列表。

3、最邻近插值（Nearest Neighbor）

最邻近插值不使用任何插值算法和平滑算法，只需在搜索椭圆中找到离中心网格节点最近的离散点，然后把该离散点的值作为网格的节点值。如果没有找到点，将该格网节点设置为NODATA值。

可以参考函数GDALGridCreate中结构体GDALGridNearestNeighborOptions的参数列表和gdal_grid工具的nearest选项列表。

三、数据指标计算

所有的指标计算都是用一个叫GDALGridDataMetricsOptions的结构体来进行控制。

1、数据最小值（Minimum Data Value）

在搜索椭圆中找到的最小值，如果没有找到点，将返回NODATA值，公式如下：

式中：

Z是返回的结果值；
Zi表示第i个点的值；
N表示在搜索椭圆中的点。

2、数据最大值（Maximum Data Value）

在搜索椭圆中找到的最大值，如果没有找到点，将返回NODATA值，公式如下：

式中：

Z是返回的结果值；
Zi表示第i个点的值；
N表示在搜索椭圆中的点。

3、数据范围（Data Range）

在搜索椭圆中最大值和最小值的差，如果没有点，返回NODATA值，公式如下：

式中：

Z是返回的结果值；

Zi表示第i个点的值；

N表示在搜索椭圆中的点。

4、搜索椭圆（Search Ellipse）

在格网插值算法中用到的搜索窗口是一个旋转的椭圆，可以用下面三个参数来进行描述：

第一半径（如果旋转角度为0，就是x轴）
第二半径（如果旋转角度为0，就是y轴）
旋转角度（逆时针方向）

只有点位于搜索椭圆之内（包括边界上）的点才用于计算。

本文转自http://blog.csdn.net/liminlu0314窗体顶端

转载自：https://blog.csdn.net/startwithdp/article/details/45563691

GDAL源码剖析

相关推荐

You may also like...

公众号

微信群

微信小程序

GDAL源码剖析

相关推荐

You may also like...

GDAL地理与投影坐标系

Python空间数据处理2： GDAL栅格图像格式转换

Gdal关于CAD转SHP格式

公众号

微信群

微信小程序