GDAL变换类重采样变形

GDALDriverH hDriver; GDALDataType eDT; GDALDatasetH hDstDS; GDALDatasetH hSrcDS; GDALRasterBandH hSrcbandDS; // Open the source file. GDALAllRegister(); hSrcDS = GDALOpen( "D://7.tif", GA_ReadOnly ); CPLAssert( hSrcDS != NULL ); // Create output with same datatype as first input band. hSrcbandDS = GDALGetRasterBand(hSrcDS,1); eDT = GDALGetRasterDataType(hSrcbandDS); // Get output driver (GeoTIFF format) hDriver = GDALGetDriverByName( "GTiff" ); CPLAssert( hDriver != NULL ); // Get Source coordinate system. const char *pszSrcWKT; char*pszDstWKT = NULL; pszSrcWKT = GDALGetProjectionRef( hSrcDS ); CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 ); // Get approximate output georeferenced bounds and resolution for file. double adfDstGeoTransform[6]; int nPixels=0, nLines=0; CPLErr eErr; nPixels=6000; nLines=6000; // Create the output file. hDstDS = GDALCreate( hDriver, "D://11.tif", nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, NULL ); CPLAssert( hDstDS != NULL ); // Write out the projection definition. // Setup output coordinate system that is UTM 11 WGS84. OGRSpatialReference oSRS; oSRS.SetUTM( 50, TRUE ); oSRS.SetWellKnownGeogCS( "WGS84" ); oSRS.exportToWkt(&pszDstWKT ); // Create a transformer that maps from source pixel/line coordinates // to destination georeferenced coordinates (not destination // pixel line). We do that by omitting the destination dataset // handle (setting it to NULL). void *hTransformArg; hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 1000, 1 ); eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines ); CPLAssert( hTransformArg != NULL ); GDALDestroyGenImgProjTransformer( hTransformArg ); adfDstGeoTransform[1]=adfDstGeoTransform[1]/2; adfDstGeoTransform[5]=adfDstGeoTransform[5]/2; CPLAssert( eErr == CE_None ); GDALSetProjection( hDstDS, pszDstWKT ); GDALSetGeoTransform( hDstDS, adfDstGeoTransform ); // Copy the color table, if required. GDALColorTableH hCT; hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) ); if( hCT != NULL ) GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT ); // Setup warp options. GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->hSrcDS = hSrcDS; psWarpOptions->hDstDS = hDstDS; psWarpOptions->nBandCount = 3; psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); psWarpOptions->panSrcBands[0] = 1; psWarpOptions->panSrcBands[1] = 2; psWarpOptions->panSrcBands[2] = 3; psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); psWarpOptions->panDstBands[0] = 1; psWarpOptions->panDstBands[1] = 2; psWarpOptions->panDstBands[2] = 3; psWarpOptions->pfnProgress = GDALTermProgress; // Establish reprojection transformer. psWarpOptions->eResampleAlg=GRA_NearestNeighbour;// psWarpOptions->pTransformerArg =// GDALCreateGenImgProjTransformer( hSrcDS,// GDALGetProjectionRef(hSrcDS),// hDstDS,// GDALGetProjectionRef(hDstDS),// FALSE, 1000.0, 3 );// psWarpOptions->pfnTransformer = GDALGenImgProjTransform; void * m_hGenTransformArg=NULL; void * m_hApproxArg=NULL; m_hGenTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 1000.0, 3 ); if( m_hGenTransformArg == NULL ) return CE_Failure; m_hApproxArg = GDALCreateApproxTransformer( GDALGenImgProjTransform, m_hGenTransformArg, 0.125 ); if ( m_hApproxArg == NULL) return CE_Failure; psWarpOptions->pfnTransformer = GDALApproxTransform; psWarpOptions->pTransformerArg = m_hApproxArg; // Initialize and execute the warp operation. GDALWarpOperation oOperation; oOperation.Initialize( psWarpOptions );// int Typesize=GDALGetDataTypeSize(GDT_Byte);// int Dstwidth,DstHeight;// int Srcwidth,Srcheight;// Dstwidth=500;// DstHeight=500;// Srcheight=1000;// Srcwidth=1000;// void *pDstBuffer = new byte[ Typesize* Srcheight* Srcwidth*GDALGetRasterCount(hDstDS) ];// memset(pDstBuffer, 0,Typesize* Srcheight* Srcwidth*GDALGetRasterCount(hDstDS));// int nJpgBandCount = GDALGetRasterCount(hDstDS);// int *panBands = new int[nJpgBandCount];// for (int i = 0; i < nJpgBandCount; i++)// {// panBands[i] = i+1;// }//// //AfxMessageBox("123");// CPLErr err = GDALDatasetRasterIO( hSrcDS, GF_Read,// 0, 0, Srcwidth, Srcheight,// pDstBuffer, Srcwidth, Srcheight,// GDT_Byte, nJpgBandCount, panBands,// 0, 0, 0 );//// oOperation.WarpRegionToBuffer(0,0,Dstwidth,DstHeight,pDstBuffer,GDT_Byte);//// eErr = GDALDatasetRasterIO( hDstDS, GF_Write,// 0, 0, Dstwidth, DstHeight,// pDstBuffer, Srcwidth, Srcheight,// GDT_Byte,// 3,// panBands,// 0, 0, 0 );// delete[] panBands; oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS)); GDALDestroyGenImgProjTransformer( m_hGenTransformArg ); GDALDestroyWarpOptions( psWarpOptions ); GDALDestroyApproxTransformer(m_hApproxArg); GDALClose( hDstDS ); GDALClose( hSrcDS ); return TRUE; BOOL CMapDivisionDoc::Interpolation(CString strFileTemp,CString strFile){ GDALDriverH hDriver; GDALDataType eDT; GDALDatasetH hDstDS; GDALDatasetH hSrcDS; GDALRasterBandH hSrcbandDS; // Open the source file. GDALAllRegister(); hSrcDS = GDALOpen(strFileTemp, GA_ReadOnly ); CPLAssert( hSrcDS != NULL ); // Create output with same datatype as first input band. hSrcbandDS = GDALGetRasterBand(hSrcDS,1); eDT = GDALGetRasterDataType(hSrcbandDS); // Get output driver (GeoTIFF format) char *pszFormat=GetOutimgInfo(strFileTemp); hDriver = GDALGetDriverByName(pszFormat); CPLAssert( hDriver != NULL ); // Get Source coordinate system. const char *pszSrcWKT; char*pszDstWKT = NULL; pszSrcWKT = GDALGetProjectionRef( hSrcDS ); CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 ); // Get approximate output georeferenced bounds and resolution for file. double adfDstGeoTransform[6]; int nPixels=0, nLines=0; CPLErr eErr; nPixels=m_nFileSize; nLines=m_nFileSize; // Create the output file. hDstDS = GDALCreate( hDriver, strFile, nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, NULL ); CPLAssert( hDstDS != NULL ); // Write out the projection definition. // Setup output coordinate system that is UTM 11 WGS84.// OGRSpatialReference oSRS;// oSRS.SetUTM( 50, TRUE );// oSRS.SetWellKnownGeogCS( "WGS84" );// oSRS.exportToWkt(&pszDstWKT ); // Create a transformer that maps from source pixel/line coordinates // to destination georeferenced coordinates (not destination // pixel line). We do that by omitting the destination dataset // handle (setting it to NULL). void *hTransformArg; hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszSrcWKT, FALSE, 1000, 1 ); eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines ); CPLAssert( hTransformArg != NULL ); GDALDestroyGenImgProjTransformer( hTransformArg ); adfDstGeoTransform[1]=adfDstGeoTransform[1]*GDALGetRasterXSize(hSrcDS)/m_nFileSize; adfDstGeoTransform[5]=adfDstGeoTransform[5]*GDALGetRasterXSize(hSrcDS)/m_nFileSize; CPLAssert( eErr == CE_None ); GDALSetProjection( hDstDS, pszSrcWKT ); GDALSetGeoTransform( hDstDS, adfDstGeoTransform ); // Copy the color table, if required. GDALColorTableH hCT; hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) ); if( hCT != NULL ) GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT ); // Setup warp options. GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->hSrcDS = hSrcDS; psWarpOptions->hDstDS = hDstDS; psWarpOptions->nBandCount = GDALGetRasterCount(hDstDS) ; psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); for (int i=0;i<GDALGetRasterCount(hDstDS);i++) { psWarpOptions->panSrcBands[i] = i+1; }// psWarpOptions->panSrcBands[0] = 1;// psWarpOptions->panSrcBands[1] = 2;// psWarpOptions->panSrcBands[2] = 3; psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); for (int m=0;m<GDALGetRasterCount(hDstDS);m++) { psWarpOptions->panDstBands[m] = m+1; }// psWarpOptions->panDstBands[0] = 1;// psWarpOptions->panDstBands[1] = 2;// psWarpOptions->panDstBands[2] = 3; psWarpOptions->pfnProgress = GDALTermProgress; // Establish reprojection transformer. psWarpOptions->eResampleAlg=GRA_Bilinear; void * m_hGenTransformArg=NULL; void * m_hApproxArg=NULL; m_hGenTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 1000.0, 3 ); if( m_hGenTransformArg == NULL ) return CE_Failure; m_hApproxArg = GDALCreateApproxTransformer( GDALGenImgProjTransform, m_hGenTransformArg, 0.125 ); if ( m_hApproxArg == NULL) return CE_Failure; psWarpOptions->pfnTransformer = GDALApproxTransform; psWarpOptions->pTransformerArg = m_hApproxArg; // Initialize and execute the warp operation. GDALWarpOperation oOperation; oOperation.Initialize( psWarpOptions ); oOperation.ChunkAndWarpMulti( 0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS)); GDALDestroyGenImgProjTransformer( m_hGenTransformArg ); GDALDestroyWarpOptions( psWarpOptions ); GDALDestroyApproxTransformer(m_hApproxArg); GDALClose( hDstDS ); GDALClose( hSrcDS ); remove(strFileTemp); return TRUE;}