Код: Выделить всё
#include <iostream>
#include "gdal_priv.h"
#include "cpl_conv.h" // for CPLMalloc()
#include "gdalwarper.h"
#include "ogr_spatialref.h"
using namespace std;
int main()
{
const char * src_filename = "./xxxx.tiff";
const char * dst_filename = "./xxxx_OUT.tiff";
CPLSetConfigOption( "CPL_LOG", "./log.txt" );
CPLSetConfigOption( "CPL_LOG_ERRORS", "ON" );
CPLSetConfigOption( "GDAL_DATA", "C:/gdal/gdal-data" );
// Open input and output files.
GDALAllRegister();
GDALDriver* hDriver;
GDALDataType eDT;
GDALDataset* hDstDS;
GDALDataset* hSrcDS;
// Open the source file.
hSrcDS = (GDALDataset*) GDALOpen( src_filename, GA_ReadOnly );
CPLAssert( hSrcDS != NULL );
// Create output with same datatype as first input band.
eDT = hSrcDS->GetRasterBand(1)->GetRasterDataType();
// Get output driver (GeoTIFF format)
hDriver = (GDALDriver*) GDALGetDriverByName( "GTiff" );
CPLAssert( hDriver != NULL );
// Get Source coordinate system.
char *pszDstWKT = NULL;
const char * pszSrcWKT = hSrcDS->GetProjectionRef();
CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );
// Setup output coordinate system that is UTM 54 WGS84.
OGRSpatialReference oSRS;
oSRS.SetUTM( 54, 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 = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 );
CPLAssert( hTransformArg != NULL );
// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
CPLErr eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );
CPLAssert( eErr == CE_None );
GDALDestroyGenImgProjTransformer( hTransformArg );
// Create the output file.
hDstDS = (GDALDataset*) GDALCreate( hDriver, dst_filename, nPixels, nLines, hSrcDS->GetRasterCount(), eDT, NULL );
/* OUTS */
CPLAssert( hDstDS != NULL );
// Write out the projection definition.
GDALSetProjection( hDstDS, pszDstWKT );
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
// Copy the color table, if required.
GDALColorTableH hCT[3];
hCT[0] = GDALGetRasterColorTable( hSrcDS->GetRasterBand(1) );
hCT[1] = GDALGetRasterColorTable( hSrcDS->GetRasterBand(2) );
hCT[2] = GDALGetRasterColorTable( hSrcDS->GetRasterBand(3) );
if( hCT[0] != NULL )
{
GDALSetRasterColorTable( hDstDS->GetRasterBand(1), hCT[0] );
GDALSetRasterColorTable( hDstDS->GetRasterBand(2), hCT[1] );
GDALSetRasterColorTable( hDstDS->GetRasterBand(3), hCT[2] );
}
// Setup warp options.
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
psWarpOptions->eWorkingDataType = eDT;
psWarpOptions->eResampleAlg = GRA_Cubic;
psWarpOptions->nBandCount = hSrcDS->GetRasterCount();
psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
for(int i = 0; i < psWarpOptions->nBandCount; ++i)
{
psWarpOptions->panSrcBands = i+1;
psWarpOptions->panDstBands = i+1;
}
psWarpOptions->pfnProgress = GDALTermProgress;
// Establish reprojection transformer.
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( hSrcDS,
GDALGetProjectionRef(hSrcDS),
hDstDS,
GDALGetProjectionRef(hDstDS),
FALSE, 0.0, 1 );
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// Initialize and execute the warp operation.
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions );
oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );
// End
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
return 0;
}Отлично работает с 8-битными геотифами, но при 16-битном выдает на выходе 16-битноке (как и нужно) но чернобелое одноканальное-изображение (т.е. RGB8 -> Gray16).
Пробовал воспользоваться утилитой gdalwarp - итог тот же: при работе с 16-битным геотифам (либо с параметром "-ot UInt16") выдает одноканал.
В чем дело? как это исправить?
