GDAL C++多边形裁剪栅格异常：命令行正常代码生成空文件

问题描述 #

尝试用GDAL通过多边形裁剪栅格时，初始化WarpOperation出现读取访问冲突；终端执行gdalwarp命令可正常完成裁剪，但C++代码运行后输出文件始终为0字节。已知条件：

• 可正常访问Shapefile特征数量及栅格投影信息
• 所有文件采用相同CRS

终端可正常执行的命令：

gdalwarp -of GTiff -cutline polygon.shp -crop_to_cutline input.tif output.tif

出现问题的C++代码：

int main()
{
    GDALAllRegister();
    const char* rasterPath = "input.tif";
    const char* vectorPath = "polygon.shp";
    const char* outputRasterPath = "output.tif";

    GDALDataset* rasterDataset = static_cast<GDALDataset*>(GDALOpen(rasterPath, GA_ReadOnly));
    GDALDataset* vectorDataset = static_cast<GDALDataset*>(GDALOpenEx(vectorPath, GDAL_OF_UPDATE, nullptr, nullptr, nullptr));

    OGRLayer* vectorLayer = vectorDataset->GetLayer(0);

    OGRSpatialReference* rasterSRS = new OGRSpatialReference(rasterDataset->GetProjectionRef());

    OGRSpatialReference* vectorSRS = vectorLayer->GetSpatialRef();
    OGRCoordinateTransformation* transform = OGRCreateCoordinateTransformation(vectorSRS, rasterSRS);

    OGRGeometry* clipGeometry = vectorLayer->GetNextFeature()->GetGeometryRef();

    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* outputDataset = driver->Create(outputRasterPath, rasterDataset->GetRasterXSize(), rasterDataset->GetRasterYSize(), 1, GDT_Float32, nullptr);
    outputDataset->SetProjection(rasterDataset->GetProjectionRef());
    double adGeoTransform[6];
    rasterDataset->GetGeoTransform(adGeoTransform);
    outputDataset->SetGeoTransform(adGeoTransform);

    OGRGeometry* cutlineGeometry = clipGeometry->clone();

    GDALWarpOptions* warpOptions = GDALCreateWarpOptions();
    warpOptions->hSrcDS = rasterDataset;
    warpOptions->hDstDS = outputDataset;
    warpOptions->nBandCount = 1;
    warpOptions->panSrcBands = (int*)CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    warpOptions->panSrcBands[0] = 1;
    warpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    warpOptions->panDstBands[0] = 1;
    warpOptions->papszWarpOptions = CSLDuplicate(nullptr);
    warpOptions->eWorkingDataType = GDT_Float32;
    warpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( rasterDataset, 
                                                                    rasterDataset->GetProjectionRef(), 
                                                                    outputDataset, 
                                                                    outputDataset->GetProjectionRef(), 
                                                                    TRUE, 1000, 1);
    warpOptions->pfnTransformer = GDALGenImgProjTransform;
    warpOptions->pfnProgress = GDALTermProgress;
    warpOptions->hCutline = cutlineGeometry;

    GDALWarpOperation warpOperation;
    warpOperation.Initialize(warpOptions);

    warpOperation.ChunkAndWarpImage(0, 0, outputDataset->GetRasterXSize(), outputDataset->GetRasterYSize());

    TransformCutlineToSource(outputDataset, reinterpret_cast<OGRGeometryH>(cutlineGeometry), &(warpOptions->papszWarpOptions), NULL);

    GDALDestroyGenImgProjTransformer(warpOptions->pTransformerArg);
    GDALDestroyWarpOptions(warpOptions);
    delete cutlineGeometry;
    delete clipGeometry;
    delete transform;
    GDALClose(outputDataset);
    GDALClose(vectorDataset);
    GDALClose(rasterDataset);
    return 0;
}

问题分析与修改方案 #

代码存在多个关键问题，逐一修正后可解决0字节输出和访问冲突：

1. 矢量数据集打开模式错误 #

用GDAL_OF_UPDATE打开Shapefile无必要，改为只读模式即可：

GDALDataset* vectorDataset = static_cast<GDALDataset*>(GDALOpenEx(vectorPath, GDAL_OF_READONLY, nullptr, nullptr, nullptr));

2. 几何对象内存管理错误 #

vectorLayer->GetNextFeature()返回的OGRFeature对象需手动销毁，直接获取的clipGeometry属于特征对象，不能直接delete：

OGRFeature* feature = vectorLayer->GetNextFeature();
OGRGeometry* clipGeometry = feature->GetGeometryRef();
// 后续使用完毕后销毁特征对象
OGRFeature::DestroyFeature(feature);

3. 冗余投影变换 #

已知CRS一致，无需创建OGRCoordinateTransformation，直接移除相关代码。

4. TransformCutlineToSource调用时机错误 #

该函数需在WarpOperation初始化前调用，确保裁剪几何与源栅格坐标匹配：

// 转换裁剪几何到源坐标系
TransformCutlineToSource(rasterDataset, reinterpret_cast<OGRGeometryH>(cutlineGeometry), &(warpOptions->papszWarpOptions), nullptr);

// 再初始化WarpOperation
GDALWarpOperation warpOperation;
CPLErr err = warpOperation.Initialize(warpOptions);

5. 输出数据集配置优化 #

创建GTiff时可添加压缩选项减少文件体积，同时设置NoData值标记裁剪外区域：

char** papszOptions = CSLSetNameValue(nullptr, "COMPRESS", "LZW");
GDALDataset* outputDataset = driver->Create(outputRasterPath, rasterDataset->GetRasterXSize(), rasterDataset->GetRasterYSize(), 1, GDT_Float32, papszOptions);
CSLDestroy(papszOptions);

// 设置NoData值
outputDataset->GetRasterBand(1)->SetNoDataValue(-9999.0);

6. 移除冗余投影转换器 #

因CRS一致，无需创建GDALGenImgProjTransformer，直接删除相关代码。

修正后的完整代码 #

#include "gdal_priv.h"
#include "cpl_conv.h"
#include "ogr_spatialref.h"
#include "ogr_geometry.h"
#include "gdalwarper.h"

int main()
{
    GDALAllRegister();
    const char* rasterPath = "input.tif";
    const char* vectorPath = "polygon.shp";
    const char* outputRasterPath = "output.tif";

    // 打开栅格数据集（只读）
    GDALDataset* rasterDataset = static_cast<GDALDataset*>(GDALOpen(rasterPath, GA_ReadOnly));
    if (!rasterDataset) {
        printf("无法打开栅格文件\n");
        return 1;
    }

    // 打开矢量数据集（只读）
    GDALDataset* vectorDataset = static_cast<GDALDataset*>(GDALOpenEx(vectorPath, GDAL_OF_READONLY, nullptr, nullptr, nullptr));
    if (!vectorDataset) {
        printf("无法打开矢量文件\n");
        GDALClose(rasterDataset);
        return 1;
    }

    OGRLayer* vectorLayer = vectorDataset->GetLayer(0);
    if (!vectorLayer) {
        printf("无法获取矢量图层\n");
        GDALClose(vectorDataset);
        GDALClose(rasterDataset);
        return 1;
    }

    // 获取裁剪几何并管理特征对象内存
    OGRFeature* feature = vectorLayer->GetNextFeature();
    if (!feature) {
        printf("无法获取矢量特征\n");
        GDALClose(vectorDataset);
        GDALClose(rasterDataset);
        return 1;
    }
    OGRGeometry* clipGeometry = feature->GetGeometryRef();
    if (!clipGeometry) {
        printf("无法获取几何对象\n");
        OGRFeature::DestroyFeature(feature);
        GDALClose(vectorDataset);
        GDALClose(rasterDataset);
        return 1;
    }
    OGRGeometry* cutlineGeometry = clipGeometry->clone();
    OGRFeature::DestroyFeature(feature);

    // 创建带压缩选项的输出数据集
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    char** papszOptions = CSLSetNameValue(nullptr, "COMPRESS", "LZW");
    GDALDataset* outputDataset = driver->Create(outputRasterPath, rasterDataset->GetRasterXSize(), rasterDataset->GetRasterYSize(), 1, GDT_Float32, papszOptions);
    CSLDestroy(papszOptions);
    if (!outputDataset) {
        printf("无法创建输出数据集\n");
        delete cutlineGeometry;
        GDALClose(vectorDataset);
        GDALClose(rasterDataset);
        return 1;
    }

    // 设置输出投影、地理变换及NoData值
    outputDataset->SetProjection(rasterDataset->GetProjectionRef());
    double adGeoTransform[6];
    rasterDataset->GetGeoTransform(adGeoTransform);
    outputDataset->SetGeoTransform(adGeoTransform);
    outputDataset->GetRasterBand(1)->SetNoDataValue(-9999.0);

    // 初始化Warp选项
    GDALWarpOptions* warpOptions = GDALCreateWarpOptions();
    warpOptions->hSrcDS = rasterDataset;
    warpOptions->hDstDS = outputDataset;
    warpOptions->nBandCount = 1;
    warpOptions->panSrcBands = (int*)CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    warpOptions->panSrcBands[0] = 1;
    warpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    warpOptions->panDstBands[0] = 1;
    warpOptions->papszWarpOptions = CSLDuplicate(nullptr);
    warpOptions->eWorkingDataType = GDT_Float32;
    warpOptions->pfnProgress = GDALTermProgress;
    warpOptions->hCutline = cutlineGeometry;

    // 转换裁剪几何到源坐标系
    TransformCutlineToSource(rasterDataset, reinterpret_cast<OGRGeometryH>(cutlineGeometry), &(warpOptions->papszWarpOptions), nullptr);

    // 初始化并执行Warp操作
    GDALWarpOperation warpOperation;
    CPLErr err = warpOperation.Initialize(warpOptions);
    if (err != CE_None) {
        printf("Warp初始化失败\n");
        GDALDestroyWarpOptions(warpOptions);
        delete cutlineGeometry;
        GDALClose(outputDataset);
        GDALClose(vectorDataset);
        GDALClose(rasterDataset);
        return 1;
    }

    err = warpOperation.ChunkAndWarpImage(0, 0, outputDataset->GetRasterXSize(), outputDataset->GetRasterYSize());
    if (err != CE_None) {
        printf("Warp执行失败\n");
    }

    // 清理资源
    GDALDestroyWarpOptions(warpOptions);
    delete cutlineGeometry;
    GDALClose(outputDataset);
    GDALClose(vectorDataset);
    GDALClose(rasterDataset);

    return 0;
}

实现-crop_to_cutline效果（可选） #

若需要输出栅格范围匹配裁剪多边形，需先计算多边形外接矩形，再调整输出栅格尺寸和地理变换：

OGREnvelope envelope;
cutlineGeometry->GetEnvelope(&envelope);
// 根据源分辨率计算输出尺寸
int outputXSize = static_cast<int>((envelope.MaxX - envelope.MinX) / adGeoTransform[1]) + 1;
int outputYSize = static_cast<int>((envelope.MaxY - envelope.MinY) / (-adGeoTransform[5])) + 1;
// 构建新的地理变换
double outputGeoTransform[6] = {envelope.MinX, adGeoTransform[1], 0, envelope.MaxY, 0, adGeoTransform[5]};
// 使用新尺寸和地理变换创建输出数据集
GDALDataset* outputDataset = driver->Create(outputRasterPath, outputXSize, outputYSize, 1, GDT_Float32, papszOptions);
outputDataset->SetGeoTransform(outputGeoTransform);

内容的提问来源于stack exchange，提问作者Alireza Rahimi