[gdal-dev] Image(GeoTiff) not rotated properly?

[Andy Canfield]

I am attempting to tranform a GeoTiff from UTM14N to WGS84 using GDAL from
VC++ 2008 and I am brand new to using GDAL.

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?

The source image is a 3 band GeoTiff if that helps at all.

Here is my complete code as I have it right now:

#include "stdafx.h"
#include "gdalwarper.h"
#include "cpl_conv.h"
#include "gdal_priv.h"
#include <ogr_spatialref.h>
#include "ogr_srs_api.h"

int _tmain(int argc, _TCHAR* argv[])
{
GDALAllRegister();

GDALDataset  *poDataset, *poDstDS;
const char *pszFormat = "GTiff";
const char *pszSrcFilename = "C:\\Example\\SampleData\\In.tif";
const char *pszDstFilename = "C:\\Example\\SampleData\\Out.tif";
GDALDriver *poDriver;
int rasterCount;

OGRSpatialReference oSrcSRS,oDstSRS;
char *pszSrcSRSWKT = NULL;
char *pszDstSRSWKT = NULL;

OGRCoordinateTransformation *transform = NULL;

double adfSrcGeoTransform[6];
double adfDstGeoTransform[6];
double xS[2];
double yS[2];

//Set destination spatial reference to WGS84 and get its WKT
oDstSRS.SetWellKnownGeogCS("WGS84");
oDstSRS.exportToWkt(&pszDstSRSWKT);

//Now set the source spatial reference to UTM14N WGS84 and get its WKT
oSrcSRS.SetProjCS( "UTM 14 (WGS84) in northern hemisphere." );
    oSrcSRS.SetWellKnownGeogCS( "WGS84" );
    oSrcSRS.SetUTM( 14, TRUE );
oSrcSRS.exportToWkt(&pszSrcSRSWKT);

//Set the coordinate transformation
transform = OGRCreateCoordinateTransformation(&oSrcSRS,&oDstSRS);

//Set the GeoTiff driver
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

// Open input and output files.
    poDataset = (GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );
    rasterCount = poDataset->GetRasterCount();

poDstDS = poDriver->CreateCopy( pszDstFilename, poDataset, FALSE, NULL,
NULL, NULL );

    // Setup warp options.
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

psWarpOptions->hSrcDS = poDataset;
psWarpOptions->hDstDS = poDstDS;


psWarpOptions->nBandCount = rasterCount;
psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) *
psWarpOptions->nBandCount );
psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) *
psWarpOptions->nBandCount );

for(int i=0;i<rasterCount;i++)
{
psWarpOptions->panSrcBands[i] = i+1;//src band number
psWarpOptions->panDstBands[i] = i+1;//dst band number
}

psWarpOptions->pfnProgress = GDALTermProgress;

// Establish reprojection transformer.

psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( poDataset,
pszSrcSRSWKT,
poDstDS,
pszDstSRSWKT,
FALSE, 0.0, 1 );

psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
psWarpOptions->eResampleAlg = GRA_Cubic;

// Initialize and execute the warp operation.
GDALWarpOperation oOperation;

oOperation.Initialize( psWarpOptions );
oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( poDstDS ),
GDALGetRasterYSize( poDstDS ) );

GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );

//Set the new image's projection tag info
poDstDS->SetProjection(pszDstSRSWKT);

//Get the original affine parameters
poDataset->GetGeoTransform(adfSrcGeoTransform);

double srcWidth = poDataset->GetRasterXSize();
double srcHeight = poDataset->GetRasterYSize();


//Get the corners
xS[0] = adfSrcGeoTransform[0];//minx
xS[1] = adfSrcGeoTransform[0] + srcWidth*adfSrcGeoTransform[1] +
srcHeight*adfSrcGeoTransform[2];//??? maxx

yS[0] = adfSrcGeoTransform[3];//maxy
yS[1] = adfSrcGeoTransform[3] + srcWidth*adfSrcGeoTransform[4] +
srcHeight*adfSrcGeoTransform[5];//??? miny


//Reproject its corners coordinate info
transform->Transform(2,xS,yS);

//Set all the new affine parameters
adfDstGeoTransform[0] = xS[0];
adfDstGeoTransform[1] = (xS[1] - xS[0]) / srcWidth;//(xMax - xMin) /
imageWidth
adfDstGeoTransform[2] = adfSrcGeoTransform[2];
adfDstGeoTransform[3] = yS[0];
adfDstGeoTransform[4] = adfSrcGeoTransform[4];
adfDstGeoTransform[5] = ((yS[0] - yS[1]) / srcHeight)*(-1);//((yMax - yMin)
/ imageHeight)*(-1)

//Set the new images affine information
poDstDS->SetGeoTransform(adfDstGeoTransform);

    GDALClose( poDstDS );
    GDALClose( poDataset );

    return 0;
}

Thank you in advance for any help you folks may be able to give me on this.

-Andy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20130325/b91ff2db/attachment.html>