C++ GDAL 计算栅格影像覆盖范围 (排除 0 值区域)

使用 C++ 和 GDAL 计算栅格影像覆盖范围 (排除 0 值区域)

本文将介绍如何使用 C++ 语言和 GDAL 库来计算栅格影像的覆盖范围，并排除 0 值区域。最终结果将以 OGRPolygon 格式返回。

代码示例

#include <gdal_priv.h>
#include <ogrsf_frmts.h>

int main(int argc, char* argv[])
{
    GDALAllRegister(); // 注册所有的驱动

    // 打开栅格影像
    GDALDataset* poDataset = (GDALDataset*)GDALOpen(argv[1], GA_ReadOnly);
    if (poDataset == NULL)
    {
        printf("无法打开栅格影像文件！\n");
        exit(1);
    }

    // 获取影像范围
    double adfGeoTransform[6];
    poDataset->GetGeoTransform(adfGeoTransform);
    double xmin = adfGeoTransform[0];
    double ymax = adfGeoTransform[3];
    double xmax = xmin + adfGeoTransform[1] * poDataset->GetRasterXSize();
    double ymin = ymax + adfGeoTransform[5] * poDataset->GetRasterYSize();

    // 创建 OGC 多边形对象
    OGRPolygon polygon;
    OGRLinearRing ring;
    ring.addPoint(xmin, ymin);
    ring.addPoint(xmin, ymax);
    ring.addPoint(xmax, ymax);
    ring.addPoint(xmax, ymin);
    ring.addPoint(xmin, ymin);
    polygon.addRing(&ring);

    // 获取栅格影像的 NoData 值
    double nodata = poDataset->GetRasterBand(1)->GetNoDataValue();

    // 遍历栅格影像像素，计算覆盖范围
    for (int j = 0; j < poDataset->GetRasterYSize(); j++)
    {
        for (int i = 0; i < poDataset->GetRasterXSize(); i++)
        {
            double value;
            poDataset->GetRasterBand(1)->RasterIO(GF_Read, i, j, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);
            if (value != nodata && value != 0.0)
            {
                double x = xmin + adfGeoTransform[1] * i;
                double y = ymax - adfGeoTransform[5] * j;
                OGRPoint point(x, y);
                polygon.intersection(&point);
            }
        }
    }

    // 输出结果
    OGRGeometry* geometry = &polygon;
    OGRSpatialReference* srs = poDataset->GetSpatialRef();
    geometry->assignSpatialReference(srs);
    char* wkt;
    geometry->exportToWkt(&wkt);
    printf('覆盖范围为：%s\n', wkt);

    // 释放资源
    CPLFree(wkt);
    GDALClose(poDataset);

    return 0;
}

代码说明

注册驱动: 使用 GDALAllRegister() 函数注册所有 GDAL 驱动。
打开栅格影像: 使用 GDALOpen() 函数打开栅格影像文件。
获取影像范围: 使用 GetGeoTransform() 函数获取影像的地理变换矩阵，并计算出影像的最小坐标 (xmin, ymin) 和最大坐标 (xmax, ymax)。
创建 OGRPolygon 对象: 创建一个 OGRPolygon 对象，并添加影像四个角点作为多边形边界。
获取 NoData 值: 使用 GetNoDataValue() 函数获取栅格影像的 NoData 值。
遍历像素: 遍历栅格影像所有像素，如果像素值不为 NoData 值且不为 0，则将其坐标转换为地理坐标，并创建 OGRPoint 对象。
计算交集: 使用 intersection() 函数计算 OGRPoint 对象与 OGRPolygon 对象的交集，从而更新 OGRPolygon 对象的边界，最终得到影像的覆盖范围。
输出结果: 将 OGRPolygon 对象转换为 WKT 格式的字符串并输出。

注意事项

该代码只考虑了单波段栅格影像。对于多波段影像，需要在遍历像素时同时遍历所有波段，以保证计算出的覆盖范围包含所有波段的有效像素。
如果栅格影像存在其他需要排除的无效值，可以在遍历像素时添加额外的判断条件。

应用场景

识别栅格影像中有效数据区域。
在进行地理空间分析时，可以使用该方法获得影像的有效范围，避免对无效数据进行处理。
可以用于对影像进行裁剪，只保留有效数据区域。

总结

本文介绍了如何使用 C++ 和 GDAL 库来计算栅格影像的覆盖范围，并排除 0 值区域。该方法可以帮助您快速识别影像中的有效数据区域，并为后续的地理空间分析提供基础数据。