[gdal-dev] TIFF write and comparison

[Akhil Jaggarwal]

Hi,

I just started using the GDAL C API.

I wrote some code that compares pixel values of 4 equally sized
GeoTIFFs, and depending on the values, writes the output to an output
TIFF.

I'd like to ask if:

1) I'm doing it the correct way, specifically if I'm writing the
pixels the correctly.
2) Is using CreateCopy() an ideal method to write projection
information? Should I write using Create() instead? tiffinfo output of
the output tiff does not state all the projection information using
CreateCopy().

Thanks.

Here's the code:

#include "/usr/include/gdal/gdal.h"
#include "/usr/include/gdal/cpl_conv.h" /* for CPLMalloc() */
#include "/usr/include/gdal/cpl_string.h"	// for createcopy()


int main(int argc, char **argv)
{
	GDALDatasetH  hDataset;

	GDALDatasetH  hDataset2;
	GDALDatasetH  hDataset3;
	GDALDatasetH  hDataset4;

	// output tiff
	GDALDatasetH hDestDS;

	GDALAllRegister();

	hDataset = GDALOpen( argv[1], GA_ReadOnly );
	hDataset2 = GDALOpen( argv[2], GA_ReadOnly );
	hDataset3 = GDALOpen( argv[3], GA_ReadOnly );
	hDataset4 = GDALOpen( argv[4], GA_ReadOnly );

	GDALDriverH   hDriver;
	GDALDriverH   hDriver2;
	GDALDriverH   hDriver3;
	GDALDriverH   hDriver4;

	// array to hold projection data
	double        adfGeoTransform[6];

	hDriver = GDALGetDatasetDriver( hDataset );
	hDriver2 = GDALGetDatasetDriver( hDataset2 );
	hDriver3 = GDALGetDatasetDriver( hDataset3 );
	hDriver4 = GDALGetDatasetDriver( hDataset4 );
	//	hDriver = GDALGetDatasetDriver( hDataset2 );
	//	hDriver = GDALGetDatasetDriver( hDataset3 );
	//	hDriver = GDALGetDatasetDriver( hDataset4 );
	printf( "Driver: %s/%s\n",
			GDALGetDriverShortName( hDriver ),
			GDALGetDriverLongName( hDriver ) );
	/* ====================== get TIFF attributes ========================= */
	printf( "Size is %dx%dx%d\n",
			GDALGetRasterXSize( hDataset ),
			GDALGetRasterYSize( hDataset ),
			GDALGetRasterCount( hDataset ) );

	/*	printf( "Size is %dx%dx%d\n",
		GDALGetRasterXSize( hDataset2 ),
		GDALGetRasterYSize( hDataset2 ),
		GDALGetRasterCount( hDataset2 ) );
		*/
	if( GDALGetProjectionRef( hDataset ) != NULL )
		printf( "Projection (dataset1) is `%s'\n", GDALGetProjectionRef( hDataset ) );

	/*	if( GDALGetProjectionRef( hDataset2 ) != NULL )
		printf( "Projection (dataset2) is `%s'\n", GDALGetProjectionRef(
hDataset2 ) );
		*/
	if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
	{
		printf( "Origin = (%.6f,%.6f)\n",
				adfGeoTransform[0], adfGeoTransform[3] );

		printf( "Pixel Size = (%.6f,%.6f)\n",
				adfGeoTransform[1], adfGeoTransform[5] );
	}
	/* ================= Fetch raster data to read on later ================*/
	// fetching raster data
	//
	GDALRasterBandH hBand;
	GDALRasterBandH hBand2;
	GDALRasterBandH hBand3;
	GDALRasterBandH hBand4;
	int             nBlockXSize, nBlockYSize;
	int             bGotMin, bGotMax;
	double          adfMinMax[2];

	hBand = GDALGetRasterBand( hDataset, 1 );
	hBand2 = GDALGetRasterBand( hDataset, 1 );
	hBand3 = GDALGetRasterBand( hDataset, 1 );
	hBand4 = GDALGetRasterBand( hDataset, 1 );
	
	GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
	printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
			nBlockXSize, nBlockYSize,
			GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
			GDALGetColorInterpretationName(
				GDALGetRasterColorInterpretation(hBand)) );

	adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
	adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
	if( ! (bGotMin && bGotMax) )
		GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );

	printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );

	if( GDALGetOverviewCount(hBand) > 0 )
		printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand));

	if( GDALGetRasterColorTable( hBand ) != NULL )
		printf( "Band has a color table with %d entries.\n",
				GDALGetColorEntryCount(
					GDALGetRasterColorTable( hBand ) ) );

	/* ================== check if format supports createcopy() ===== */

	const char *pszFormat = "GTiff";
	hDriver = GDALGetDriverByName( pszFormat );
	char ** papszMetadata;

	if (hDriver == NULL )
		exit(1);

	papszMetadata = GDALGetMetadata( hDriver, NULL );
	
	if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) ) {
			printf( "Driver %s supports Create() method.\n", pszFormat );

	}
		
	if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) ) {
		printf( "Driver %s supports CreateCopy() method.\n", pszFormat );
	}
	/* ================== Read raster data ========================*/

	// for class tiff
	unsigned char * pafScanline;
	// for slope tiff
	short int * pafScanline2;
	// for elevation tiff
	short int * pafScanline3;
	// for wb1 tiff
	short int * pafScanline4;

	const char *pszDestFileName = "out_142029.tif";
	hDestDS = GDALCreateCopy(hDriver, pszDestFileName, hDataset, FALSE,
NULL, NULL, NULL );
	int nxSize = GDALGetRasterBandXSize( hBand );

	int nySize = GDALGetRasterBandYSize( hBand );
	printf("raster band xSize: %d\n", nxSize);
	printf("raster band ySize: %d\n", nySize);

	int rule1 = 5;
	int rule2 = 3;
        // pafScanline is a void pointer
	void *change1 = &rule1;
	void *change2 = &rule2;

	// allocate buffer to hold pixel data
	//
	// classification
	pafScanline = (unsigned char *)CPLMalloc(sizeof(unsigned char)*nxSize);
	// slope
	pafScanline2 = (short int *)CPLMalloc(sizeof(short int)*nxSize);
	// dem-elevation
	pafScanline3 = (short int *)CPLMalloc(sizeof(short int)*nxSize);
	// wb1
	pafScanline4 = (short int *)CPLMalloc(sizeof(short int)*nxSize);
	// erwflag, nxoff, nyoff, nxsize, nysize, pdata, nbufxsize, nbufysize,
	// ebuftype, npixelspace, nlinespace
	//

			//nBlockXSize, nBlockYSize,
	int i,j;
        // compare and read datasets, change values accordingly
	for(i = 0; i < nxSize; i++) {
		for(j = nySize; j < 5; j++) {
			GDALRasterIO(hBand, GF_Read, 0,0, i, j,	pafScanline,
			i, j,
			GDT_Byte,
			0, 0);
			GDALRasterIO(hBand, GF_Read, 0,0, i, j,	pafScanline2,
			i, j,
			GDT_Int16,
			0, 0);
			GDALRasterIO(hBand, GF_Read, 0,0, i, j,	pafScanline3,
			i, j,
			GDT_Int16,
			0, 0);
			GDALRasterIO(hBand, GF_Read, 0,0, i, j,	pafScanline4,
			i, j,
			GDT_Int16,
			0, 0);

			// read datasets, compare data and write
			if ( (pafScanline2[j] > 6 || pafScanline3[j] > 950) &&
					pafScanline[j] == 2) {
				GDALRasterIO( hDestDS, GF_Write, 0,0, i,j,
						change1, i, i, GDT_Byte, 0,0);
			}
			else if( (pafScanline4[j] < -700) &&
					pafScanline[j] == 6 ) {
				GDALRasterIO( hDestDS, GF_Write, 0,0, i,j,
						change2, i, i, GDT_Byte, 0,0);
			}
			else if( (pafScanline2[j] > 6) && pafScanline[j] == 3 ) {
				GDALRasterIO( hDestDS, GF_Write, 0,0, i,j,
						change2, i, i, GDT_Byte, 0,0);
			}
			else if((pafScanline3[j]>1500) && pafScanline[j] == 6) {
				GDALRasterIO( hDestDS, GF_Write, 0,0, i,j,
						change1, i, i, GDT_Byte, 0,0);
			}
			else
				GDALRasterIO( hDestDS, GF_Write, 0,0, i,j,
						&pafScanline[j], i, i, GDT_Byte, 0,0);
		//	printf("%u\n", pafScanline[j]);
		}
	}

	//	printf("width of pixels in this band: %d\n", );

	CPLFree(pafScanline);
	CPLFree(pafScanline2);
	CPLFree(pafScanline3);
	CPLFree(pafScanline4);
	
	GDALClose(hDataset);
	GDALClose(hDataset2);
	GDALClose(hDataset3);
	GDALClose(hDataset4);	
	GDALClose(hDestDS);

//	CSLDestroy(argv);
	GDALDestroyDriverManager();
	
	return 0;

}