<html><div style='background-color:'><DIV class=RTE>
<P>Hello Frank,</P>
<P> Thank you for your very quick response, I changed the values like you said and the output file size is indeed almost the same as the original, just a couple hunderds kilobytes bigger, which makes sense since the transition from NAD83 to WGS84 should not have much effect. The squishing effect, however, is still present and is still greyscale.</P>
<P> For my experimentations I use visual C++ 6 and my processing is triggered by just one button click, as you will notice it is mostly your tutorial example from the main GDAL website, here is the whole code:</P>
<P>void CTest01Dlg::OnButtonProcess() <BR>{<BR> // ------------ Example of initializing an image into the GDALDataset object<BR> <BR> GDALDataset *poDataset;</P>
<P> // -------------- Reprojection example</P>
<P> GDALDatasetH hSrcDS, hDstDS;</P>
<P> // Open input and output files. </P>
<P> GDALAllRegister();</P>
<P> hSrcDS = GDALOpen( fileNameIn, GA_ReadOnly );</P>
<P> poDataset = (GDALDataset *) hSrcDS;<BR> if( poDataset == NULL )<BR> {<BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)strcat("The input file does not exist: ", fileNameIn));<BR> exit(1);<BR> }</P>
<P> GDALDriverH hDriver;<BR> GDALDataType eDT;</P>
<P> // Open the source file. </P>
<P> CPLAssert( hSrcDS != NULL );<BR> <BR> // Create output with same datatype as first input band. </P>
<P> eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));</P>
<P> //CString str;<BR> //str.Format("%d\n",eDT);<BR> //AfxMessageBox(str);</P>
<P> // Get output driver (GeoTIFF format)</P>
<P> hDriver = GDALGetDriverByName( "GTiff" );<BR> CPLAssert( hDriver != NULL );</P>
<P> // Get Source coordinate system. </P>
<P> const char *pszSrcWKT = NULL;<BR> char *pszDstWKT;</P>
<P> pszSrcWKT = GDALGetProjectionRef( hSrcDS );<BR> CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );</P>
<P> // Setup output coordinate system that is UTM 11 WGS84. </P>
<P> OGRSpatialReference oSRS;</P>
<P> oSRS.SetUTM( 17, TRUE );<BR> oSRS.SetWellKnownGeogCS( "WGS84" );</P>
<P> oSRS.exportToWkt( &pszDstWKT );</P>
<P> // Create a transformer that maps from source pixel/line coordinates<BR> // to destination georeferenced coordinates (not destination <BR> // pixel line). We do that by omitting the destination dataset<BR> // handle (setting it to NULL). </P>
<P> void *hTransformArg;</P>
<P> hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 );<BR> CPLAssert( hTransformArg != NULL );</P>
<P> // Get approximate output georeferenced bounds and resolution for file. </P>
<P> double adfDstGeoTransform[6];<BR> int nPixels = 0, nLines = 0;<BR> CPLErr eErr;</P>
<P> eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );<BR> CPLAssert( eErr == CE_None );</P>
<P> GDALDestroyGenImgProjTransformer( hTransformArg );</P>
<P> // Create the output file. <BR> hDstDS = GDALCreate( hDriver, fileNameOut, nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, NULL );<BR> <BR> CPLAssert( hDstDS != NULL );</P>
<P> // Write out the projection definition. </P>
<P> GDALSetProjection( hDstDS, pszDstWKT );<BR> GDALSetGeoTransform( hDstDS, adfDstGeoTransform );</P>
<P> // Copy the color table, if required.</P>
<P> GDALColorTableH hCT;</P>
<P> // Get the color table of the first (1) band of the source dataset<BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,1));<BR> <BR> // --------------- For debugging-----------------------<BR> //str.Format("DEBUG INFO: %d\n",GDALGetRasterCount(hDstDS));<BR> //AfxMessageBox(str);<BR> // ----------------------------------------------------<BR> <BR> // If the color table of the first (1) band of the source dataset is not<BR> // NULL then assign it to the first (1) band of the destination dataset<BR> //if(hCT != NULL)<BR> //{<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);<BR> //}</P>
<P> // Get the color table of the second (2) band of the source dataset<BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,2));</P>
<P> // If the color table of the second (2) band of the source dataset is not<BR> // NULL then assign it to the second (2) band of the destination dataset<BR> //if(hCT != NULL)<BR> //{<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,2), hCT);<BR> //}</P>
<P> // Get the color table of the third (3) band of the source dataset<BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,3));</P>
<P> // If the color table of the third (3) band of the source dataset is not<BR> // NULL then assign it to the third (3) band of the destination dataset<BR> //if(hCT != NULL)<BR> //{<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,3), hCT);<BR> //}</P>
<P> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Setup warp options ...");</P>
<P> // Setup warp options. <BR> <BR> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();</P>
<P> psWarpOptions->hSrcDS = hSrcDS;<BR> psWarpOptions->hDstDS = hDstDS;<BR> psWarpOptions->nBandCount = 3;<BR> psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);<BR> psWarpOptions->panSrcBands[0] = 1;<BR> psWarpOptions->panSrcBands[1] = 2;<BR> psWarpOptions->panSrcBands[2] = 3;<BR> psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);<BR> psWarpOptions->panDstBands[0] = 1;<BR> psWarpOptions->panDstBands[1] = 2;<BR> psWarpOptions->panDstBands[2] = 3;</P>
<P> psWarpOptions->pfnProgress = GDALTermProgress;<BR> <BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Establishing transformer ...");</P>
<P> // Establish reprojection transformer. </P>
<P> psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 0.0, 1);<BR> psWarpOptions->pfnTransformer = GDALGenImgProjTransform;</P>
<P> // Initialize and execute the warp operation.<BR> <BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Executing ...");</P>
<P> GDALWarpOperation oOperation;</P>
<P> // Initialize transformation operation<BR> oOperation.Initialize(psWarpOptions);</P>
<P> // Perform the transformation<BR> oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));</P>
<P> // Clean up the intermediate objects<BR> GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);<BR> GDALDestroyWarpOptions(psWarpOptions);</P>
<P> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Done ...");</P>
<P> // Close the input and output files<BR> GDALClose(hDstDS);<BR> GDALClose(hSrcDS); <BR>}</P>
<P>void CTest01Dlg::OnButtonCleanup() <BR>{<BR> remove(fileNameOut);<BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Cleanup complete ...");<BR>}</P>
<P>I also tried using the precompiled tools from the GDAL website, specifically the warp tool and it also produces the squished effect. Perhaps the GeoTIFF I am using as a source is special in some sense and needs to be processed in a special way, I am not sure. Here is the information extracted using gdalinfo.exe about the image, it is around 75 MB: </P>
<P>Driver: GTiff/GeoTIFF<BR>Size is 5010, 5009<BR>Coordinate System is:<BR>PROJCS["unnamed",<BR> GEOGCS["NAD83",<BR> DATUM["North_American_Datum_1983",<BR> SPHEROID["GRS 1980",6378137,298.2572221010002,<BR> AUTHORITY["EPSG","7019"]],<BR> AUTHORITY["EPSG","6269"]],<BR> PRIMEM["Greenwich",0],<BR> UNIT["degree",0.0174532925199433],<BR> AUTHORITY["EPSG","4269"]],<BR> PROJECTION["Transverse_Mercator"],<BR>
PARAMETER["latitude_of_origin",0],<BR> PARAMETER["central_meridian",-81],<BR> PARAMETER["scale_factor",0.9996],<BR> PARAMETER["false_easting",500000],<BR> PARAMETER["false_northing",0],<BR> UNIT["metre",1,<BR> AUTHORITY["EPSG","9001"]],<BR> AUTHORITY["EPSG","26917"]]<BR>Origin = (608999.000000,4845001.000000)<BR>Pixel Size = (0.20000000,-0.20000000)<BR>Corner Coordinates:<BR>Upper Left ( 608999.000, 4845001.000) ( 79d38'46.13"W, 43d45'0.13"N)<BR>Lower Left ( 608999.000, 4843999.200) ( 79d38'46.86"W, 43d44'27.66"N)<BR>Upper Right ( 610001.000, 4845001.000) ( 79d38'1.34"W, 43d44'59.59"N)<BR>Lower Right ( 610001.000, 4843999.200) ( 79d38'2.08"W, 43d44'27.13"N)<BR>Center (
609500.000, 4844500.100) ( 79d38'24.10"W, 43d44'43.63"N)<BR>Band 1 Block=5010x1 Type=Byte, ColorInterp=Red<BR>Band 2 Block=5010x1 Type=Byte, ColorInterp=Green<BR>Band 3 Block=5010x1 Type=Byte, ColorInterp=Blue<BR><BR> Thank you for your time and I really appreciate your help. Peter says hi, I am working as an undergrad student in GeoICT.</P>
<P> Ilya S.</P>
<P> </P>
<P> </P>
<P> </P></DIV>
<DIV></DIV>>From: Frank Warmerdam <fwarmerdam@gmail.com>
<DIV></DIV>>Reply-To: warmerdam@pobox.com
<DIV></DIV>>To: Ilya Serebrianik <predator2448@hotmail.com>
<DIV></DIV>>CC: gdal-dev@xserve.flids.com
<DIV></DIV>>Subject: Re: [Gdal-dev] Warping GeoTIFF colour image - beginner question(s)
<DIV></DIV>>Date: Mon, 2 May 2005 14:25:53 -0400
<DIV></DIV>>
<DIV></DIV>>On 5/2/05, Ilya Serebrianik <predator2448@hotmail.com> wrote:
<DIV></DIV>> > psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) *
<DIV></DIV>> > psWarpOptions->nBandCount);
<DIV></DIV>> > psWarpOptions->panSrcBands[0] = 1;
<DIV></DIV>> > psWarpOptions->panSrcBands[1] = 1;
<DIV></DIV>> > psWarpOptions->panSrcBands[2] = 1;
<DIV></DIV>> > psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) *
<DIV></DIV>> > psWarpOptions->nBandCount);
<DIV></DIV>> > psWarpOptions->panDstBands[0] = 1;
<DIV></DIV>> > psWarpOptions->panDstBands[1] = 1;
<DIV></DIV>> > psWarpOptions->panDstBands[2] = 1;
<DIV></DIV>>
<DIV></DIV>>Ilya,
<DIV></DIV>>
<DIV></DIV>>I think you want the src and dst bands to be:
<DIV></DIV>>
<DIV></DIV>> > psWarpOptions->panSrcBands[0] = 1;
<DIV></DIV>> > psWarpOptions->panSrcBands[1] = 2;
<DIV></DIV>> > psWarpOptions->panSrcBands[2] = 3;
<DIV></DIV>> > psWarpOptions->panDstBands[0] = 1;
<DIV></DIV>> > psWarpOptions->panDstBands[1] = 2;
<DIV></DIV>> > psWarpOptions->panDstBands[2] = 3;
<DIV></DIV>>
<DIV></DIV>>This would explain why you are getting greyscale output.
<DIV></DIV>>I'm less clear on this squishing effect you are seeing.
<DIV></DIV>>
<DIV></DIV>>Best regards,
<DIV></DIV>>
<DIV></DIV>>--
<DIV></DIV>>---------------------------------------+--------------------------------------
<DIV></DIV>>I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
<DIV></DIV>>light and sound - activate the windows | http://pobox.com/~warmerdam
<DIV></DIV>>and watch the world go round - Rush | Geospatial Programmer for Rent
<DIV></DIV></div></html>