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

[Andy Canfield]

Please disregard I found what I needed to know in here:
http://www.gdal.org/warptut.html

and in here:
http://trac.osgeo.org/gdal/browser/trunk/gdal/alg/gdalwarper.cpp

and I now have it coming out correctly tilted.

Thank You,
Andy Canfield


On Mon, Mar 25, 2013 at 5:51 PM, Andy Canfield <andy.canfield at gmail.com>wrote:

> 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/20130326/91b8151e/attachment-0001.html>