GDAL对影像进行划分网格，并通过像素坐标转换为地理坐标

GDAL对影像进行划分网格，并通过像素坐标转换为地理坐标

#include "gdal_priv.h"
#include<iostream>  
#include<fstream>
using namespace std;

//这里通过结构体来定义从图像中获取到的像素点
typedef struct imgpoint {
	int id;
	int leftPointX;
	int leftPointY;
	double w;
	double h;
};

void ImageProcess(const char* pszSrcFile, const char* pszFormat)
{
	

	FILE *fp;
	fp = fopen("D:/dcj/test.txt", "r+");
	//打开一个txt文件，并将点的像素坐标和地理坐标写到文件中

	//注册GDAL驱动
	GDALAllRegister();

	//pszSrcFile = "G:/3Datas/wgs841_10w/H49D001004.tif";
	
	//打开输入图像
	GDALDataset *poSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
	
	//获取输入图像的宽高波段书
	//获取图像波段 
	GDALRasterBand *poBand;//*poBand1, *poBand2
	int bandcount;
	bandcount = poSrcDS->GetRasterCount();	// 获取波段数:DEM为1个波段
	cout << "波段数：" << bandcount << endl;
	
	
	

	//获取输入图像仿射变换6个参数
	double adfGeotransform[6] = { 0 };
	poSrcDS->GetGeoTransform(adfGeotransform);
	//获取输入图像空间参考
	const char* pszProj = poSrcDS->GetProjectionRef();

	GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
	if (poBand == NULL)    //获取输入文件中的波段失败
	{
		GDALClose((GDALDatasetH)poSrcDS);
		//return false;
	}

	

	int nBlockSize = 256;     //分块大小，可以自定义

	//分配输入分块缓存
	int nBands = poSrcDS->GetRasterCount();
	unsigned char *pSrcData = new unsigned char[nBlockSize*nBlockSize*nBands];
	

	//定义读取输入图像波段顺序
	int *pBandMaps = new int[nBands];
	for (int b = 0; b < nBands; b++)
		pBandMaps[b] = b + 1;

	//循环分块并进行处理
	vector<imgpoint> block;
	imgpoint sig;
	int nYSize = poSrcDS->GetRasterYSize();
	int id = 0;
	for (int i = 0; i < nYSize; i += nBlockSize)
	{
		int Yid = 0;
		
		int nXSize = poSrcDS->GetRasterXSize();
		for (int j = 0; j <nXSize; j += nBlockSize)
		{
			int xid = 0;
		
			
			//定义两个变量来保存分块大小
			int nXBK = nBlockSize;
			int nYBK = nBlockSize;

			//如果最下面和最右边的块不够256，剩下多少读取多少
			if (i + nBlockSize > nYSize)     //最下面的剩余块
				nYBK = nYSize - i;
			if (j + nBlockSize > nXSize)     //最右侧的剩余块
				nXBK = nXSize - j;

			//读取原始图像块
			poSrcDS->RasterIO(GF_Read, j, i, nXBK, nYBK, pSrcData, nXBK, nYBK, GDT_Byte, nBands, pBandMaps, 0, 0, 0, NULL);

			sig.id = id;
			sig.leftPointX = i;
			sig.leftPointY = j;
			sig.w = nBlockSize;
			sig.h = nBlockSize;
			block.push_back(sig);
			imgpoint sigs = block[id];

			double Xgeo, Ygeo;
			Xgeo = adfGeotransform[0] + sig.leftPointX * adfGeotransform[1] + sig.leftPointY * adfGeotransform[2];
			Ygeo = adfGeotransform[3] + sig.leftPointX * adfGeotransform[4] + sig.leftPointY * adfGeotransform[5];
			//把坐标信息存到文件中，写法不是很标准，需要修改
			fprintf(fp, "%d\t",id);
			fprintf(fp, "%d\t", sig.leftPointX);
			fprintf(fp, "%d\t", sig.leftPointY);
			fprintf(fp, "%lf\t", Xgeo);
			fprintf(fp, "%lf\n", Ygeo);
			cout << id << ":" << sig.leftPointX << "," << sig.leftPointY << "," << Xgeo << "," << Ygeo << endl;
			//fprintf(fp, );
			id++;

;			
		}
		
		
	}

	

	//释放申请的内存
	delete[]pSrcData;
	delete[]pBandMaps;

	//关闭原始图像和结果图像
	GDALClose((GDALDatasetH)poSrcDS);

	//return true;
}


int main()
{
	const char* pszFile;
	const char* outFile;
	GDALAllRegister();
	pszFile = "G:/3Datas/wgs841_10w/H49D001004.tif";//影像路径
	GDALDataset *poDataset = (GDALDataset*)GDALOpen(pszFile, GA_ReadOnly);
	GDALRasterBand *poBand = poDataset->GetRasterBand(1);
	int xsize = poBand->GetXSize();//这里是获取影像大小
	int ysize = poBand->GetYSize();
	cout << xsize << endl;
	cout << ysize << endl;

	const char *pszFormat = "Tif";

	//double adfGeoTransform[6];//仿射参数
	//poDataset->GetGeoTransform(adfGeoTransform);//获取参数
	//for (int i = 0; i < 6; i++) {
	//	cout << adfGeoTransform[i]<<endl;
	//}
	
	ImageProcess(pszFile,pszFormat);

	

	
	return 0;
}