<div dir="ltr"><div>I am attempting to tranform a GeoTiff from UTM14N to WGS84 using GDAL from VC++ 2008 and I am brand new to using GDAL.</div><div><br></div><div>When I do the transformation everything seems to go off correctly but when I compare my transformed image in Quantum GIS to the source image projected on the fly by Quantum GIS their image is slightly rotated and mine is straight. Theirs is correct, mine is not(I believe), yet I am not sure what I am doing incorrectly. Is it my affine translation parameters, is it that I am warping it improperly, did I calculate corners improperly or something else entirely?</div>
<div><br></div><div>The source image is a 3 band GeoTiff if that helps at all. </div><div><br></div><div>Here is my complete code as I have it right now:</div><div><br></div><div>#include "stdafx.h"</div><div>#include "gdalwarper.h"</div>
<div>#include "cpl_conv.h"</div><div>#include "gdal_priv.h"</div><div>#include <ogr_spatialref.h></div><div>#include "ogr_srs_api.h"</div><div><br></div><div>int _tmain(int argc, _TCHAR* argv[])</div>
<div>{</div><div><span class="" style="white-space:pre"> </span>GDALAllRegister();</div><div><br></div><div><span class="" style="white-space:pre"> </span>GDALDataset *poDataset, *poDstDS;</div><div><span class="" style="white-space:pre"> </span>const char *pszFormat = "GTiff";</div>
<div><span class="" style="white-space:pre"> </span>const char *pszSrcFilename = "C:\\Example\\SampleData\\In.tif";</div><div><span class="" style="white-space:pre"> </span>const char *pszDstFilename = "C:\\Example\\SampleData\\Out.tif";</div>
<div><span class="" style="white-space:pre"> </span>GDALDriver *poDriver;</div><div><span class="" style="white-space:pre"> </span>int rasterCount;</div><div><br></div><div><span class="" style="white-space:pre"> </span>OGRSpatialReference oSrcSRS,oDstSRS;</div>
<div><span class="" style="white-space:pre"> </span>char *pszSrcSRSWKT = NULL;</div><div><span class="" style="white-space:pre"> </span>char *pszDstSRSWKT = NULL;</div><div><br></div><div><span class="" style="white-space:pre"> </span>OGRCoordinateTransformation *transform = NULL;</div>
<div><br></div><div><span class="" style="white-space:pre"> </span>double adfSrcGeoTransform[6];</div><div><span class="" style="white-space:pre"> </span>double adfDstGeoTransform[6];</div><div><span class="" style="white-space:pre"> </span>double xS[2];</div>
<div><span class="" style="white-space:pre"> </span>double yS[2];</div><div><br></div><div><span class="" style="white-space:pre"> </span>//Set destination spatial reference to WGS84 and get its WKT</div><div><span class="" style="white-space:pre"> </span>oDstSRS.SetWellKnownGeogCS("WGS84");</div>
<div><span class="" style="white-space:pre"> </span>oDstSRS.exportToWkt(&pszDstSRSWKT);</div><div><br></div><div><span class="" style="white-space:pre"> </span>//Now set the source spatial reference to UTM14N WGS84 and get its WKT</div>
<div><span class="" style="white-space:pre"> </span>oSrcSRS.SetProjCS( "UTM 14 (WGS84) in northern hemisphere." );</div><div> <span class="" style="white-space:pre"> </span>oSrcSRS.SetWellKnownGeogCS( "WGS84" );</div>
<div> <span class="" style="white-space:pre"> </span>oSrcSRS.SetUTM( 14, TRUE );</div><div><span class="" style="white-space:pre"> </span>oSrcSRS.exportToWkt(&pszSrcSRSWKT);</div><div><br></div><div><span class="" style="white-space:pre"> </span>//Set the coordinate transformation</div>
<div><span class="" style="white-space:pre"> </span>transform = OGRCreateCoordinateTransformation(&oSrcSRS,&oDstSRS);</div><div><br></div><div><span class="" style="white-space:pre"> </span>//Set the GeoTiff driver</div>
<div><span class="" style="white-space:pre"> </span>poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);</div><div><br></div><div><span class="" style="white-space:pre"> </span>// Open input and output files. </div>
<div> <span class="" style="white-space:pre"> </span>poDataset = (GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );</div><div> <span class="" style="white-space:pre"> </span>rasterCount = poDataset->GetRasterCount();</div>
<div><br></div><div><span class="" style="white-space:pre"> </span>poDstDS = poDriver->CreateCopy( pszDstFilename, poDataset, FALSE, NULL, NULL, NULL );</div><div><br></div><div><span class="" style="white-space:pre"> </span></div>
<div> <span class="" style="white-space:pre"> </span>// Setup warp options. </div><div><span class="" style="white-space:pre"> </span>GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();</div><div><br></div><div><span class="" style="white-space:pre"> </span>psWarpOptions->hSrcDS = poDataset;</div>
<div><span class="" style="white-space:pre"> </span>psWarpOptions->hDstDS = poDstDS;</div><div><br></div><div><br></div><div><span class="" style="white-space:pre"> </span>psWarpOptions->nBandCount = rasterCount;</div>
<div><span class="" style="white-space:pre"> </span>psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );</div><div><span class="" style="white-space:pre"> </span>psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );</div>
<div><br></div><div><span class="" style="white-space:pre"> </span>for(int i=0;i<rasterCount;i++)</div><div><span class="" style="white-space:pre"> </span>{</div><div><span class="" style="white-space:pre"> </span>psWarpOptions->panSrcBands[i] = i+1;//src band number</div>
<div><span class="" style="white-space:pre"> </span>psWarpOptions->panDstBands[i] = i+1;//dst band number</div><div><span class="" style="white-space:pre"> </span>}</div><div><br></div><div><span class="" style="white-space:pre"> </span>psWarpOptions->pfnProgress = GDALTermProgress; </div>
<div><br></div><div><span class="" style="white-space:pre"> </span>// Establish reprojection transformer. </div><div><br></div><div><span class="" style="white-space:pre"> </span>psWarpOptions->pTransformerArg = </div>
<div><span class="" style="white-space:pre"> </span>GDALCreateGenImgProjTransformer( poDataset, </div><div><span class="" style="white-space:pre"> </span>pszSrcSRSWKT, </div><div><span class="" style="white-space:pre"> </span>poDstDS,</div>
<div><span class="" style="white-space:pre"> </span>pszDstSRSWKT, </div><div><span class="" style="white-space:pre"> </span>FALSE, 0.0, 1 );</div><div><br></div><div><span class="" style="white-space:pre"> </span>psWarpOptions->pfnTransformer = GDALGenImgProjTransform;</div>
<div><span class="" style="white-space:pre"> </span>psWarpOptions->eResampleAlg = GRA_Cubic;</div><div><br></div><div><span class="" style="white-space:pre"> </span>// Initialize and execute the warp operation. </div><div>
<span class="" style="white-space:pre"> </span>GDALWarpOperation oOperation;</div><div><br></div><div><span class="" style="white-space:pre"> </span>oOperation.Initialize( psWarpOptions );</div><div><span class="" style="white-space:pre"> </span>oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( poDstDS ), GDALGetRasterYSize( poDstDS ) );</div>
<div><br></div><div><span class="" style="white-space:pre"> </span>GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );</div><div><span class="" style="white-space:pre"> </span>GDALDestroyWarpOptions( psWarpOptions );</div>
<div><span class="" style="white-space:pre"> </span></div><div><br></div><div><span class="" style="white-space:pre"> </span>//Set the new image's projection tag info </div><div><span class="" style="white-space:pre"> </span>poDstDS->SetProjection(pszDstSRSWKT);</div>
<div><br></div><div><span class="" style="white-space:pre"> </span>//Get the original affine parameters</div><div><span class="" style="white-space:pre"> </span>poDataset->GetGeoTransform(adfSrcGeoTransform);</div><div>
<br></div><div><span class="" style="white-space:pre"> </span>double srcWidth = poDataset->GetRasterXSize();</div><div><span class="" style="white-space:pre"> </span>double srcHeight = poDataset->GetRasterYSize();</div>
<div><br></div><div><br></div><div><span class="" style="white-space:pre"> </span>//Get the corners</div><div><span class="" style="white-space:pre"> </span>xS[0] = adfSrcGeoTransform[0];//minx</div><div><span class="" style="white-space:pre"> </span>xS[1] = adfSrcGeoTransform[0] + srcWidth*adfSrcGeoTransform[1] + srcHeight*adfSrcGeoTransform[2];//??? maxx </div>
<div><span class="" style="white-space:pre"> </span> </div><div><span class="" style="white-space:pre"> </span>yS[0] = adfSrcGeoTransform[3];//maxy</div><div><span class="" style="white-space:pre"> </span>yS[1] = adfSrcGeoTransform[3] + srcWidth*adfSrcGeoTransform[4] + srcHeight*adfSrcGeoTransform[5];//??? miny</div>
<div><br></div><div><br></div><div><span class="" style="white-space:pre"> </span>//Reproject its corners coordinate info</div><div><span class="" style="white-space:pre"> </span>transform->Transform(2,xS,yS);</div><div>
<br></div><div><span class="" style="white-space:pre"> </span>//Set all the new affine parameters</div><div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[0] = xS[0];</div><div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[1] = (xS[1] - xS[0]) / srcWidth;//(xMax - xMin) / imageWidth </div>
<div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[2] = adfSrcGeoTransform[2];</div><div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[3] = yS[0];</div><div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[4] = adfSrcGeoTransform[4];</div>
<div><span class="" style="white-space:pre"> </span>adfDstGeoTransform[5] = ((yS[0] - yS[1]) / srcHeight)*(-1);//((yMax - yMin) / imageHeight)*(-1) </div><div><br></div><div><span class="" style="white-space:pre"> </span>//Set the new images affine information</div>
<div><span class="" style="white-space:pre"> </span>poDstDS->SetGeoTransform(adfDstGeoTransform);</div><div><br></div><div> <span class="" style="white-space:pre"> </span>GDALClose( poDstDS );</div><div> <span class="" style="white-space:pre"> </span>GDALClose( poDataset );</div>
<div><br></div><div> <span class="" style="white-space:pre"> </span>return 0;</div><div>}</div><div><br></div><div style>Thank you in advance for any help you folks may be able to give me on this.</div><div style><br></div>
<div style>-Andy</div></div>