利用gdal实现重采样与裁剪

最近来到新的环境里，分配的第一个任务是学习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;} 

计算采用的四张影像，由于保密需要，这里就不贴上来啦，但从名字上来看，也可以猜出是怎样的区域范围。

有不懂的地方请在下面评论，科研不忙的时候，我会回复大家的啦~，毕竟帮助别人，快乐自己嘛，嘿嘿~