物流网站建设推广,摄影网站开发综述,烟台城乡建设学校官方网站,天津滨海新区网站建设本文介绍基于C 语言的GDAL模块#xff0c;按照给定的像元行数与列数#xff0c;批量裁剪大量多波段栅格遥感影像文件#xff0c;并将所得到的裁剪后新的多波段遥感影像文件保存在指定路径中的方法。 在之前的文章中#xff0c;我们多次介绍了在不同平台#xff0c;或基于不… 本文介绍基于C 语言的GDAL模块按照给定的像元行数与列数批量裁剪大量多波段栅格遥感影像文件并将所得到的裁剪后新的多波段遥感影像文件保存在指定路径中的方法。 在之前的文章中我们多次介绍了在不同平台或基于不同代码语言对栅格遥感影像加以裁剪、批量裁剪的方法主要包括Python中ArcPy基于矢量范围批量裁剪大量栅格遥感影像https://blog.csdn.net/zhebushibiaoshifu/article/details/128307676、Google Earth Engine谷歌地球引擎GEE矢量数据裁剪栅格数据https://blog.csdn.net/zhebushibiaoshifu/article/details/117390431、基于ENVI实现栅格遥感影像按图层行列号与像元数量划定矩形研究区域并裁剪https://blog.csdn.net/zhebushibiaoshifu/article/details/118978851等而本文我们就介绍一下基于**C**语言的GDAL模块实现批量裁剪需求的方法。 首先来看一下本文的需求。现在有一个文件夹如下图所示其中具有多个.tiff格式的多波段遥感影像文件为了方便我们这里文件夹内就只有2个文件但实际上一般我们批量处理的需求肯定远远大于这个数量。 我们希望实现的就是基于这个文件夹内每一景遥感影像将其左上角100 * 100像元的这一部分给裁剪下来如下图所示并分别保存为新的遥感影像文件其中新的文件名称就在原有文件名称后加一个_C后缀即可并存放在另一个指定的结果文件夹中。我们希望裁剪后的遥感影像和原有的遥感影像对比起来呈现如下图所示的情况。 本文所用代码如下。
#include iostream
#include string
#include gdal/gdal.h
#include gdal/gdal_priv.h
using namespace std;int main()
{GDALAllRegister();string inputFolder /home/cppGDAL/TIF/WFV1_2021_STB/test;string outputFolder /home/cppGDAL/TIF/WFV1_2021_STB_C;CPLStringList fileList;fileList VSIReadDir(inputFolder.c_str());for (int i 0; i fileList.size(); i){string inputImagePath fileList[i];if (!EQUAL(CPLGetExtension(inputImagePath.c_str()), tiff)){continue;}string full_path inputFolder / inputImagePath;GDALDataset *poDataset (GDALDataset *)GDALOpen(full_path.c_str(), GA_ReadOnly);int width poDataset-GetRasterXSize();int height poDataset-GetRasterYSize();int xOffset 0;int yOffset 0;int xSize 100;int ySize 100;cout full_path;string outputImageName (outputFolder / inputImagePath.substr(0, inputImagePath.size() - 5) _C.tiff);cout outputImageName;GDALDriver *poDriver GetGDALDriverManager()-GetDriverByName(GTiff);GDALDataset *poOutputDataset poDriver-Create(outputImageName.c_str(), xSize, ySize, 4, GDT_Float32, NULL);poOutputDataset-SetProjection(poDataset-GetProjectionRef());double adfGeoTransform[6];poDataset-GetGeoTransform(adfGeoTransform);// adfGeoTransform[1] adfGeoTransform[1] / (width / xSize);// adfGeoTransform[5] adfGeoTransform[5] / (height / ySize);poOutputDataset-SetGeoTransform(adfGeoTransform);for (int bandIndex 1; bandIndex 4; bandIndex){float *buffer new float[xSize * ySize];GDALRasterBand *poBand poDataset-GetRasterBand(bandIndex);GDALRasterBand *poOutputBand poOutputDataset-GetRasterBand(bandIndex);poBand-RasterIO(GF_Read, xOffset, yOffset, xSize, ySize, buffer, xSize, ySize, GDT_Float32, 0, 0, NULL);poOutputBand-RasterIO(GF_Write, 0, 0, xSize, ySize, buffer, xSize, ySize, GDT_Float32, 0, 0, NULL);delete[] buffer;}GDALClose(poOutputDataset);GDALClose(poDataset);}GDALDestroyDriverManager();return 0;
}我们来介绍一下上述代码的具体内容。 首先我们需要通过GDALAllRegister();来注册所有的GDAL驱动器。同时我们定义了输入和输出文件夹路径——inputFolder就是存储输入遥感影像待裁剪的遥感影像的文件夹路径outputFolder则是存储结果遥感影像的文件夹路径。 其次我们通过CPLStringList fileList;定义一个字符串列表用于存储文件夹中的文件列表并使用VSIReadDir函数读取输入文件夹中的所有文件并将结果存储在fileList中。 接下来我们使用循环迭代处理每个文件。首先通过string inputImagePath fileList[i];获取当前文件的路径如果文件的扩展名不是tiff则跳过该文件。接下来对于文件的扩展名是tiff的我们构建完整的输入文件路径并使用GDALOpen函数打开输入文件返回一个GDALDataset对象存储在poDataset中。 接下来我们即可获取输入文件的宽度和高度并定义裁剪区域的偏移量左上角像元的位置、宽度和高度。前面提到了我这里就是需要在原本遥感影像的最左上角所以偏移量均为0裁剪下来100 * 100像元的这一部分。其次构建输出文件的路径并使用GetGDALDriverManager()-GetDriverByName函数获取GTiff驱动器对象存储在poDriver中。随后我们使用poDriver-Create函数创建输出文件返回一个GDALDataset对象存储在poOutputDataset中。 接下来这个部分需要稍微注意一下。首先我们使用poOutputDataset-SetProjection设置输出文件的投影信息即与输入文件相同的投影其次使用poDataset-GetGeoTransform获取输入文件的地理变换参数存储在adfGeoTransform数组中。由于在我这里裁剪后遥感影像的像元大小即单个像元的长度与宽度没有改变且裁剪前后栅格遥感影像的左上角像元没有发生变化所以新的栅格遥感影像的地理变换参数和老的栅格遥感影像比起来无需有任何改变但是如果大家的裁剪需求不是这样的话比如像元大小发生变化了或者是裁剪并不是从左上角像元开始的那么就需要调整这6个地理变换参数——至于怎么变这就比较复杂了我也还没完全搞清楚大家就结合自己的实际需求到GDAL官网查阅即可。最后我们使用poOutputDataset-SetGeoTransform设置输出文件的地理变换参数在我这里就是与输入文件完全相同的地理变换参数。 代码最后我们使用循环迭代处理每个波段我这里每一个遥感影像都是共4个波段。首先创建一个大小为xSize * ySize的浮点型缓冲区并使用poBand-RasterIO从输入文件中读取对应波段的像元数据到缓冲区接下来使用poOutputBand-RasterIO将缓冲区中的数据写入到输出文件对应波段中。随后即可释放缓冲区内存并关闭输出文件和输入文件。 运行上述代码我们即可在结果文件夹中看到已经裁剪好的遥感影像文件且新的文件的文件名称也符合我们的要求如下图所示。 至此大功告成。
欢迎关注疯狂学习GIS
