使用GDAL对遥感影像进行投影转换

1.在进行投影转换前需要做的工作：

首先编译GDAL和proj4库,GDAL中进行投影转换依赖的开源库是proj4，可以到网上下载已经编译好的包含proj4的GDAL，CSDN上就有很多资源。接下来将编译好的GDAL和PROJ4的lib和include路径分别加在VS2010的“附加库目录”，和“附加包含目录”中；在连接器->附加依赖项中添加编译好的proj.lib文件和gdal_i.lib。然后将proj.dll拷到相应的debug文件下。

环境设置好之后，接下来进行投影转换的第一步。

2.投影转换：

我做的是同一椭球体基准面下的转换，将投影坐标转化为地理坐标。只要把原理和环境弄清楚了，相信其他的转换也很简单。

OGRSpatialReference Raster_spf;

OGRSpatialReference *Dest_Raster_spf;

Raster_spf=m_pDataset->GetProjectionRef();

if (Raster_spf.IsProjected())

{

char*pszSrcWKT=NULL;

Raster_spf.exportToWkt(&pszSrcWKT);

Dest_Raster_spf=Raster_spf.CloneGeogCS();

char *pszDestWKT=NULL;

Dest_Raster_spf->exportToPrettyWkt(&pszDestWKT);

OGRCoordinateTransformation *poCT;

poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);//反过来就是地理坐标转投影坐标

if (NULL==poCT)

{

QMessageBox::warning(this, tr("Failed"), tr("poCT is null!"), tr("OK"));

}

int ndbFlag=poCT->Transform(4,dbX,dbY,NULL);

if (ndbFlag)

{

OGRCoordinateTransformation::DestroyCT(poCT);

}

其中dbX和dbY存储的是影像四个顶点坐标。

3,转换后处理

转换完之后需要重新计算影像的地理六参数和行列号，然后再利用GDAL中gdalwarp这个接口进行重采样。这部分的内容是借鉴李民录的GDAL书籍和zhouxuguang236的CSDN博客，地址为/article/1386058.html。

首先，重新计算地理六参数和行列号

dbRes = dbRes / 111000; //投影坐标转地理坐标分辨率需除以111000m,反过来就得乘以111000m.原理参考上述博客

double dbMinx = 0;

double dbMaxx = 0;

double dbMiny = 0;

double dbMaxy = 0;

dbMinx = min(min(min(dbX[0],dbX[1]),dbX[2]),dbX[3]);

dbMaxx = max(max(max(dbX[0],dbX[1]),dbX[2]),dbX[3]);

dbMiny = min(min(min(dbY[0],dbY[1]),dbY[2]),dbY[3]);

dbMaxy = max(max(max(dbY[0],dbY[1]),dbY[2]),dbY[3]);

//估算行列号

padfTransform[0] = dbMinx; //左上角点坐标

padfTransform[3] = dbMaxy;

//padfTransform[1] 像素宽度, padfTransform[5]像素高度

padfTransform[1] = dbRes;

padfTransform[5] = -dbRes;

//估算行列数

int nPixels = ceil(fabs(dbMaxx-dbMinx)/dbRes);

int nLines = ceil(fabs(dbMaxy-dbMiny)/dbRes);

接下来，利用gdalwarp来重采样并将其写入新的影像中。

GDALDriver* pGDalDriver = NULL;//用于creat的驱动

GDALDataset * poDataset; //结果影像

GDALResampleAlg eResampleMethod=GRA_Bilinear; //默认设置采样方式为双线性插值，时间会稍短一些

int iDstWidth=nPixels;

int iDstHeight=nLines;

const char *pszFormat;

pszFormat= m_pDataset->GetDriver()->GetDescription();

pGDalDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

QString FilenameD;

FilenameD=QFileDialog::getSaveFileName(this,tr("Save Image"),"",tr("Images (*.tif)"));

if(FilenameD.isEmpty())

{ return; }

poDataset=pGDalDriver->Create(FilenameD.toStdString().c_str(),iDstWidth,iDstHeight,dataBands,dataType,NULL);

poDataset->SetGeoTransform(padfTransform);

poDataset->SetProjection(pszDestWKT);

//构造坐标转换关系

void *hTransformArg=NULL;

hTransformArg=GDALCreateGenImgProjTransformer((GDALDatasetH)m_pDataset,pszSrcWKT,(GDALDatasetH)poDataset,

pszDestWKT,FALSE,0.0,1);

GDALTransformerFunc pfnTransformer=GDALGenImgProjTransform;

//构造gdalwarp的变化选项

GDALWarpOptions *psWO=GDALCreateWarpOptions();

psWO->papszWarpOptions=CSLDuplicate(NULL);

psWO->eWorkingDataType=dataType;

psWO->eResampleAlg=eResampleMethod;

psWO->hSrcDS=(GDALDatasetH)m_pDataset;

psWO->hDstDS=(GDALDatasetH)poDataset;

psWO->pfnTransformer=pfnTransformer;

psWO->pTransformerArg=hTransformArg;

psWO->nBandCount=dataBands;

psWO->panSrcBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));

psWO->panDstBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));

psWO->panSrcBands[0]=1;

psWO->panDstBands[0]=1;

//创建GDALWarp执行对象并使用GDALWarpOptions来进行初始化

GDALWarpOperation oWo;

oWo.Initialize(psWO);

//执行处理,转换会耗点时间，请耐心等待

oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);

//释放和关闭文件

GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);

GDALDestroyWarpOptions(psWO);

GDALClose((GDALDatasetH)m_pDataset);

GDALClose((GDALDatasetH)poDataset);

到此，投影转换的整个过程就结束了。之前遇到的问题是在转换过程中执行poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);得到的poCT为NULL，这肯定就是跟proj4库有关了，后来才发现是环境设置的问题。还有在调试时，执行到oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);就不动了，原以为是转换除了问题，但跟踪发现结果影像中相关的信息都得到了，应该不会有错。后来发现是重采样过程比较慢，需要时间，但最邻近和双线性插值可能需要的时间稍短些。