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() 函数来完成图像的变换。

二、简单的影像重投影示例

首先我们以一个图像重投影的例子来入手，需要注意的是，这里假设输出结果文件已经存在，同时这里没有对错误信息进行检查，仅仅演示的最正常的处理流程。

#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（2x2双线性内插采样）和 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。

#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的时候提高处理效率。

1. 增加变换API使用的内存，可以使程序执行的更快。这个参数叫GDALWarpOptions::dfWarpMemoryLimit。理论上，越大的块对于数据I/O效率越高，并且变换的效率也会越高。然而，变换所需的内存和GDAL缓存应该小于机器的内存，最多是内存的2/3左右。
2. 增加GDAL的缓存。这个尤其对于在处理非常大的输入图像很有用。如果所有的输入输出图像的数据频繁的读取会极大的降低运行效率。使用函数GDALSetCacheMax()来设置GDAL使用的最大缓存。
3.  使用近似的变换代替精确的变换（精确的变换是每个象元都会计算一次）。下面的代码用来说明近似变换的使用方式，近似变换要求指定一个变换的误差dfErrorThreshold，这个误差单位是输出图像的象元个数。

4. hTransformArg =
       GDALCreateApproxTransformer(GDALGenImgProjTransform,
   hGenImgProjArg, dfErrorThreshold );
   pfnTransformer = GDALApproxTransform;
   

5. 当写入一个空的输出文件，在GDALWarpOptions::papszWarpOptions 对象中使用INIT_DEST选项来进行设置。这样可以很大的减少磁盘的IO操作。
6. 输入和输出文件使用分块存储的文件格式。分块文件格式可以快速的访问一块数据，相比来说，普通的大数据量的顺序存储文件格式在IO操作中要花费更多的时间。
7. 在一次调用中处理所有的波段。一次处理所有的波段比每次处理一个波段要保险。
8. 使用GDALWarpOperation::ChunkAndWarpMulti()方法来代替GDALWarpOperation::ChunkAndWarpImage()方法。这个使用多个独立的线程来进行IO操作和变换影像处理，能够提高CPU和IO的效率。执行这个操作，GDAL需要多线程的支持（从GDAL1.8.0开始，默认在Win32平台、Unix平台是支持的，对于之前的版本，需要在配置中进行修改才行）。
9. 重采样方式，最邻近象元执行速度最快，双线性内插次之，三次立方卷次最慢。不要使用根据复杂的重采样函数。
10. 避免使用复杂的掩码选项。比如，最邻近采样在处理没有掩码的8bit数据要比一般的效率高很多。

六、关于掩码选项

GDALWarpOptions包含了一些处理复杂的掩码的能力，比如掩码的有效性，对输入和输出数据的掩码。这些功能还没有做足够的测试。其他每个波段的有效的掩码在使用的时候要小心。