利用GDAL提取正射影像中的感兴趣区域

以前知道有GDAL用于遥感影像的处理，这回第一次真正的用上了，任务如标题。当然，专门的遥感软件解决这个问题是分分钟的事，这里利用C++进行了实现。废话不多说，代码如下：

/* * domPath：DOM输入路径 * leftTx，leftTy左上角坐标 * rightBx，rightBy右下角坐标 * roiPath：roi输出路径 */bool GetROI(const string domPath, const double leftTx, const double leftTy, const double rightBx, const double rightBy, const string roiPath){GDALAllRegister();CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");//打开数据pSrcDS = (GDALDataset *)GDALOpen(domPath.data(), GA_ReadOnly);if(pSrcDS==NULL)return false;//获取格式GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");if (pDriver == NULL){GDALClose((GDALDatasetH) pSrcDS);return false;}int dXSize = pSrcDS->GetRasterXSize();int dYSize = pSrcDS->GetRasterYSize();int nBandCount = pSrcDS->GetRasterCount();//波段数//获取第一波段的数据类型GDALDataType dataType = pSrcDS->GetRasterBand(1)->GetRasterDataType();//获取投影信息char *pszWkt = const_cast<char *>(pSrcDS->GetProjectionRef());//如果没有投影，人为设置一个if(strlen(pszWkt)<=0){OGRSpatialReference oSRS;oSRS.SetUTM(50,true);   //北半球  东经120度oSRS.SetWellKnownGeogCS("WGS84");oSRS.exportToWkt(&pszWkt);}void *hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH) pSrcDS,pszWkt,NULL,pszWkt,false,0.0,1);if (hTransformArg == NULL){GDALClose((GDALDatasetH) pSrcDS);return false;}GDALDestroyGenImgProjTransformer(hTransformArg);double adfGeoTransform[6] = {0};double newGeoTransform[6] = {0};if(pSrcDS->GetGeoTransform( adfGeoTransform ) != CE_None)return false;memcpy(newGeoTransform, adfGeoTransform, sizeof(double)*6);//根据输入的感兴趣区域平面坐标转换为像点坐标//计算左顶点int leftT_u = 0, leftT_v = 0, rightB_u = 0, rightB_v = 0;Projection2Image(adfGeoTransform, leftTx, leftTy, leftT_u, leftT_v);double leftTop_proX, leftTop_proY;Image2Projection(adfGeoTransform, leftT_u, leftT_v, leftTop_proX, leftTop_proY);newGeoTransform[0] = leftTop_proX;newGeoTransform[3] = leftTop_proY;//计算右底点Projection2Image(adfGeoTransform, rightBx, rightBy, rightB_u, rightB_v);//宽度和高度int roiWidth = rightB_u - leftT_u;int roiHeigh = rightB_v - leftT_v;//创建输出文件GDALDataset* pDstDS = pDriver->Create(roiPath.data(), roiWidth, roiHeigh, nBandCount, dataType, NULL);if (pDstDS == NULL){GDALClose((GDALDatasetH) pSrcDS);return false;}    pDstDS->SetProjection(pszWkt);    pDstDS->SetGeoTransform(newGeoTransform);    //从原图像上读取感兴趣区域    BYTE *ppafScan= new BYTE [roiWidth * roiHeigh *nBandCount];    pSrcDS->RasterIO(GF_Read, leftT_u, leftT_v, roiWidth, roiHeigh, ppafScan, roiWidth, roiHeigh, dataType, nBandCount,0,0,0,0 );    //将缓存图像写入新图像    pDstDS->RasterIO(GF_Write, 0, 0, roiWidth, roiHeigh, ppafScan, roiWidth, roiHeigh, dataType, nBandCount,0,0,0,0 );    GDALFlushCache( pDstDS );    GDALClose((GDALDatasetH) pSrcDS);    GDALClose((GDALDatasetH) pDstDS);return true;}void Projection2Image(double *adfGeoTransform, double proX, double proY, int &iCol, int &iRow){   try   {double dTemp = adfGeoTransform[1]*adfGeoTransform[5] - adfGeoTransform[2]*adfGeoTransform[4];double dCol = 0.0, dRow = 0.0;dCol = (adfGeoTransform[5]*(proX - adfGeoTransform[0]) - adfGeoTransform[2]*(proY - adfGeoTransform[3])) / dTemp + 0.5;dRow = (adfGeoTransform[1]*(proY - adfGeoTransform[3]) - adfGeoTransform[4]*(proX - adfGeoTransform[0])) / dTemp + 0.5;iCol = static_cast<int>(dCol);iRow = static_cast<int>(dRow);return;}catch(...){return;}}void Image2Projection(double *adfGeoTransform, int iCol, int iRow, double &proX, double &proY){try{proX = adfGeoTransform[0] + adfGeoTransform[1] * iCol + adfGeoTransform[2] * iRow;proY = adfGeoTransform[3] + adfGeoTransform[4] * iCol + adfGeoTransform[5] * iRow;return;}catch(...){return;}}