给定一张TIF格式的影像，给定一套格网重投影坐标，输出一张新的TIF影像（要求是分块存储格式Tiled TIF）; 实现为C++动态库，输出函数的接口定义如下： imgPath是输入的TIF影像路径; outImgPath是输出的TIF影像路径，要求以Tiled TIF格式存储，分块大小为256; gridWidth是原影像的格网矫正宽度，通常为64; gridPoints是原影像的格网矫正坐标，此数组中坐标点的顺序是逐行顺序存储。假设影像宽高分别为W，H，那么gridPoints中的点数为ceil（W/64）*H; 由于矫正坐标的最小值不一定为（0，0），因而输出影像是以矫正坐标的覆盖范围来决定影像的宽高，同时要求将矫正坐标的最小值作为偏移坐标值写入到TIF数据的GeoTransform头信息中。功能要求： 影像可能较大，因而不能完全加载全部影像到内存进行处理，所以只能分块处理; 速度较快，一张1GB的TIF数据，处理时间不能超过20秒; 自编一个程序测试上述接口;3）环境要求： C++， GDAL3.0.2，VS 2017 推荐使用GDAL3.0.2读写影像

实现步骤：

1. 打开输入影像并读取地理参考信息、影像大小等信息；
2. 计算输出影像的大小和偏移坐标值；
3. 创建输出影像并设置其地理参考信息、分块大小等参数；
4. 分块读取输入影像数据，进行重投影，并写入输出影像中。

接口定义如下：

bool reprojectTif(const char* imgPath, const char* outImgPath, double gridWidth, double* gridPoints);

其中，imgPath为输入影像路径，outImgPath为输出影像路径，gridWidth为原影像格网矫正宽度，gridPoints为原影像格网矫正坐标数组。

下面是实现代码：

#include "gdal_priv.h"
#include "cpl_conv.h" // for CPLMalloc()

bool reprojectTif(const char* imgPath, const char* outImgPath, double gridWidth, double* gridPoints)
{
    GDALDataset* poSrcDS;
    GDALDataset* poDstDS;
    GDALDataType eDT;
    int nBands;

    // 打开输入影像
    poSrcDS = (GDALDataset*)GDALOpen(imgPath, GA_ReadOnly);
    if (poSrcDS == NULL) {
        return false;
    }

    // 读取输入影像信息
    double adfGeoTransform[6];
    poSrcDS->GetGeoTransform(adfGeoTransform);

    int nXSize = poSrcDS->GetRasterXSize();
    int nYSize = poSrcDS->GetRasterYSize();

    eDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
    nBands = poSrcDS->GetRasterCount();

    // 计算输出影像大小和偏移坐标值
    double minX = adfGeoTransform[0] + gridWidth * (int)(gridPoints[0] / gridWidth);
    double minY = adfGeoTransform[3] + gridWidth * (int)(gridPoints[1] / gridWidth);
    double maxX = adfGeoTransform[0] + gridWidth * ((int)(gridPoints[(nXSize / 64) * 2 - 2] / gridWidth) + 1);
    double maxY = adfGeoTransform[3] + gridWidth * ((int)(gridPoints[nYSize * 2 - 1] / gridWidth) + 1);

    int nOutXSize = (int)((maxX - minX) / adfGeoTransform[1] + 0.5);
    int nOutYSize = (int)((maxY - minY) / fabs(adfGeoTransform[5]) + 0.5);

    // 创建输出影像
    char** papszOptions = NULL;
    papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES");
    papszOptions = CSLSetNameValue(papszOptions, "BLOCKXSIZE", "256");
    papszOptions = CSLSetNameValue(papszOptions, "BLOCKYSIZE", "256");

    poDstDS = (GDALDataset*)GDALCreateCopy(GDALGetDriverByName("GTiff"), outImgPath, poSrcDS, FALSE, papszOptions, NULL, NULL);
    if (poDstDS == NULL) {
        GDALClose(poSrcDS);
        return false;
    }

    // 设置输出影像地理参考信息
    double adfGeoTransformOut[6] = { minX, adfGeoTransform[1], 0.0, maxY, 0.0, adfGeoTransform[5] };
    poDstDS->SetGeoTransform(adfGeoTransformOut);

    // 分块处理影像
    int nXBlocks = (nOutXSize + 255) / 256;
    int nYBlocks = (nOutYSize + 255) / 256;
    int nBlockXSize = 256;
    int nBlockYSize = 256;
    int bufSize = nBlockXSize * nBlockYSize * nBands;

    void* pBuf = CPLMalloc(bufSize);

    for (int i = 0; i < nYBlocks; i++) {
        for (int j = 0; j < nXBlocks; j++) {
            int nXOff = j * nBlockXSize;
            int nYOff = i * nBlockYSize;

            if (nXOff + nBlockXSize > nOutXSize) {
                nBlockXSize = nOutXSize - nXOff;
            }

            if (nYOff + nBlockYSize > nOutYSize) {
                nBlockYSize = nOutYSize - nYOff;
            }

            double* pProjX = (double*)CPLMalloc(nBlockXSize * sizeof(double));
            double* pProjY = (double*)CPLMalloc(nBlockYSize * sizeof(double));

            // 计算当前分块的矩形范围
            double xMin = adfGeoTransformOut[0] + j * 256 * adfGeoTransformOut[1];
            double yMax = adfGeoTransformOut[3] - i * 256 * fabs(adfGeoTransformOut[5]);
            double xMax = adfGeoTransformOut[0] + (j + 1) * 256 * adfGeoTransformOut[1];
            double yMin = adfGeoTransformOut[3] - (i + 1) * 256 * fabs(adfGeoTransformOut[5]);

            // 重投影
            for (int k = 0; k < nBlockXSize; k++) {
                pProjX[k] = xMin + k * adfGeoTransformOut[1];
            }

            for (int k = 0; k < nBlockYSize; k++) {
                pProjY[k] = yMax - k * fabs(adfGeoTransformOut[5]);
            }

            GDALRasterIO(poSrcDS->GetRasterBand(1), GF_Read, (int)((pProjX[0] - adfGeoTransform[0]) / adfGeoTransform[1] + 0.5), (int)((pProjY[0] - adfGeoTransform[3]) / adfGeoTransform[5] + 0.5), nBlockXSize, nBlockYSize, pBuf, nBlockXSize, nBlockYSize, eDT, nBands, NULL, 0, 0, 0);

            for (int k = 0; k < nBlockXSize * nBlockYSize; k++) {
                double x = adfGeoTransform[0] + (nXOff + k % nBlockXSize) * adfGeoTransform[1];
                double y = adfGeoTransform[3] - (nYOff + k / nBlockXSize) * fabs(adfGeoTransform[5]);

                double dfGeoX, dfGeoY;
                GDALApplyGeoTransform(adfGeoTransformOut, x, y, &dfGeoX, &dfGeoY);

                int nProjX = (int)((dfGeoX - adfGeoTransform[0]) / adfGeoTransform[1] + 0.5);
                int nProjY = (int)((dfGeoY - adfGeoTransform[3]) / adfGeoTransform[5] + 0.5);

                GDALRasterIO(poDstDS->GetRasterBand(1), GF_Write, nProjX, nProjY, 1, 1, (char*)pBuf + k * nBands * GDALGetDataTypeSize(eDT) / 8, 1, 1, eDT, nBands, NULL, bufSize, 1, GDALGetDataTypeSize(eDT) / 8);
            }

            CPLFree(pProjX);
            CPLFree(pProjY);
        }
    }

    CPLFree(pBuf);

    // 关闭数据集
    GDALClose(poSrcDS);
    GDALClose(poDstDS);

    return true;
}

测试程序如下：

#include <iostream>
#include "gdal_priv.h"

using namespace std;

bool reprojectTif(const char* imgPath, const char* outImgPath, double gridWidth, double* gridPoints);

int main()
{
    GDALAllRegister();

    double gridPoints[] = { 0, 0, 64, 0, 128, 0, 192, 0, 256, 0, 320, 0, 384, 0, 448, 0, 0, 64, 64, 64, 128, 64, 192, 64, 256, 64, 320, 64, 384, 64, 448, 64, 0, 128, 64, 128, 128, 128, 192, 128, 256, 128, 320, 128, 384, 128, 448, 128, 0, 192, 64, 192, 128, 192, 192, 192, 256, 192, 320, 192, 384, 192, 448, 192, 0, 256, 64, 256, 128, 256, 192, 256, 256, 256, 320, 256, 384, 256, 448, 256, 0, 320, 64, 320, 128, 320, 192, 320, 256, 320, 320, 320, 384, 320, 448, 320, 0, 384, 64, 384, 128, 384, 192, 384, 256, 384, 320, 384, 384, 384, 448, 384, 0, 448, 64, 448, 128, 448, 192, 448, 256, 448, 320, 448, 384, 448, 448 };
    double gridWidth = 64.0;

    reprojectTif("input.tif", "output.tif", gridWidth, gridPoints);

    return 0;
}