GDAL使用（三）
GDALWarp的功能比较强大，可用于影像投影/重投影、图像镶嵌、重采样、规则裁切、图像校正等。
GDALWarp API定义在头文件gdalwarper.h中。主要有两部分构成：
GDALWarpOptions和 GDALWarpOperation
如字面意思，GDALWarpOperation主要用来执行变换操作，GDALWarpOptions主要是对变换的一些选项进行设置。
1、GDALWarpOptions
这是个结构体，在头文件中的定义如下：
点击查看结构体
/************************************************************************/
/* GDALWarpOptions */
/************************************************************************/
/** Warp control options for use with GDALWarpOperation::Initialize() */
typedef struct {
char **papszWarpOptions;
/*! In bytes, 0.0 for internal default */
double dfWarpMemoryLimit;
/*! Resampling algorithm to use */
GDALResampleAlg eResampleAlg;
/*! data type to use during warp operation, GDT_Unknown lets the algorithm
select the type */
GDALDataType eWorkingDataType;
/*! Source image dataset. */
GDALDatasetH hSrcDS;
/*! Destination image dataset - may be NULL if only using GDALWarpOperation::WarpRegionToBuffer(). */
GDALDatasetH hDstDS;
/*! Number of bands to process, may be 0 to select all bands. */
int nBandCount;
/*! The band numbers for the source bands to process (1 based) */
int *panSrcBands;
/*! The band numbers for the destination bands to process (1 based) */
int *panDstBands;
/*! The source band so use as an alpha (transparency) value, 0=disabled */
int nSrcAlphaBand;
/*! The dest. band so use as an alpha (transparency) value, 0=disabled */
int nDstAlphaBand;
/*! The "nodata" value real component for each input band, if NULL there isn't one */
double *padfSrcNoDataReal;
/*! The "nodata" value imaginary component - may be NULL even if real
component is provided. This value is not used to flag invalid values.
Only the real component is used. */
double *padfSrcNoDataImag;
/*! The "nodata" value real component for each output band, if NULL there isn't one */
double *padfDstNoDataReal;
/*! The "nodata" value imaginary component - may be NULL even if real
component is provided. Note that warp operations only use real component
for flagging invalid data.*/
double *padfDstNoDataImag;
/*! GDALProgressFunc() compatible progress reporting function, or NULL
if there isn't one. */
GDALProgressFunc pfnProgress;
/*! Callback argument to be passed to pfnProgress. */
void *pProgressArg;
/*! Type of spatial point transformer function */
GDALTransformerFunc pfnTransformer;
/*! Handle to image transformer setup structure */
void *pTransformerArg;
/** Unused. Must be NULL */
GDALMaskFunc *papfnSrcPerBandValidityMaskFunc;
/** Unused. Must be NULL */
void **papSrcPerBandValidityMaskFuncArg;
/** Unused. Must be NULL */
GDALMaskFunc pfnSrcValidityMaskFunc;
/** Unused. Must be NULL */
void *pSrcValidityMaskFuncArg;
/** Unused. Must be NULL */
GDALMaskFunc pfnSrcDensityMaskFunc;
/** Unused. Must be NULL */
void *pSrcDensityMaskFuncArg;
/** Unused. Must be NULL */
GDALMaskFunc pfnDstDensityMaskFunc;
/** Unused. Must be NULL */
void *pDstDensityMaskFuncArg;
/** Unused. Must be NULL */
GDALMaskFunc pfnDstValidityMaskFunc;
/** Unused. Must be NULL */
void *pDstValidityMaskFuncArg;
/** Unused. Must be NULL */
CPLErr (*pfnPreWarpChunkProcessor)( void *pKern, void *pArg );
/** Unused. Must be NULL */
void *pPreWarpProcessorArg;
/** Unused. Must be NULL */
CPLErr (*pfnPostWarpChunkProcessor)( void *pKern, void *pArg);
/** Unused. Must be NULL */
void *pPostWarpProcessorArg;
/*! Optional OGRPolygonH for a masking cutline. */
void *hCutline;
/*! Optional blending distance to apply across cutline in pixels, default is zero. */
double dfCutlineBlendDist;
} GDALWarpOptions;下面对结构体的成员进行说明：
papszWarpOptions
char** 类型，控制warp操作，是格式为“name=value”的字符串数组，可以使用CSLSetNameValue()来构造该字符数组。选项的具体说明见如下链接：https://www.osgeo.cn/gdal/api/gdalwarp_cpp.html#_CPPv4N15GDALWarpOptions16papszWarpOptionsE
dfMemoryLimit
double 类型，用于指定图像变换时使用的内存大小，单位为字节，为0.0时系统内部自动分配
eResampleAlg
GDALResampleAlg 类型，指示变换时的重采样方法，包含了大部分常用的重采样方法，具体方法见以下链接：
https://www.osgeo.cn/gdal/api/gdalwarp_cpp.html?highlight=gdalwarpoptions#_CPPv415GDALResampleAlg
eWorkingDataType
GDALDataType 类型，指示在图像变换时使用的数据类型。设置为GDT_Unknown时算法自动选择数据类型。类型的具体说明见以下链接：
https://www.osgeo.cn/gdal/api/raster_c_api.html#_CPPv412GDALDataType
hSrcDS
GDALDatasetH 格式，原始影像数据集。
hDstDS
GDALDatasetH格式，目标数据集。
nBandCount
int 类型，要处理的波段数目，如果为0，将处理所有波段。
panSrcBands
int* 类型，指定要处理的原始图像的波段序号，序号从1开始。
panDstBands
int* 类型，指定要处理的原结果图像的波段序号，序号从1开始。
nSrcAlphaBand
int 类型，指定原始图像所 使用的Alpha 通道（ 透明 通道）序号， 默认值为0，表示 没有。
nDstAlphaBand
int 类型，指定结果图像所使用的 Alpha 通道（ 透明 通道）序号，默认值 为 0，表示没有。
padfSrcNoDataReal
double* 类型，该变量指定原始图像中实部的NoData值，NULL表示没有（对于复数图像来说）。
padfSrcNoDataImag
double* 类型，该变量指定原始图像中虚部的NoData值，NULL表示没有（对于复数图像来说）。
padfDstNoDataReal
double* 类型，该变量指定结果图像中实部的NoData值，NULL表示没有（对于复数图像来说）。
padfDstNoDataImag
double* 类型，该变量指定结果图像中虚部的NoData值，NULL表示没有（对于复数图像来说）。
pfnProgress
GDALProgressFunc 类型，进度条信息函数指针
pProgressArg
void* 类型，进度条函数的参数指针。
pfnTransformer
GDALTransformerFunc 类型，图像坐标变换函数指针，这是核心函数，用来确定图像变换的方式。
pTransformerArg
void* 类型，图像坐标变换函数所需要的结构体指针参数。
hCutline
void* 类型，该变量使用OGRPolygonH指定 一个裁切范围，作为掩码区域。默认为 NULL，表示 处理全图。
dfCutlineBlendDist
double 类型，该变量指定在使用hCutline作为掩码区域时的羽化距离。默认为0。
最后，这个结构体指针可以使用GDALCreateWarpOption()函数进行创建，这个函数的返回值直接就是GDALWarpOptions*类型。
2、GDALWarpOperation
是负责对图像进行读取、分块和写入操作的类。
函数：Initialize
原型：CPLErr Initialize(constGDALWarpOptions *psNewOptions)
参数：psNewOptions：constGDALWarpOptions *类型，就是上面一节说的结构体指针。
返回值：是否初始化成功
函数：GetOptions
原型：const GDALWarpOptions *GetOptions()
返回值：返回上一节说的结构体指针
函数：ChunkAndWarpImage
原型：CPLErr ChunkAndWarpImage(int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize)
参数：nDstXOff：结果数据X起始坐标，字面意思就是X方向的偏移量，nDstYOff：结果数据Y起始坐标，nDstXSize：结果数据宽度，nDstYSize：结果数据高度
返回值：是否执行成功
函数：ChunkAndWarpMulti
原型：CPLErr ChunkAndWarpMulti(int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize)
参数与返回值与第三条中的一样，功能也相同，不过内部是多线程操作的
其它函数参考下面链接：
https://www.osgeo.cn/gdal/api/gdalwarp_cpp.html#classGDALWarpOperation
3、重投影代码
这里的例子是将一幅未投影影像投影到投影坐标下。
原始数据：海洋一号C卫星（HY-1C）L1C数据，此原始数据为WGS84坐标系，一共4个波段。
目标数据：这里目标是将原始数据投影到投影坐标系下，目标空间参考为 Asia North Albers Equal Area Conic，此坐标系为WGS84椭球下的等积圆锥投影，且数据的分辨率约为50m，这里让GDAL自己计算合适的分辨率。
#include <iostream>
#include "gdal_priv.h"
#include "cpl_conv.h"
#include <cpl_string.h>
#include <gdalwarper.h>
int main()
{
clock_t start = clock();
GDALAllRegister();
GDALDriverH hDriver;
GDALDataType eDT;
GDALDatasetH hDstDS;
GDALDatasetH hSrcDS;
// 打开源文件
hSrcDS = GDALOpen("E:\testdata\CZIOrigin.tiff", GA_ReadOnly);
CPLAssert(hSrcDS != NULL);
// 获取数据类型
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
// 获取输入数据的栅格驱动 (GeoTIFF格式)
hDriver = GDALGetDriverByName("GTiff");
CPLAssert(hDriver != NULL);
// 获取源数据坐标系统
const char* pszSrcWKT;
char* pszDstWKT = NULL;
pszSrcWKT = GDALGetProjectionRef(hSrcDS);
CPLAssert(pszSrcWKT != NULL && strlen(pszSrcWKT) > 0);
// 设置输出坐标系为
OGRSpatialReference oSRS;
oSRS.importFromProj4("+proj=aea +lat_1=15 +lat_2=65 +lat_0=30 +lon_0=95 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs");
oSRS.exportToWkt(&pszDstWKT);
printf("%s
", 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 = 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);
// 创建输出文件
hDstDS = GDALCreate(hDriver, "O:\CZIproject.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);
//设置warp选项
GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;//设置源数据集
psWarpOptions->hDstDS = hDstDS;//设置目的数据集
psWarpOptions->nBandCount = 0;//要处理的波段数，0的话是选择所有波段
psWarpOptions->pfnProgress = GDALTermProgress;//控制台进度条函数
//建立重投影变换
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 0.0, 1);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
//初始化并执行warp操作
GDALWarpOperation oOperation;
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));//这是实际执行warp的函数
//释放资源
GDALDestroyGenImgProjTransformer(hTransformArg);
GDALClose(hDstDS);
GDALClose(hSrcDS);
clock_t end = clock();
double endtime = (double)(end - start) / CLOCKS_PER_SEC;
printf("time elapse: %f
", endtime);
system("pause");
return 0;
}下面的是投影前后的空间参考信息，可以看到GDAL自己计算出的像元大小约为46.2250m
以上部分内容参考：
[1] 李民录. GDAL源码剖析与开发指南（异步图书）. 人民邮电出版社. Kindle 版本.
[2] https://www.osgeo.cn/gdal/index.html
鲁公网安备 37148202000241号