<html><div style='background-color:'><DIV class=RTE>Hello,</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE> I have a problem with my reprojection example program, I am new to GDAL and still getting familiar with the its functionality. A couple of days back I could not get the program output a coloured image, but it turned out that I should have set INTERLEAVE=PIXEL, thanks to Frank I fixed that part. I thought that the reprojection occured and the output image was no different from input because the transition between NAD83 to WGS84 causes almost no distortion on larger scale. Today, however, I wanted to verify this by using OpenEV tool written, as I understand, in Python. I took my source image file, which is NAD83, then I took another proper geotiff which is in WGS84 and extracted the projection information from it to use in my transformation, here are the input and output WKTs:</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>INPUT:</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>PROJCS["NAD83 / UTM zone 17N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","26917"]]</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>OUTPUT:</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235630016,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE> I used the above in the Compose Dataset function in OpenEV and using the geotransform reprojection I converted the original geotiff from NAD83 to WGS84, the result came out to be what I expected from the beginning, which is a slightly rotated image and slightly squished from top to bottom. I do not know how OpenEV uses GDAL to make the correct transformation, but in my case it does not work. If someone can help me I would really appreciate it. Thank you very much ahead. I pasted the code below. I could make screenshots if that would help.</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE> Ilya.</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>P.S. The following is the code of my button click on a dialog form in vc++ 6, very simple:</DIV>
<DIV class=RTE> </DIV>
<DIV class=RTE>void CTest01Dlg::OnButtonProcess() <BR>{<BR> // ------------ Example of initializing an image into the GDALDataset object<BR> <BR> GDALDataset *poDataset;</DIV>
<DIV class=RTE> // -------------- Reprojection example</DIV>
<DIV class=RTE> GDALDatasetH hSrcDS, hDstDS, hRefDS;</DIV>
<DIV class=RTE> // Open input file. </DIV>
<DIV class=RTE> GDALAllRegister();</DIV>
<DIV class=RTE> hSrcDS = GDALOpen( fileNameIn, GA_ReadOnly );</DIV>
<DIV class=RTE>// Open the reference file, from which projection WKT would be used for output file.<BR> hRefDS = GDALOpen( fileNameRef, GA_ReadOnly );</DIV>
<DIV class=RTE><BR> //AfxMessageBox(str.Format("%d\n",hSrcDS->GetGCPCount()));</DIV>
<DIV class=RTE> 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> }</DIV>
<DIV class=RTE> GDALDriverH hDriver;<BR> GDALDataType eDT;</DIV>
<DIV class=RTE> // Open the source file. </DIV>
<DIV class=RTE> CPLAssert( hSrcDS != NULL );<BR> <BR> // Create output with same datatype as first input band. </DIV>
<DIV class=RTE> eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));</DIV>
<DIV class=RTE> //CString str;<BR> //str.Format("%d\n",eDT);<BR> //AfxMessageBox(str);</DIV>
<DIV class=RTE> // Get output driver (GeoTIFF format)</DIV>
<DIV class=RTE> hDriver = GDALGetDriverByName("GTiff");<BR> CPLAssert( hDriver != NULL );</DIV>
<DIV class=RTE> // Get Source coordinate system. </DIV>
<DIV class=RTE> const char *pszSrcWKT = NULL;<BR> const char *pszRefWKT = NULL;<BR> char *pszDstWKT;</DIV>
<DIV class=RTE> pszSrcWKT = GDALGetProjectionRef( hSrcDS );<BR> CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );<BR> <BR> pszRefWKT = GDALGetProjectionRef( hRefDS );<BR> CPLAssert( pszRefWKT != NULL && strlen(pszRefWKT) > 0 );</DIV>
<DIV class=RTE> // Setup output coordinate system that is UTM 17 WGS84. </DIV>
<DIV class=RTE> OGRSpatialReference oSRS;</DIV>
<DIV class=RTE> //oSRS.SetUTM(17, TRUE);<BR> //oSRS.SetWellKnownGeogCS("WGS84");</DIV>
<DIV class=RTE> oSRS.exportToWkt( &pszDstWKT );</DIV>
<DIV class=RTE> // 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). </DIV>
<DIV class=RTE> void *hTransformArg;</DIV>
<DIV class=RTE> hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszRefWKT, FALSE, 0, 1 );<BR> CPLAssert( hTransformArg != NULL );</DIV>
<DIV class=RTE> // Get approximate output georeferenced bounds and resolution for file. </DIV>
<DIV class=RTE> double adfDstGeoTransform[6];<BR> int nPixels = 0, nLines = 0;<BR> CPLErr eErr;</DIV>
<DIV class=RTE> eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );<BR> CPLAssert( eErr == CE_None );</DIV>
<DIV class=RTE> GDALDestroyGenImgProjTransformer( hTransformArg );</DIV>
<DIV class=RTE> // Create the output file.<BR> char *apszOptions[] = {"INTERLEAVE=PIXEL", NULL};<BR> hDstDS = GDALCreate(hDriver, fileNameOut, nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, apszOptions);<BR> <BR> CPLAssert( hDstDS != NULL );</DIV>
<DIV class=RTE> // Write out the projection definition. </DIV>
<DIV class=RTE> GDALSetProjection(hDstDS, pszRefWKT);<BR> GDALSetGeoTransform(hDstDS, adfDstGeoTransform);</DIV>
<DIV class=RTE> // Copy the color table, if required.</DIV>
<DIV class=RTE> // --------------- For debugging-----------------------<BR> //str.Format("DEBUG INFO: %d\n",GDALGetRasterCount(hDstDS));<BR> //AfxMessageBox(str);<BR> // ----------------------------------------------------</DIV>
<DIV class=RTE> GDALColorTableH hCT;</DIV>
<DIV class=RTE> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,1));<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);</DIV>
<DIV class=RTE> //str.Format("DEBUG INFO: \n %d \n ",GDALGetRasterCount(hDstDS));<BR> //AfxMessageBox(eRBCI.*/<BR> <BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,2));<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,2), hCT);</DIV>
<DIV class=RTE> <BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,3));<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,3), hCT);</DIV>
<DIV class=RTE><BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Setup warp options ...");</DIV>
<DIV class=RTE> // Setup warp options. </DIV>
<DIV class=RTE> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();</DIV>
<DIV class=RTE> 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;</DIV>
<DIV class=RTE> psWarpOptions->pfnProgress = GDALTermProgress;<BR> <BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Establishing transformer ...");</DIV>
<DIV class=RTE> // Establish reprojection transformer. </DIV>
<DIV class=RTE> psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 0.0, 1);<BR> psWarpOptions->pfnTransformer = GDALGenImgProjTransform;</DIV>
<DIV class=RTE> // Initialize and execute the warp operation.<BR> <BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Executing ...");</DIV>
<DIV class=RTE> GDALWarpOperation oOperation;</DIV>
<DIV class=RTE> // Initialize transformation operation<BR> oOperation.Initialize(psWarpOptions);</DIV>
<DIV class=RTE> // Perform the transformation<BR> oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));</DIV>
<DIV class=RTE> // Clean up the intermediate objects<BR> GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);<BR> GDALDestroyWarpOptions(psWarpOptions);</DIV>
<DIV class=RTE> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Done ...");</DIV>
<DIV class=RTE> // Close the input and output files<BR> GDALClose(hDstDS);<BR> GDALClose(hSrcDS);<BR> GDALClose(hRefDS);<BR>}</DIV>
<DIV class=RTE>void CTest01Dlg::OnButtonCleanup() <BR>{<BR> remove(fileNameOut);<BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Cleanup complete ...");<BR>}</DIV></div></html>