<html><div style='background-color:'><DIV class=RTE>
<P>Hello Frank,</P>
<P> It does not fail in a sense that it crashes, it just doesn't do anything to the image as was the case before I added projection to be latlong. The pszSrcWKT, and pszDstWKT seem to correct, and the functions that compose the latter do not complain either. Here is the code I use, I apologize for the mess, don't want to erase experimental lines:</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; //, hRefDS;</P>
<P> // Open input and output files. </P>
<P> GDALAllRegister();</P>
<P> hSrcDS = GDALOpen( fileNameIn, GA_ReadOnly );<BR> //hRefDS = GDALOpen( fileNameRef, GA_ReadOnly );</P>
<P><BR> //AfxMessageBox(str.Format("%d\n",hSrcDS->GetGCPCount()));</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> //const char *pszRefWKT = NULL;<BR> char *pszDstWKT;</P>
<P> pszSrcWKT = GDALGetProjectionRef( hSrcDS );<BR> CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );<BR> <BR> //pszRefWKT = GDALGetProjectionRef( hRefDS );<BR> //CPLAssert( pszRefWKT != NULL && strlen(pszRefWKT) > 0 );</P>
<P> // Setup output coordinate system that is UTM 17 WGS84. </P>
<P> OGRSpatialReference oSRS;</P>
<P> //oSRS.SetUTM(17, TRUE);</P>
<P> //oSRS.SetWellKnownGeogCS("WGS84");<BR> oSRS.importFromProj4("+proj=latlong +datum=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> char *apszOptions[] = {"INTERLEAVE=PIXEL", NULL};<BR> hDstDS = GDALCreate(hDriver, fileNameOut, nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, apszOptions);<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> // --------------- For debugging-----------------------<BR> //str.Format("DEBUG INFO: %d\n",GDALGetRasterCount(hDstDS));<BR> //AfxMessageBox(str);<BR> // ----------------------------------------------------</P>
<P> GDALColorTableH hCT;</P>
<P> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,1));<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);</P>
<P> //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);</P>
<P> <BR> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,3));<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,3), hCT);</P>
<P><BR> CTest01Dlg::SendDlgItemMessage(IDC_LIST_STATUS, LB_ADDSTRING, 0, (LPARAM)"Setup warp options ...");</P>
<P> // Setup warp options. </P>
<P> 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> // GDALClose(hRefDS);<BR>}</P>
<P><BR> Thank you.</P>
<P> Ilya.</P>
<P> <BR></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] UTM to Lat/Long projection
<DIV></DIV>>Date: Wed, 4 May 2005 15:10:44 -0400
<DIV></DIV>>
<DIV></DIV>>On 5/4/05, Ilya Serebrianik <predator2448@hotmail.com> wrote:
<DIV></DIV>> >
<DIV></DIV>> > Hello Frank,
<DIV></DIV>> >
<DIV></DIV>> > Thanks for clearing this up. I am not sure how to project from UTM
<DIV></DIV>> > NAD83 to Lat/Long WGS84. In your example application, namely gdalwrap, there
<DIV></DIV>> > is a option "+proj=latlong" which would be specified for the -t_srs in order
<DIV></DIV>> > to go from whatever the projection of the source into lat and long, I tried
<DIV></DIV>> > that and it works, but I do not know how to set this up inside my code, I
<DIV></DIV>> > tried using importFromProj4("+proj=longlat"), which should
<DIV></DIV>> > automatically set the lat/long option into the srs object I am about to
<DIV></DIV>> > assign to the destination image, like this:
<DIV></DIV>>
<DIV></DIV>>Ilya,
<DIV></DIV>>
<DIV></DIV>>And did the example you gave not work? In what way does it fail?
<DIV></DIV>>Do the pszSrcWKT and pszDstWKT seem OK?
<DIV></DIV>>
<DIV></DIV>> > OGRSpatialReference oSRS;
<DIV></DIV>> > //oSRS.SetUTM(17, TRUE);
<DIV></DIV>> > //oSRS.SetWellKnownGeogCS("WGS84");
<DIV></DIV>> > oSRS.importFromProj4("+proj=longlat");
<DIV></DIV>>
<DIV></DIV>>I would encourage you to be specific about the datum. So
<DIV></DIV>>you could do:
<DIV></DIV>>
<DIV></DIV>> oSRS.importFromProj4( "+proj=latlong +datum=WGS84" );
<DIV></DIV>>
<DIV></DIV>>Best regards,
<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>