GDAL Warp API教程（重新投影…）

概述

GDAL Warp API（在 gdalwarper.h )使用应用程序提供的几何变换函数（GDALTransformerFunc）、各种重采样内核和各种遮罩选项，为高性能图像扭曲提供服务。文件大得多，可以保存在内存中可以扭曲。

本教程演示如何使用Warp API实现应用程序。假定C++中的实现为C，Python绑定对于WARP API来说是不完整的。它还假设您熟悉 栅格数据模型 ，以及通用的GDAL API。

应用程序通常通过初始化 GDALWarpOptions 使用要使用的选项构造，基于这些选项实例化gdalwarp操作，然后调用 GDALWarpOperation::ChunkAndWarpImage() 方法在内部使用 GDALWarpKernel 班级。

一个简单的重审案件

首先，我们将构造一个相对简单的示例，用于重新投影图像，假设已经存在一个适当的输出文件，并且错误检查最少。

#include "gdalwarper.h"

int main()
{
    GDALDatasetH  hSrcDS, hDstDS;

    // Open input and output files.
    GDALAllRegister();
    hSrcDS = GDALOpen( "in.tif", GA_ReadOnly );
    hDstDS = GDALOpen( "out.tif", GA_Update );

    // Setup warp options.
    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;

    // Establish reprojection transformer.
    psWarpOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer( hSrcDS,
                                        GDALGetProjectionRef(hSrcDS),
                                        hDstDS,
                                        GDALGetProjectionRef(hDstDS),
                                        FALSE, 0.0, 1 );
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;
    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）。一个 GDALWarpOptions 结构已分配 (GDALCreateWarpOptions() 为内容设置许多合理的默认值，总是将其用于默认内容），并设置输入和输出文件句柄以及带区列表。panSrcBands和panDstBands列表在此处动态分配，并将通过 GDALDestroyWarpOptions() . 安装了简单终端输出进度监视器（GDALTermProgress），用于向用户报告完成进度。

GDALCreateGenImgProjTransformer() 用于初始化源图像和目标图像之间的重新投影转换。我们假设它们已经有了合理的边界和坐标系集。GCP的使用被禁用。

一旦options结构就绪，就可以使用它们实例化gdalwarp操作，而warp实际上是用 GDALWarpOperation::ChunkAndWarpImage() . 然后清理transformer、warp选项和数据集。

通常在打开文件、设置重新投影转换器（失败时返回空值）和初始化warp之后需要进行错误检查。

其他扭曲选项

gdalwarppoptions结构包含许多可以设置为控制扭曲行为的项。一些特别感兴趣的是：

• GDALWarpOptions::dfWarpMemoryLimit -设置在选择要操作的图像块大小时gdalwarp操作要使用的最大内存量。该值以字节为单位，默认值可能是保守的（小）。在某些情况下，增大块大小可能会有很大帮助，但应注意确保此大小加上GDAL缓存大小加上GDAL的工作集，应用程序和操作系统小于RAM的大小，否则过度交换可能会影响性能。在具有256MB RAM的系统上，至少64MB（大约64000000字节）的值是合理的。请注意，该值不包括GDAL用于低级块缓存的内存。

• GDALWarpOpations::eResampleAlg -GRA_nearest neighbor（默认且最快）、GRA_bilinar（2x2双线性重采样）或GRA_Cubic之一。GRA_nearest neighbor类型通常应用于主题或彩色映射图像。其他的重采样类型可以为专题图像提供更好的结果，特别是在分辨率发生实质性变化时。

• GDALWarpOptions::padfSrcNoDataReal -如果希望避免将某个背景值的像素复制到目标图像，则可以为每个波段设置一个“nodata”值来设置该数组（每个正在处理的波段一个条目）。

• GDALWarpOptions::papszWarpOptions -这是传递给整经器的NAME=VALUE选项的字符串列表。见 GDALWarpOptions::papszWarpOptions 所有选项的文档。支持的值包括：

  • INIT_DEST= [价值] 或in it_DEST=NO_DATA：此选项强制将目标图像初始化为指示值（对于所有频带），或指示应将其初始化为padfdstnodeatareal/padfdstnodeataimag中的NO_DATA值。如果未设置此值，则将读取目标图像并覆盖源扭曲。

  • WRITE_FLUSH=YES/NO：此选项在处理每个数据块后强制对磁盘刷新数据。在某些情况下，这有助于确保输出数据的串行写入，否则每次为输入缓冲区读取数据块时，可能会将数据块写入磁盘，从而导致在磁盘周围进行大量额外的查找，并降低IO吞吐量。此时的默认值是NO。

创建输出文件

在前一种情况下，已假定存在适当的输出文件。现在我们将讨论一个在新坐标系中创建具有适当边界的新文件的情况。此操作与warp API无关。它只是在使用转换API。

#include "gdalwarper.h"
#include "ogr_spatialref.h"
...
GDALDriverH hDriver;
GDALDataType eDT;
GDALDatasetH hDstDS;
GDALDatasetH hSrcDS;

// Open the source file.
hSrcDS = GDALOpen( "in.tif", GA_ReadOnly );
CPLAssert( hSrcDS != NULL );

// Create output with same datatype as first input band.
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

// Get output driver (GeoTIFF format)
hDriver = GDALGetDriverByName( "GTiff" );
CPLAssert( hDriver != NULL );

// Get Source coordinate system.
const char *pszSrcWKT, *pszDstWKT = NULL;
pszSrcWKT = GDALGetProjectionRef( hSrcDS );
CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );

// Setup output coordinate system that is UTM 11 WGS84.
OGRSpatialReference oSRS;
oSRS.SetUTM( 11, TRUE );
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszDstWKT );

// Create a transformer that maps from source pixel/line coordinates
// to destination georeferenced coordinates (not destination
// pixel line).  We do that by omitting the destination dataset
// handle (setting it to NULL).
void *hTransformArg;
hTransformArg =
    GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT,
                                     FALSE, 0, 1 );
CPLAssert( hTransformArg != NULL );

// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
CPLErr eErr;
eErr = GDALSuggestedWarpOutput( hSrcDS,
                                GDALGenImgProjTransform, hTransformArg,
                                adfDstGeoTransform, &nPixels, &nLines );
CPLAssert( eErr == CE_None );
GDALDestroyGenImgProjTransformer( hTransformArg );

// Create the output file.
hDstDS = GDALCreate( hDriver, "out.tif", nPixels, nLines,
                     GDALGetRasterCount(hSrcDS), eDT, NULL );
CPLAssert( hDstDS != NULL );

// Write out the projection definition.
GDALSetProjection( hDstDS, pszDstWKT );
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

// Copy the color table, if required.
GDALColorTableH hCT;
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
if( hCT != NULL )
    GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
... proceed with warp as before ...

关于这个逻辑的一些注释：

• 我们需要创建转换器以输出坐标，这样转换器的输出是地理参考坐标，而不是像素线坐标，因为我们使用转换器将源图像周围的像素映射到目标地理参考坐标。

• 这个 GDALSuggestedWarpOutput() 函数将返回一个adfDstGeoTransform、nPixels和nLines，用于描述输出图像大小和地理参考范围，这些范围应包含源图像中的所有像素。该分辨率旨在与源相当，但无论输入像素的形状如何，输出像素始终是方形的。

• 整经器需要一个输出文件，其格式可以“随机”写入。这通常限制为具有Create（）方法实现的未压缩格式（与CreateCopy（）相反）。要转换为压缩格式或CreateCopy（）样式格式，必须以性能更好的格式生成图像的完整临时副本，然后CreateCopy（）将其转换为所需的最终格式。

• Warp API只复制像素。应用程序必须将所有颜色映射、地理参考和其他元数据复制到目标位置。

性能优化

可以做很多事情来优化warp API的性能：

• 增加可用于Warp API分块的内存量，以便一次可以操作较大的分块。这是 GDALWarpOptions::dfWarpMemoryLimit 参数。理论上，块大小越大，I/O策略的效率就越高，近似转换的效率也就越高。然而，扭曲内存和GDAL缓存的总和应该小于RAM大小，可能是RAM大小的2/3左右。

• 增加GDAL缓存的内存量。这在处理面向扫描线的非常大的输入和输出图像时尤为重要。如果必须为每个块重新读取所有输入或输出扫描线，则它们相交的性能可能会大大降低。使用 GDALSetCacheMax() 控制GDAL中可用于缓存的内存量。

• 对于要变换的每个像素，使用近似变换而不是精确的重投影。此代码说明如何基于重投影变换创建近似变换，但具有给定的错误阈值（输出像素中的dfErrorThreshold）。

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

• 当写入空白输出文件时，请使用 GDALWarpOptions::papszWarpOptions 使输出块初始化为固定值，而不是从输出中读取。这可以大大减少不必要的IO工作。

• 使用平铺输入和输出格式。平铺格式允许访问给定的源图像和目标图像块，而无需接触大量额外的图像数据。面向扫描线的大文件会导致大量额外IO的浪费。

• 在一次呼叫中处理所有波段。这可以确保不必为每个频带执行变换计算。

• 使用 GDALWarpOperation::ChunkAndWarpMulti() 方法而不是 GDALWarpOperation::ChunkAndWarpImage() . 它为IO和实际的映像扭曲操作使用单独的线程，从而更有效地利用CPU和IO带宽。要使其工作，GDAL需要使用多线程支持构建（Win32上的默认值，Unix上的默认值，对于以前的版本，configure中需要使用线程）。

• 重采样核的变化从最近邻最小到双线性再到立方。不要使用比所需更复杂的重采样内核。

• 避免使用深奥的屏蔽选项，以便将特殊的简化逻辑案例用于常见的特殊案例。例如，与一般情况相比，对8bit数据没有掩蔽的近邻重采样是高度优化的。

其他掩蔽选项

gdalwarppoptions包括一系列深奥的屏蔽功能，用于有效性屏蔽，以及输入和输出的密度屏蔽。其中一些还没有实现，其他的已经实现，但是测试很差。除了每个频带的有效性掩码，建议此时小心使用这些功能。