GDAL影像合并

// merge.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

#include "../3dparty/include/gdal_priv.h"
#include "../3dparty/include/gdal_alg.h"
#include "../3dparty/include/gdalwarper.h"

#include <cassert>

#ifdef _DEBUG
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#else
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#endif

#define SIZEOFARRAY(array) sizeof(array)/sizeof(array[0])

static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double dfXRes=0.0, dfYRes=0.0;

static int bTargetAlignedPixels = FALSE;
static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
static int bVRT = FALSE;

static char *pszDstFileName = "../data/out.tif";
static char *apszSrcFiles[] = 
{
	/*"../data/H50G038013DOM.tif",
	"../data/H50G038014DOM.tif",
	"../data/H50G038015DOM.tif",
	"../data/H50G039012DOM.tif",
	"../data/H50G039013DOM.tif",
	"../data/H50G039014DOM.tif",*/
	"../data/H50G039015DOM.tif",
	"../data/H50G039016DOM.tif"
};

GDALDataType eWorkingType = GDT_Byte;

static GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
static GDALDatasetH GDALWarpCreateOutput( char **papszSrcFiles, 
										 int nFilesCount, 
										 const char *pszFilename, 
										 const char *pszFormat, 
										 GDALDataType eDT );

static void CreateMask(char* const* paFilePath, const int nFilesCount);

int _tmain(int argc, _TCHAR* argv[])
{
	GDALAllRegister();

	GDALDatasetH hSrcDS = NULL;
	GDALDatasetH hDstDS = NULL;

	int nFilesCount = SIZEOFARRAY(apszSrcFiles);
	hDstDS = GDALWarpCreateOutput(apszSrcFiles, nFilesCount,
		pszDstFileName, "GTiff", GDT_Byte);

	char **papszSrcFiles = apszSrcFiles;
	CreateMask(papszSrcFiles, nFilesCount);


	for( int iSrc = 0; iSrc < nFilesCount; iSrc++ )
	{
		GDALDatasetH hSrcDS;

		hSrcDS = GDALOpen( apszSrcFiles[iSrc], GA_ReadOnly );
		GDALDataset *pSrcDS = (GDALDataset *)hSrcDS;
		assert(pSrcDS->GetRasterBand(1)->GetMaskBand() != NULL);

		if( hSrcDS == NULL )
			exit( 2 );

		if ( GDALGetRasterCount(hSrcDS) == 0 )
		{
			fprintf(stderr, "Input file %s has no raster bands.\n", apszSrcFiles[iSrc] );
			exit( 1 );
		}

		if( !bQuiet )
			printf( "Processing input file %s.\n", apszSrcFiles[iSrc] );

		if ( eResampleAlg != GRA_NearestNeighbour &&
			GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
		{
			if( !bQuiet )
				fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
				"bad results when using a resampling method other than "
				"nearest neighbour. Converting the dataset prior to 24/32 bit "
				"is advised.\n", apszSrcFiles[iSrc] );
		}

		if( GDALGetRasterColorInterpretation( 
			GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) 
			== GCI_AlphaBand 
			&& !bEnableSrcAlpha )
		{
			bEnableSrcAlpha = TRUE;
			if( !bQuiet )
				printf( "Using band %d of source image as alpha.\n", 
				GDALGetRasterCount(hSrcDS) );
		}
		void *hTransformArg = NULL;

		hTransformArg = GD