This article introduces the GDAL module based on the C + + language. According to the given number of pixel rows and number of columns strong>, batch crop a large number of multi-band raster remote sensing image files, and save the resulting cropped new multi-band remote sensing image files Method at the specified path.
In previous articles, we have repeatedly introduced methods of cropping and batch cropping of raster remote sensing images on different platforms or based on different coding languages. Mainly include ArcPy in Python batch cropping of a large number of raster remote sensing images based on vector range (https://blog.csdn.net/zhebushibiaoshifu/article/details/128307676), Google Earth Engine Google Earth Engine GEE vector data cropping raster data (https ://blog.csdn.net/zhebushibiaoshifu/article/details/117390431), based on ENVI to implement raster remote sensing images to delimit and crop the rectangular research area according to the layer row number and number of pixels (https://blog.csdn. net/zhebushibiaoshifu/article/details/118978851), etc.; in this article, we will introduce the GDAL module based on **C++** language to implement batch cropping requirements.
First, let’s take a look at the requirements of this article. There is now a folder, as shown in the figure below, which contains multiple multi-band remote sensing image files in .tiff format (for convenience, we only have 2 files, but in fact our batch processing needs are generally much greater than this number).
What we hope to achieve is based on each remote sensing image in this folder, and its upper left corner100 * 100 pixel This part is cropped (as shown in the figure below) and saved as new remote sensing image files (the new file name just adds a _C suffix to the original file name), and Store in another specified results folder. We hope that the cropped remote sensing image, compared with the original remote sensing image, will appear as shown in the figure below.
The code used in this article is as follows.
#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;
}
Let’s introduce the specific content of the above code.
First, we need to register all GDAL drivers through GDALAllRegister();. At the same time, we defined the input and output folder paths – inputFolder is the folder path that stores input remote sensing images (remote sensing images to be cropped), outputFolder is the folder path where result remote sensing images are stored.
Secondly, we define a string list through CPLStringList fileList; to store the file list in the folder; and use the VSIReadDir function to read the file in the input folder all files and store the results in fileList.
Next, we use a loop to iterate through each file. First, get the path of the current file through string inputImagePath = fileList[i];; if the file extension is not tiff, skip the file. Next, for files with a file extension of tiff, we construct the complete input file path and use the GDALOpen function to open the input file and return a GDALDataset code> object, stored in poDataset.
Next, we can get the width and height of the input file and define the offset (the position of the upper left corner pixel), width and height of the cropping area. As mentioned before, what I need here is to crop this part of the 100 * 100 pixel from the upper left corner of the original remote sensing image (so the offset is 0) . Secondly, build the path to the output file and use the GetGDALDriverManager()->GetDriverByName function to obtain the GTiff driver object and store it in poDriver. Subsequently, we use the poDriver->Create function to create the output file, return a GDALDataset object, and store it in poOutputDataset.
This next part requires a little attention. First, we use poOutputDataset->SetProjection to set the projection information of the output file, which is the same projection as the input file; secondly, we use poDataset->GetGeoTransform to obtain the geographical transformation of the input file Parameters, stored in the adfGeoTransform array. Because in my case, the pixel size (i.e. the length and width of a single pixel) of the cropped remote sensing image has not changed, and the upper left corner pixel of the raster remote sensing image before and after cropping strong>There has been no change, so compared with the geographic transformation parameters of the new raster remote sensing image and the old raster remote sensing image, there is no need to Any changes; but if your cropping needs are not like this (for example, the pixel size has changed, or the cropping does not start from the upper left corner pixel), then You need to adjust these 6 geographical transformation parameters – as for how to change it, this is more complicated. I haven’t fully figured it out yet. Everyone can combine their actual needs and go to GDALYou can check it on the official website. Finally, we use poOutputDataset->SetGeoTransform to set the geographical transformation parameters of the output file. In my case, they are exactly the same as the input file.
At the end of the code, we use a loop to iteratively process each band (each remote sensing image here has a total of 4 bands). First, create a floating-point buffer with a size of xSize * ySize, and use poBand->RasterIO to read the pixel data of the corresponding band from the input file into the buffer area; next, use poOutputBand->RasterIO to write the data in the buffer to the corresponding band of the output file. The buffer memory can then be released and the output and input files can be closed.
Run the above code, we can see the cropped remote sensing image file in the result folder, and the file name of the new file also meets our requirements; as shown in the figure below.
At this point, you’re done.
Welcome to follow: Crazy Learning GIS