最近来到新的环境里，分配的第一个任务是学习gdal，对不同地区、不同影响分辨率的tif图像重采样再剪切生成相同大小、分辨率的tif影像。我用的编译器是VS2010，环境配置请自行参考网上的教程，这里就不在讨论。
对于相同分辨率的遥感影像，很容易对相同像元进行处理，难得是对不同分辨率的影像处理。
强大的gdal当然封装了可用的方法来实现，聪明的你肯定想到了，RasterIO函数可以实现啊，RasterIO函数怎样使用，请转到李明录老师的博客，老师写的很详细。
下面贴上我自己写的函数，完全是源码，写的足够详细，也配有注释，不足是没有考虑效率问题，大家耐住性子看一看吧
// GDAl_ResampleAndCut.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include "gdal_priv.h" #include "ogrsf_frmts.h" #include "gdalwarper.h" #include <iostream> using namespace std; int gdalResample(vector <char*> pszSrcFile,vector<char*> pszOutFile) { GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO"); int iCount = pszSrcFile.size(); pszSrcFile.resize(iCount); int i=0; int j=0; vector<GDALDataset* > pDSrc; pDSrc.resize(iCount); for (i=0;i<iCount;i++) { pDSrc[i] = (GDALDataset*)GDALOpen(pszSrcFile[i],GA_ReadOnly); } pszOutFile.resize(iCount); GDALDriver *poDriver=GetGDALDriverManager()->GetDriverByName("GTiff"); if (poDriver==NULL) { for ( j =0;j< iCount ; j++ ) { GDALClose((GDALDatasetH)pDSrc[j]); } return -2; } int *width=new int[iCount]; int *heigh=new int[iCount]; int *nBandCount = new int[iCount]; GDALDataType *dataType = new GDALDataType[iCount]; double (*dGeoTrans)[6]=new double[iCount][6]; vector<char *> pszWKT; pszWKT.resize(iCount); for ( j = 0 ; j < iCount ; j++ ) { dataType[j]=pDSrc[j]->GetRasterBand(1)->GetRasterDataType(); width[j]=pDSrc[j]->GetRasterXSize(); heigh[j]=pDSrc[j]->GetRasterYSize(); nBandCount[j]=pDSrc[j]->GetRasterCount(); pDSrc[j]->GetGeoTransform(dGeoTrans[j]);//获取放射变换系数 pszWKT[j] = const_cast<char*>(pDSrc[j]->GetProjectionRef()); } //重采样到栅格单元的最大分辨率 //横向分辨率。。。 double maxX = dGeoTrans[0][1]; for (i = 0;i<iCount;i++) { if (maxX<dGeoTrans[i][1]) { maxX = dGeoTrans[i][1]; } } //纵向分辨率 double maxY = dGeoTrans[0][5]; for (i=0;i<iCount;i++) { if ( maxY>dGeoTrans[i][5] ) { maxY = dGeoTrans[i][5]; } //纵向分辩率是负数，求最小值 } //求转换参数 double *transformer_X = new double[iCount]; double *transformer_Y = new double[iCount]; for (i = 0 ; i<iCount;i++) { transformer_X[i] = dGeoTrans[i][1]/maxX; transformer_Y[i] = dGeoTrans[i][5]/maxY; } //求重采样后的图像的尺寸 int *newWidth = new int [iCount]; int *newHeigh = new int [iCount]; for (i = 0;i<iCount;i++) { newWidth[i] = static_cast<int>(width[i] * transformer_X[i] + 0.5); newHeigh[i] = static_cast<int>(heigh[i] * transformer_Y[i] + 0.5); } vector<GDALDataset*> poDstDs ; i=0; poDstDs.resize(iCount); float *pdata = new float[10000*10000]; for (i = 0;i<iCount;i++) { poDstDs[i] = poDriver->Create(pszOutFile[i], newWidth[i], newHeigh[i],nBandCount[i],dataType[i],NULL); /*dGeoTrans[i][1] = maxX; dGeoTrans[i][5] = maxY * (-1.0);*/ //转换后的分辨率 dGeoTrans[i][1] = dGeoTrans[i][1]/transformer_X[i]; dGeoTrans[i][5] = dGeoTrans[i][5]/transformer_Y[i]; poDstDs[i]->SetProjection(pszWKT[i]); poDstDs[i]->SetGeoTransform(dGeoTrans[i]); pDSrc[i]->RasterIO(GF_Read,0,0, width[i] , heigh[i] , pdata , newWidth[i] , newHeigh[i],dataType[i],nBandCount[i],0,0,0,0 ); poDstDs[i]->RasterIO(GF_Write,0,0 , newWidth[i] , newHeigh[i] , pdata , newWidth[i] , newHeigh[i] , dataType[i],nBandCount[i],0,0,0,0); } for (i = 0;i<iCount;i++) { if (pDSrc[i]!= NULL) { GDALClose((GDALDatasetH)pDSrc[i]); } if (poDstDs[i]!=NULL) { GDALClose((GDALDatasetH)poDstDs[i]); } /*if (pszOutFile[i]!=NULL) { GDALClose((GDALDatasetH)pszOutFile[i]); }*/ } cout<<"Rsample Completed"<<endl; } int gdalCut(vector<char*> pszOutFile,vector<char*> pszCutFile) { GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO"); int iCount = pszOutFile.size(); pszOutFile.resize(iCount); int i = 0; int j = 0; GDALDriver *poDriver=GetGDALDriverManager()->GetDriverByName("GTiff"); //求得icount个影像的角点的 地理坐标 //左上角点 double *GeoLeftTop_X = new double[iCount]; double *GeoLeftTop_Y = new double[iCount]; //左下角点 double *GeoLeftButtom_X = new double[iCount]; double *GeoLeftButtom_Y = new double[iCount]; //右上 double *GeoRightTop_X = new double[iCount]; double *GeoRightTop_Y = new double[iCount]; double (*dGeoTrans)[6]=new double[iCount][6]; vector<char *> pszWKT; pszWKT.resize(iCount); int *width=new int[iCount]; int *heigh=new int[iCount]; int *nBandCount = new int[iCount]; GDALDataType *dataType = new GDALDataType[iCount]; vector<GDALDataset* > poDstDs; poDstDs.resize(iCount); for (i=0;i<iCount;i++) { poDstDs[i] = (GDALDataset*)GDALOpen(pszOutFile[i],GA_ReadOnly); } for ( j = 0 ; j < iCount ; j++ ) { dataType[j]=poDstDs[j]->GetRasterBand(1)->GetRasterDataType(); width[j]=poDstDs[j]->GetRasterXSize(); heigh[j]=poDstDs[j]->GetRasterYSize(); nBandCount[j]=poDstDs[j]->GetRasterCount(); poDstDs[j]->GetGeoTransform(dGeoTrans[j]);//获取放射变换系数 pszWKT[j] = const_cast<char*>(poDstDs[j]->GetProjectionRef()); } //由图像行列坐标求地理坐标 for (i = 0;i < iCount;i ++) { GeoLeftTop_X[i] = dGeoTrans[i][0]; GeoLeftTop_Y[i] = dGeoTrans[i][3]; GeoRightTop_X[i] = (dGeoTrans[i][0] + dGeoTrans[i][1] * width[i] + dGeoTrans[i][2] * 0); GeoRightTop_Y[i] = (dGeoTrans[i][3] + dGeoTrans[i][4] * width[i] + dGeoTrans[i][5] * 0); GeoLeftButtom_X[i] = (dGeoTrans[i][0] + dGeoTrans[i][1] * 0 + dGeoTrans[i][2] * heigh[i]); GeoLeftButtom_Y[i] = (dGeoTrans[i][3] + dGeoTrans[i][4] * 0 + dGeoTrans[i][5] * heigh[i]); } //求各个角点的（x,y）坐标最小值 double GeoLeftTop_minX = GeoLeftTop_X[0],GeoLeftTop_minY =GeoLeftTop_Y[0]; double GeoLeftButtom_minX = GeoLeftButtom_X[0],GeoLeftButtom_minY = GeoLeftButtom_Y[0]; double GeoRightTop_minX = GeoRightTop_X[0],GeoRightTop_minY = GeoRightTop_Y[0]; for (int i=1;i<iCount;i++) { //左边的X找最大值 if ( GeoLeftTop_X[i] > GeoLeftTop_minX ) { GeoLeftTop_minX = GeoLeftTop_X[i]; } //上边的Y找最小值 if ( GeoLeftTop_Y[i] < GeoLeftTop_minY ) { GeoLeftTop_minY = GeoLeftTop_Y[i]; } //左边的X找最大值 if (GeoLeftButtom_X[i] > GeoLeftButtom_minX) { GeoLeftButtom_minX = GeoLeftButtom_X[i]; } //下边的Y找最大值 if (GeoLeftButtom_Y[i] > GeoLeftButtom_minY) { GeoLeftButtom_minY = GeoLeftButtom_Y[i]; } //右边的X找最小值 if (GeoRightTop_X[i] < GeoRightTop_minX) { GeoRightTop_minX = GeoRightTop_X[i]; } //上边的Y找最小值 if (GeoRightTop_Y[i] < GeoRightTop_minY) { GeoRightTop_minY = GeoRightTop_Y[i]; } } //地理坐标最小值（minx,miny）反算至图像行列坐标 /*double translator = 0.0; int tifLeftTop_X = 0,tifLeftTop_Y=0,tifLeftButtom_X = 0,tifLeftButtom_Y = 0,tifRightTop_X = 0,tifRightTop_Y = 0;*/ /*translator = dGeoTrans[0][1]*dGeoTrans[0][5] - dGeoTrans[0][2]*dGeoTrans[0][4]; tifLeftTop_X = (int)(( dGeoTrans[0][5] * (GeoLeftTop_minX-dGeoTrans[0][0]) - dGeoTrans[0][2]*( GeoLeftTop_minY - dGeoTrans[0][3]) )/translator + 0.5); tifLeftTop_Y = (int)(( dGeoTrans[0][1] * (GeoLeftTop_minY-dGeoTrans[0][3]) - dGeoTrans[0][4]*( GeoLeftTop_minX - dGeoTrans[0][0]) )/translator + 0.5);*/ /*tifLeftButtom_X = (int)(( dGeoTrans[0][5] * (GeoLeftButtom_minX-dGeoTrans[0][0]) - dGeoTrans[0][2]*( GeoLeftButtom_minY - dGeoTrans[0][3]) )/translator + 0.5); tifLeftButtom_Y = (int)(( dGeoTrans[0][1] * (GeoLeftButtom_minY-dGeoTrans[0][3]) - dGeoTrans[0][4]*( GeoLeftButtom_minX - dGeoTrans[0][0]) )/translator + 0.5); tifRightTop_X = (int)(( dGeoTrans[0][5] * (GeoRightTop_minX-dGeoTrans[0][0]) - dGeoTrans[0][2]*( GeoRightTop_minY- dGeoTrans[0][3]) )/translator + 0.5); tifRightTop_Y = (int)(( dGeoTrans[0][1] * (GeoRightTop_minY-dGeoTrans[0][3]) - dGeoTrans[0][4]*( GeoRightTop_minX - dGeoTrans[0][0]) )/translator + 0.5);*/ /*int widthCut = (int)(tifRightTop_X - tifLeftTop_X); int heighCut = (int)(tifLeftButtom_Y - tifLeftTop_Y);*/ double *translator = new double[iCount]; int *tifLeftTop_X = new int[iCount]; int *tifLeftTop_Y = new int[iCount]; int *tifLeftButtom_X = new int[iCount]; int *tifLeftButtom_Y = new int[iCount]; int *tifRightTop_X = new int[iCount]; int *tifRightTop_Y = new int[iCount]; for (i=0;i<iCount;i++) { translator[i] = dGeoTrans[i][1] * dGeoTrans[i][5] - dGeoTrans[i][2] * dGeoTrans[i][4]; tifLeftTop_X[i] = (int)(( dGeoTrans[i][5] * (GeoLeftTop_minX-dGeoTrans[i][0]) - dGeoTrans[i][2]*( GeoLeftTop_minY - dGeoTrans[i][3]) )/translator[i] + 0.5); tifLeftTop_Y[i] = (int)(( dGeoTrans[i][1] * (GeoLeftTop_minY-dGeoTrans[i][3]) - dGeoTrans[i][4]*( GeoLeftTop_minX - dGeoTrans[i][0]) )/translator[i] + 0.5); tifLeftButtom_X[i] = (int)(( dGeoTrans[i][5] * (GeoLeftButtom_minX-dGeoTrans[i][0]) - dGeoTrans[i][2]*( GeoLeftButtom_minY - dGeoTrans[i][3]) )/translator[i] + 0.5); tifLeftButtom_Y[i] = (int)(( dGeoTrans[i][1] * (GeoLeftButtom_minY-dGeoTrans[i][3]) - dGeoTrans[i][4]*( GeoLeftButtom_minX - dGeoTrans[i][0]) )/translator[i] + 0.5); tifRightTop_X[i] = (int)(( dGeoTrans[i][5] * (GeoRightTop_minX-dGeoTrans[i][0]) - dGeoTrans[i][2]*( GeoRightTop_minY- dGeoTrans[i][3]) )/translator[i] + 0.5); tifRightTop_Y[i] = (int)(( dGeoTrans[i][1] * (GeoRightTop_minY-dGeoTrans[i][3]) - dGeoTrans[i][4]*( GeoRightTop_minX - dGeoTrans[i][0]) )/translator[i] + 0.5); } int *widthCut = new int[iCount]; int *heighCut = new int[iCount]; for (i = 0;i<iCount;i++) { widthCut[i] = (int)(tifRightTop_X[i] - tifLeftTop_X[i]); heighCut[i] = (int)(tifLeftButtom_Y[i] - tifLeftTop_Y[i]); } vector<GDALDataset*> pszResampleDst; pszResampleDst.resize(iCount); for (i = 0;i<iCount;i++) { pszResampleDst[i] = poDriver->Create(pszCutFile[i],widthCut[i],heighCut[i],nBandCount[i],dataType[i],NULL); pszResampleDst[i]->SetGeoTransform(dGeoTrans[i]); pszResampleDst[i]->SetProjection(pszWKT[i]); } //GDALDataset *pszCutDstDs = poDriver->Create(pszCutFile,widthCut,heighCut,1,dataType[0],NULL); //pszCutDstDs->SetGeoTransform(dGeoTrans[0]); //pszCutDstDs->SetProjection(pszWKT[0]); float *pdata2 = new float[10000*10000]; for (i = 0;i<iCount;i++) { poDstDs[i]->RasterIO(GF_Read,(int)tifLeftTop_X[i],(int)tifLeftTop_Y[i],widthCut[i],heighCut[i],pdata2,widthCut[i],heighCut[i],dataType[i],nBandCount[i],0,0,0,0); pszResampleDst[i]->RasterIO(GF_Write,0,0,widthCut[i],heighCut[i],pdata2,widthCut[i],heighCut[i],dataType[i],nBandCount[i],0,0,0,0); } //poDstDs[0]->RasterIO(GF_Read,(int)tifLeftTop_X,(int)tifLeftTop_Y,widthCut,heighCut,pdata2,widthCut,heighCut,dataType[0],nBandCount[0],0,0,0,0); //pszCutDstDs->RasterIO(GF_Write); cout<<"Cut Completed!"<<endl; return 0; } int _tmain(int argc, _TCHAR* argv[]) { vector<char *> pszSrcFile(4); pszSrcFile[0]="D:\\Globe.tif"; pszSrcFile[1]="D:\\Asia.tif"; pszSrcFile[2]="D:\\South Sea.tif"; pszSrcFile[3]="D:\\A.tif"; int count=pszSrcFile.size(); vector<char *> pszOutFile; pszOutFile.resize(4); pszOutFile[0]="D:\\Resample_Globe.tif" ; pszOutFile[1]="D:\\Resample_Asia.tif" ; pszOutFile[2]="D:\\Resample_South Sea.tif" ; pszOutFile[3]="D:\\Resample_A.tif"; gdalResample(pszSrcFile,pszOutFile); vector<char *> pszResampleFile; pszResampleFile.resize(4); pszResampleFile[0] = "D:\\Cut_Globe.tif"; pszResampleFile[1] = "D:\\Cut_Asia.tif"; pszResampleFile[2] = "D:\\Cut_South Sea.tif"; pszResampleFile[3] = "D:\\Cut_A.tif"; gdalCut(pszOutFile,pszResampleFile); system("pause"); return 0; } 计算采用的四张影像，由于保密需要，这里就不贴上来啦 ， 但从名字上来看，也可以猜出是怎样的区域范围。
有不懂的地方请在下面评论，科研不忙的时候，我会回复大家的啦~，毕竟帮助别人，快乐自己嘛，嘿嘿~
