GISAlgorithm学习笔记第2篇：GDAL影像裁剪

星星也在雾里

已于 2024-07-09 16:14:02 修改

阅读量958

点赞数 17

分类专栏： GISAlgorithm学习笔记文章标签： visual studio arcgis

于 2024-07-06 00:54:50 首次发布

本文链接：https://blog.csdn.net/weixin_62934777/article/details/140213443

版权

GISAlgorithm学习笔记专栏收录该内容

3 篇文章 0 订阅

订阅专栏

内容摘要：

本文介绍利用GDAL开源库对tif格式的影像进行裁剪的操作，并介绍其原理涉及的多个重要函数。

往期回顾：

GISAlgorithm学习笔记第1篇：VS2017配置GDAL库

一、开发环境

二、实现原理

三、完整代码

一、开发环境

1.语言：C++

2.编程环境：Visual Studio2017

3.开源库：GDAL

二、实现原理

1.打开源图像文件：

GDALOpen函数是GDAL库中用于打开地理空间数据文件的关键函数。它可以处理各种栅格数据格式，包括常见的GeoTIFF、JPEG、PNG、HDF、NetCDF等。通过调用GDALOpen，可以获取一个表示数据集的指针，该指针可以用于后续的读写操作。

函数定义：

GDALDatasetH GDALOpen(const char * pszFilename, GDALAccess eAccess);

参数说明：

pszFilename：要打开的文件的路径（字符串）。
eAccess：文件的访问模式，类型为GDALAccess。可以是以下两种值之一：
- GA_ReadOnly：只读模式打开文件。
- GA_Update：读写模式打开文件。

本文示例：

使用GDALOpen函数打开源图像文件，以只读模式（GA_ReadOnly）获取一个指向GDALDataset对象的指针。

GDALDataset *poSrcDS = (GDALDataset *)GDALOpen(pSrcFilename, GA_ReadOnly);
if (poSrcDS == NULL) {
    cout << "Cannot open file " << pSrcFilename << endl;
    return 1;
}

2. 定义裁剪区域：

裁剪区域的起始坐标（xOff和yOff）和大小（dstWidth和dstHeight）可以根据需要进行设置。

本文示例：

  // 定义裁剪区域的起始点（左上角）和宽度、高度
    int xOff = 100; // 裁剪区域的起始x坐标
    int yOff = 100; // 裁剪区域的起始y坐标
    int dstWidth = 200; // 裁剪区域的宽度
    int dstHeight = 200; // 裁剪区域的高度

3. 读取裁剪区域的数据：

RasterIO函数是 GDAL 中用于读写栅格数据的核心函数之一，提供了灵活且高效的方法来读取和写入数据集的栅格数据。RasterIO函数支持任意大小的读写操作，可以对整个栅格或栅格的任意部分进行操作，同时支持多种数据类型，并且在读写时可以进行数据类型转换。

函数定义：

CPLErr GDALDataset::RasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
                             void * pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
                             int nBandCount, int * panBandMap, 
                             GSpacing nPixelSpace, GSpacing nLineSpace, GSpacing nBandSpace)

参数说明：

eRWFlag：指定是读操作 (GF_Read) 还是写操作 (GF_Write)。
nXOff, nYOff：要读取或写入的图像区域的左上角像素的 X 和 Y 偏移量。
nXSize, nYSize：要读取或写入的图像区域的宽度和高度。
pData：指向用于保存读取数据或包含要写入数据的缓冲区的指针。
nBufXSize, nBufYSize：缓冲区的宽度和高度。如果与 nXSize 和 nYSize不同，则数据将被重采样。
eBufType: 缓冲区中数据的类型（如 GDT_Byte，GDT_UInt16，GDT_Float32等）。
nBandCount：要处理的波段数。
panBandMap：波段索引数组（波段从 1 开始计数）。
nPixelSpace, nLineSpace, nBandSpace：缓冲区内数据的步长，分别表示像素间距、行间距和波段间距。

本文示例：

 // 从源波段中读取裁剪区域的数据到分配的内存中
    poBand->RasterIO(GF_Read, xOff, yOff, dstWidth, dstHeight, pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

4. 创建目标图像文件：

GDALDriver类和相关函数是 GDAL 中用于管理不同栅格数据格式的关键部分。GDALDriver提供了创建、删除和识别数据集的功能。每个驱动程序负责一个或多个特定的文件格式，例如GeoTIFF、JPEG、PNG 等。

相关函数如下。

（1）获取一个特定格式的驱动程序，使用GetGDALDriverManager获取驱动程序管理器，然后调用 GetDriverByName函数：

GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

（2）创建数据集时可以使用GDALDriver类的Create函数创建一个新的数据集：

函数定义：

GDALDataset *GDALDriver::Create(const char *pszFilename, int nXSize, int nYSize, int nBands, GDALDataType eType, char **papszOptions);

参数说明：

pszFilename：要创建的数据集的文件名。
nXSize：数据集的宽度（像素）。
nYSize：数据集的高度（像素）。
nBands：数据集中的波段数。
eType：波段的数据类型（如 GDT_Byte, GDT_UInt16, GDT_Float32 等）。
papszOptions：创建选项（可选），如压缩类型等。

本文示例：

GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == NULL) {
    CPLFree(pafScanline);
    GDALClose(poSrcDS);
    return 1;
}

GDALDataset *poDstDS = poDriver->Create(pDstFilename, dstWidth, dstHeight, 1, GDT_Float32, NULL);
if (poDstDS == NULL) {
    cout << "Could not create dataset" << endl;
    CPLFree(pafScanline);
    GDALClose(poSrcDS);
    return 1;
}

5. 设置地理信息：

GetGeoTransform函数用于获取源数据集的地理变换参数。

SetGeoTransform函数用于设置栅格数据集的地理变换参数。

函数定义：

CPLErr GDALDataset::GetGeoTransform(double * padfTransform);//获取
CPLErr GDALDataset::SetGeoTransform(double *padfTransform);//设置

参数说明：

padfTransform: 一个包含 6 个双精度值的数组，用于存储地理变换参数。这些参数定义了图像像素坐标与地理坐标之间的仿射变换关系。
```
adfGeoTransform[0] += xOff * adfGeoTransform[1] + yOff * adfGeoTransform[2];
adfGeoTransform[3] += xOff * adfGeoTransform[4] + yOff * adfGeoTransform[5];
```

GetProjectionRef函数用于获取栅格数据集投影信息。

SetProjectionRef函数用于获取栅格数据集投影信息。

函数定义：

const char* GDALDataset::GetProjectionRef();//获取
CPLErr GDALDataset::SetProjection(const char *pszProjection);//设置

参数说明：

pszProjection: 一个指向表示投影信息的 WKT 格式字符串的指针。
WKT（Well-Known Text）是一种用于表示地理信息系统中的几何对象和地理坐标系统的文本格式，由 Open Geospatial Consortium (OGC) 标准化，广泛应用于各种 GIS 软件和数据库中。

本文示例：

double adfGeoTransform[6];
if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None) {
    adfGeoTransform[0] += xOff * adfGeoTransform[1] + yOff * adfGeoTransform[2];
    adfGeoTransform[3] += xOff * adfGeoTransform[4] + yOff * adfGeoTransform[5];
    poDstDS->SetGeoTransform(adfGeoTransform);
}

const char* p_WKT = poSrcDS->GetProjectionRef();
if (p_WKT != NULL) {
    poDstDS->SetProjection(p_WKT);
}

6. 写入裁剪后的数据：

使用RasterIO函数将读取到的裁剪区域数据写入目标图像文件的对应位置。

GDALRasterBand *poDstBand = poDstDS->GetRasterBand(1);
poDstBand->RasterIO(GF_Write, 0, 0, dstWidth, dstHeight, (void*)pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

三、完整代码

裁剪结果在ArcMap或ENVI软件中显示更加明显。

注意：文件路径必须为\\双斜杠！！！

#include <gdal_priv.h>
#include <cpl_conv.h>
#include <iostream>
#include <gdal.h>

using namespace std;

int main() {
    GDALAllRegister(); // 注册所有GDAL驱动

    const char* pSrcFilename = "D:\\xx.tif"; // 修改为你的源文件路径
    const char* pDstFilename = "D:\\xxafter.tif"; // 修改为你的目标文件路径,保证修改前和修改后文件名不一样

    // 打开源栅格文件
    GDALDataset *poSrcDS = (GDALDataset *)GDALOpen(pSrcFilename, GA_ReadOnly);
    if (poSrcDS == NULL) { // 若打开文件失败，输出错误信息
        cout << "Cannot open file " << pSrcFilename << endl;
        return 1;
    }

    // 获取源文件中的第一个波段（波段从1开始计数）
    GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);

    // 获取源文件的宽度和高度
    int srcWidth = poBand->GetXSize();
    int srcHeight = poBand->GetYSize();

    // 定义裁剪区域的起始点（左上角）和宽度、高度
    int xOff = 100; // 裁剪区域的起始x坐标
    int yOff = 100; // 裁剪区域的起始y坐标
    int dstWidth = 200; // 裁剪区域的宽度
    int dstHeight = 200; // 裁剪区域的高度

    // 为扫描线分配内存
    float *pafScanline = (float *)CPLMalloc(sizeof(float) * dstWidth * dstHeight);

    // 从源波段中读取裁剪区域的数据到分配的内存中
    poBand->RasterIO(GF_Read, xOff, yOff, dstWidth, dstHeight, pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

    // 获取GTiff驱动，用于创建目标文件
    GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

    // 若获取失败，退出程序
    if (poDriver == NULL) {
        CPLFree(pafScanline);
        GDALClose(poSrcDS);
        return 1;
    }

    // 使用GTiff驱动创建一个新的栅格文件
    GDALDataset *poDstDS = poDriver->Create(pDstFilename, dstWidth, dstHeight, 1, GDT_Float32, NULL);
    if (poDstDS == NULL) {
        cout << "Could not create dataset" << endl;
        CPLFree(pafScanline);
        GDALClose(poSrcDS);
        return 1;
    }

    // 获取源文件的地理变换参数, 投影信息
    double adfGeoTransform[6];
    if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None) {
        // 设置新的仿射变换矩阵
        adfGeoTransform[0] += xOff * adfGeoTransform[1] + yOff * adfGeoTransform[2];
        adfGeoTransform[3] += xOff * adfGeoTransform[4] + yOff * adfGeoTransform[5];
        poDstDS->SetGeoTransform(adfGeoTransform);
    }

    const char* p_WKT = poSrcDS->GetProjectionRef();
    if (p_WKT != NULL) {
        // 设置目标文件的投影信息
        poDstDS->SetProjection(p_WKT);
    }

    // 获取目标文件中的第一个波段
    GDALRasterBand *poDstBand = poDstDS->GetRasterBand(1);

    // 将读取的数据写入目标波段
    poDstBand->RasterIO(GF_Write, 0, 0, dstWidth, dstHeight, (void*)pafScanline, dstWidth, dstHeight, GDT_Float32, 0, 0);

    // 释放扫描线内存
    CPLFree(pafScanline);

    // 关闭目标文件和源文件
    GDALClose((GDALDatasetH)poDstDS);
    GDALClose((GDALDatasetH)poSrcDS);

    return 0;
}

欢迎交流🌹🌹