GDAL影像读写操作

读取影像并获取相关影像参数

#include "gdal_priv.h"
#include<iostream>
#include "cpl_conv.h" // 用于 CPLGetLastErrorMsg() 等函数
#include "cpl_string.h"

int main()
{

    const char* pszFilename = "E:\\gdal\\osgeopy-data\\osgeopy-data-landsat-washington\\osgeopy-data\\Landsat\\Washington\\p047r027_7t20000730_z10_nn10.tif";
    // GDALDataset唯一类型数据指针
    // GDAL数据集的唯一指针类型。适用于以非共享模式打开的数据集，并且其引用计数器未被手动修改的情况。
    GDALDatasetUniquePtr poDataset;
    // GDALDataset poDataset; // GDALDataset数据类型
    // Register all known configured GDAL drivers
    GDALAllRegister();
    const GDALAccess eAccess = GA_ReadOnly;
    // GDALDatasetH C语言绑定类型：“Opaque type”用于C++ GDALDataset类的C语言绑定。
    // Convert a GDALDatasetH to a GDALDataset*
    // GDALDataset::FromHandle
    // GDALDataset 的唯一指针类型。
    // GDALDatasetUniquePtr// 智能指针类型
    // 适用于以非共享模式打开的数据集，且其引用计数器未被手动修改的情况。
    // Once the drivers are registered, the application should call the free standing GDALOpen() 
    //function to open a dataset, passing the name of the dataset and the access desired (GA_ReadOnly:只读 or GA_Update：读写).
    poDataset = GDALDatasetUniquePtr(GDALDataset::FromHandle(GDALOpen(pszFilename, eAccess)));
    if (!poDataset)
    {
        fprintf(stderr, "GDALOpen failed - %s\n", CPLGetLastErrorMsg()); // handle error
    }
    // 获取变换参数
    double  adfGeoTransform[6];
    //数据驱动类型
    printf("Driver: %s/%s\n",
        poDataset->GetDriver()->GetDescription(),
        poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
    // 获取栅格信息
    printf("Size is %dx%dx%d\n",
        poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
        poDataset->GetRasterCount());
    // 获取投影信息
    if (poDataset->GetProjectionRef() != NULL)
        printf("Projection is `%s'\n", poDataset->GetProjectionRef());
    // 获取地理变换信息：adfGeoTransform；如果数据获取失败，返回一个CPLErr错误类型
    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None)
    {
        // 左上角坐标点
        printf("Origin = (%.6f,%.6f)\n",
            adfGeoTransform[0], adfGeoTransform[3]);
        // 图像大小
        printf("Pixel Size = (%.6f,%.6f)\n",
            adfGeoTransform[1], adfGeoTransform[5]);
    }

    // 栅格波段类型 GDALRasterBand：获取影像波段
    GDALRasterBand* poBand;
    int             nBlockXSize, nBlockYSize;
    int             bGotMin, bGotMax;
    double          adfMinMax[2];
    // 获取影像波段，及波段相应的属性
    poBand = poDataset->GetRasterBand(1);
    poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
    printf("Block=%dx%d Type=%s, ColorInterp=%s\n",
        nBlockXSize, nBlockYSize,
        GDALGetDataTypeName(poBand->GetRasterDataType()),
        GDALGetColorInterpretationName(
            poBand->GetColorInterpretation()));
    adfMinMax[0] = poBand->GetMinimum(&bGotMin);
    adfMinMax[1] = poBand->GetMaximum(&bGotMax);
    if (!(bGotMin && bGotMax))
        // 强制转换为C绑定类型：Opaque type used for the C bindings of the C++ GDALRasterBand class
        // 返回图像的最大值与最小值
        GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
    printf("Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1]);
    // 可获取的波段数
    if (poBand->GetOverviewCount() > 0)
        // Return the number of overview layers available.
        printf("Band has %d overviews.\n", poBand->GetOverviewCount());
    // Fetch the color table associated with band.
    if (poBand->GetColorTable() != NULL)
        printf("Band has a color table with %d entries.\n",
            poBand->GetColorTable()->GetColorEntryCount());
    // 读取栅格影像像元信息
    float* pafScanline;
    int   nXSize = poBand->GetXSize();
    int   nYSize = poBand->GetYSize();
    // 申请内存空间,读取栅格数据
    pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize* nYSize);
    poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize,
        pafScanline, nXSize, nYSize, GDT_Float32,
        0, 0);
    // int sizes = pafScanline.size();
    // GDALDataset*数据集需要显示调用GDALClose(hds);
    // GDALDatasetUniquePtr指向的数据无需显式调用 GDALClose()，因为 GDALDatasetUniquePtr 会在离开作用域时自动处理


    const char* pszFormat = "GTiff";
    GDALDriver* poDriver;
    char** papszMetadata;
    poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (poDriver == NULL)
        exit(1);
    papszMetadata = poDriver->GetMetadata();
    if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
        printf("Driver %s supports Create() method.\n", pszFormat);
    if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
        printf("Driver %s supports CreateCopy() method.\n", pszFormat);
    /* 通过复制创建新的文件 */ 
    const char* newFilename = "E:\\gdal\\osgeopy-data\\osgeopy-data-landsat-washington\\osgeopy-data\\Landsat\\Washington\\new.tif";
    // GDALDataset类型
    // 若读取文件时，读取了一个空文件，则poSrcDS未被正确初始化
    GDALDataset* poSrcDS =(GDALDataset*)GDALOpen(pszFilename, GA_ReadOnly);
    GDALDataset* poDstDS;
    // 将数据复制到新的驱动文件中
    /*poDstDS = poDriver->CreateCopy(newFilename, poSrcDS, FALSE,
        NULL, NULL, NULL);*/
    /* Once we're done, close properly the dataset */
    // 写入数据
    /*if (poDstDS != NULL)
        GDALClose((GDALDatasetH)poDstDS);
    GDALClose((GDALDatasetH)poSrcDS);*/
    /*设置元数据，传递给驱动器/新的数据集*/
    // 申请元数据指针
    char** papszOptions = NULL;
    papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES");
    papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "PACKBITS");
    poDstDS = poDriver->CreateCopy(newFilename, poSrcDS, FALSE,
        papszOptions, GDALTermProgress, NULL);
    /* Once we're done, close properly the dataset */
    if (poDstDS != NULL)
        GDALClose((GDALDatasetH)poDstDS);
    // 释放指针
    CSLDestroy(papszOptions);

    
    /* GDALCreate 方式创建新的数据集*/
    // 需要将参数逐个写入
    GDALDataset* hDstDS;// 数据集
    char** papszOptions_new = NULL;
    // 获取驱动句柄  
    GDALAllRegister();
    // 创建文件名称
    const char* pszDstFilename = "E:\\gdal\\osgeopy-data\\osgeopy-data-landsat-washington\\osgeopy-data\\Landsat\\Washington\\a.tif ";
    // 驱动器
    GDALDriver* hDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); // 假设我们想要创建一个GeoTIFF文件
    // 申请元数据指针
    char** papszOption = NULL;
    papszOption = CSLSetNameValue(papszOption, "TILED", "YES");
    papszOption = CSLSetNameValue(papszOption, "COMPRESS", "PACKBITS");

    // 构造数据集内存冲突，即指针未被正确初始化或已被释放；
    hDstDS = poDriver->Create(pszDstFilename, 512, 512, 1, GDT_Byte,
        papszOption);
    // 地理变换参数
    double adfGeoTransform_new[6] = { 444720, 30, 0, 3751320, 0, -30 };
    // 空间参考
    OGRSpatialReference oSRS;
    char* pszSRS_WKT = NULL;
    // 波段数据类型
    GDALRasterBand* poBand_new = NULL;
    GByte abyRaster[512 * 512];
    // GDALDatasetH指向一个指针类型
    // GDALDataset* 指针类型；如果使用GDALDatasetH数据类型来指向一个GDALDataset*的类函数
    // 会提示：表达式必须包含指向类的指针类型
    hDstDS->SetGeoTransform(adfGeoTransform_new);
    // 设置投影与坐标参考
    oSRS.SetUTM(11, TRUE);
    oSRS.SetWellKnownGeogCS("NAD27");
    // 将坐标系与投影导出为WKT格式
    oSRS.exportToWkt(&pszSRS_WKT);
    // 设置投影与坐标参考系
    hDstDS->SetProjection(pszSRS_WKT);
    // CPLFree 是GDAL/OGR库中用于释放内存的函数;
    CPLFree(pszSRS_WKT);
    poBand_new = hDstDS->GetRasterBand(1);
    poBand_new->RasterIO(GF_Write, 0, 0, 512, 512,
        abyRaster, 512, 512, GDT_Byte, 0, 0);
    /* Once we're done, close properly the dataset */
    GDALClose((GDALDatasetH)hDstDS);
    return 0;
}
/*
Driver: GTiff/GeoTIFF
Size is 8849x8023x1
Projection is `PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'
Origin = (343724.250000,5369585.250000)
Pixel Size = (28.500000,-28.500000)
Block=8849x1 Type=Byte, ColorInterp=Gray
Min=0.000d, Max=255.000
Band has 5 overviews.
*/
/*
Driver: GTiff/GeoTIFF
Size is 8849x8023x1
Projection is `PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]'
Origin = (343724.250000,5369585.250000)
Pixel Size = (28.500000,-28.500000)
Block=8849x1 Type=Byte, ColorInterp=Gray
Min=0.000d, Max=255.000
Band has 5 overviews.
*/

// GDALDatasetH指向一个指针类型, GDALDataset* 指针类型；
如果使用GDALDatasetH数据类型来指向一个GDALDataset*的类函数会提示：表达式必须包含指向类的指针类型
可以适用强制类型转换，对两种类型进行转换

// c++
GDALDataset *poDstDS;
char **papszOptions = NULL;
poDstDS = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte,
                            papszOptions );
// c
GDALDatasetH hDstDS;
char **papszOptions = NULL;
hDstDS = GDALCreate( hDriver, pszDstFilename, 512, 512, 1, GDT_Byte,
                    papszOptions );