[Gdal-dev] Using the warp api
Matt Lynch
matt at terraEngine.com
Wed Aug 27 18:27:35 EDT 2003
Hi,
I have been working my way through the warp api tutorial
(http://www.remotesensing.org/gdal/warptut.html), but have run into
trouble.
My program takes an image and reprojects it to a new file out.tif. When
I run it, I don't get any errors. However, out.tif is essentially
black, with a few lines. If my description of the output is not
sufficient, I can email the file out.tif to you.
Any help would be most appreciated.
Regards,
Matt
This is the run of my code:
C:\matt\gdal-cvs-030403\apps>test
Creating output file is that 9358P x 9357L.
Output file was successfully created.
Warp options are set.
transforming from
PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.
2572235630016,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM
["Greenw
ich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJ
ECTION["
Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["centr
al_merid
ian",-115],PARAMETER["scale_factor",0.9992],PARAMETER["false_easting",50
0000],PA
RAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
to
PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.
257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EP
SG","632
6"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174
53292519
9433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHO
RITY["EP
SG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standa
rd_paral
lel_1",35],PARAMETER["standard_parallel_2",55],PARAMETER["latitude_of_or
igin",45
],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAME
TER["fal
se_northing",0]]
Transformer is built.
Progress:0...10...20...30...40...50...60...70...80...90...100 - done.
C:\matt\gdal-cvs-030403\apps>
I also did gdalinfo on the two images in case that helps:
C:\matt\gdal-cvs-030403\apps>gdalinfo daImage.tif
Driver: GTiff/GeoTIFF
Size is 8000, 8000
Coordinate System is:
PROJCS["unnamed",
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"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-115],
PARAMETER["scale_factor",0.9992],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (644613.948217,5517483.066294)
Pixel Size = (20.000000,-20.000000)
Metadata:
TIFFTAG_DOCUMENTNAME=JeffPan8000.tif
TIFFTAG_SOFTWARE=ImageMagick 5.5.2 12/01/02 Q8
http://www.imagemagick.org
Corner Coordinates:
Upper Left ( 644613.948, 5517483.066) (112d59'20.72"W, 49d48'43.97"N)
Lower Left ( 644613.948, 5357483.066) (113d 2'47.06"W, 48d22'23.60"N)
Upper Right ( 804613.948, 5517483.066) (110d46'8.23"W, 49d45'8.47"N)
Lower Right ( 804613.948, 5357483.066) (110d53'21.22"W, 48d18'58.69"N)
Center ( 724613.948, 5437483.066) (111d55'24.30"W, 49d 4'7.55"N)
Band 1 Block=8000x1 Type=Byte, ColorInterp=Gray
C:\matt\gdal-cvs-030403\apps>gdalinfo out.tif
Driver: GTiff/GeoTIFF
Size is 9358, 9357
Coordinate System is `'
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 9357.0)
Upper Right ( 9358.0, 0.0)
Lower Right ( 9358.0, 9357.0)
Center ( 4679.0, 4678.5)
Band 1 Block=9358x1 Type=Byte, ColorInterp=Gray
I have the following code:
#include <stdio.h>
#include "cpl_string.h"
#include "ogr_srs_api.h"
#include "gdalwarper.h"
char *SanitizeSRS( const char *pszUserInput );
static GDALDatasetH GDALWarpCreateOutput( GDALDatasetH hSrcDS, const
char *pszFilename,
const char *pszFormat, const
char *pszSourceSRS,
const char *pszTargetSRS, int
nOrder,
char **papszCreateOptions );
static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double dfXRes=0.0, dfYRes=0.0;
static int nForcePixels=0, nForceLines=0;
int main( int argc, char *argv[] )
{
char *pszFilename = "daImage.tif";
char *pszTargetSRS = "+proj=lcc +lon_0=-100.0 +lat_0=45.0
+lat_1=35.0 +lat_2=55.0 +x_0=0.0 +y_0=0.0 +datum=WGS84";
char *pszSourceSRS = NULL; //"+proj=tmerc +lon_0=-115.0
+lon_0=-115.0 +k=0.9992 +x_0=500000.0 +y_0=0 +datum=NAD83";
int checkValue = NULL;
double adfTargetGeoTransform[6];
int nPixels=0, nLines=0;
void *hTransformArg;
char **papszWarpOptions = NULL;
char **papszCreateOptions = NULL;
CPLErr eErr;
GDALDatasetH hDataset=NULL,hOutDataset=NULL;
GDALDataType eDT;
GDALDriverH hDriver;
GDALWarpOptions *psWarpOptions;
GDALWarpOperation oOperation;
// register raster drivers
GDALAllRegister();
hDataset = GDALOpen( pszFilename, GA_Update );
if( hDataset != NULL )
{
if( pszSourceSRS == NULL )
{
if( GDALGetProjectionRef( hDataset ) != NULL
&& strlen(GDALGetProjectionRef( hDataset )) > 0 )
pszSourceSRS = CPLStrdup(GDALGetProjectionRef( hDataset
));
else if( GDALGetGCPProjection( hDataset ) != NULL
&& strlen(GDALGetGCPProjection(hDataset)) > 0
&& GDALGetGCPCount( hDataset ) > 1 )
pszSourceSRS = CPLStrdup(GDALGetGCPProjection( hDataset
));
else
pszSourceSRS = CPLStrdup("");
}
if( pszTargetSRS == NULL )
pszTargetSRS = CPLStrdup(pszSourceSRS);
//pszSourceSRS = SanitizeSRS( pszSourceSRS );
pszTargetSRS = SanitizeSRS( pszTargetSRS );
hOutDataset = GDALWarpCreateOutput( hDataset,
"out.tif",
"GTiff", // format
pszSourceSRS,
pszTargetSRS,
0, //order
NULL
);
papszWarpOptions = CSLSetNameValue( papszWarpOptions, "INIT",
"0" );
CSLDestroy( papszCreateOptions );
papszCreateOptions = NULL;
if( hOutDataset == NULL )
{
fprintf(stderr,"\nERROR: could not create the output file.");
exit( 1 );
}
else
{
fprintf(stderr,"\nOutput file was successfully created.");
}
// setup and warp
psWarpOptions = GDALCreateWarpOptions();
//psWarpOptions->dfWarpMemoryLimit = 1024.0*1024.0;
psWarpOptions->papszWarpOptions = papszWarpOptions;
psWarpOptions->hSrcDS = hDataset;
psWarpOptions->hDstDS = hOutDataset;
psWarpOptions->nBandCount = 1;
psWarpOptions->panSrcBands = (int *) CPLMalloc(
sizeof(int)*psWarpOptions->nBandCount);
psWarpOptions->panDstBands = (int *) CPLMalloc(
sizeof(int)*psWarpOptions->nBandCount);
psWarpOptions->panSrcBands[0] = 1;
psWarpOptions->panDstBands[0] = 1;
psWarpOptions->pfnProgress = GDALTermProgress;
fprintf(stderr,"\nWarp options are set.");
fprintf(stderr,"\ntransforming from \n%s\n to
\n%s\n",pszSourceSRS,pszTargetSRS);
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer( hDataset,
pszSourceSRS,
hOutDataset, //??
pszTargetSRS, //??
FALSE,
0.0,
1
);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
fprintf(stderr,"\nTransformer is built.");
if( oOperation.Initialize( psWarpOptions ) == CE_None )
{
fprintf(stderr,"\nProgress");
oOperation.ChunkAndWarpImage(0,0,GDALGetRasterXSize(hOutDataset),GDALGet
RasterYSize(hOutDataset) );
}
else
{
fprintf(stderr,"\nERROR: failed to initialize the warp.");
}
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg
);
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDataset );
}
else
{
fprintf( stderr,"\nERROR: could not open %s for
updating.",pszFilename );
}
}
/***********************************************************************
*/
/* SanitizeSRS()
*/
/***********************************************************************
*/
char *SanitizeSRS( const char *pszUserInput )
{
OGRSpatialReferenceH hSRS;
char *pszResult = NULL;
hSRS = OSRNewSpatialReference( NULL );
if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
OSRExportToWkt( hSRS, &pszResult );
OSRDestroySpatialReference( hSRS );
return pszResult;
}
/***********************************************************************
*/
/* GDALWarpCreateOutput()
*/
/*
*/
/* Create the output file based on various commandline options,
*/
/* and the input file.
*/
/***********************************************************************
*/
static GDALDatasetH
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename,
const char *pszFormat, const char *pszSourceSRS,
const char *pszTargetSRS, int nOrder,
char **papszCreateOptions )
{
GDALDriverH hDriver;
GDALDatasetH hDstDS;
void *hTransformArg;
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
GDALColorTableH hCT;
/* --------------------------------------------------------------------
*/
/* Find the output driver.
*/
/* --------------------------------------------------------------------
*/
hDriver = GDALGetDriverByName( pszFormat );
if( hDriver == NULL
|| GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) ==
NULL )
{
int iDr;
printf( "Output driver `%s' not recognised or does not
support\n",
pszFormat );
printf( "direct output file creation. The following format
drivers are configured\n"
"and support direct output:\n" );
for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
{
GDALDriverH hDriver = GDALGetDriver(iDr);
if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) !=
NULL )
{
printf( " %s: %s\n",
GDALGetDriverShortName( hDriver ),
GDALGetDriverLongName( hDriver ) );
}
}
printf( "\n" );
exit( 1 );
}
/* --------------------------------------------------------------------
*/
/* Create a transformation object from the source to
*/
/* destination coordinate system.
*/
/* --------------------------------------------------------------------
*/
hTransformArg =
GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS,
NULL, pszTargetSRS,
TRUE, 1000.0, nOrder );
if( hTransformArg == NULL )
return NULL;
/* --------------------------------------------------------------------
*/
/* Get approximate output definition.
*/
/* --------------------------------------------------------------------
*/
if( GDALSuggestedWarpOutput( hSrcDS,
GDALGenImgProjTransform, hTransformArg,
adfDstGeoTransform, &nPixels, &nLines )
!= CE_None )
return NULL;
GDALDestroyGenImgProjTransformer( hTransformArg );
/* --------------------------------------------------------------------
*/
/* Did the user override some parameters?
*/
/* --------------------------------------------------------------------
*/
if( dfXRes != 0.0 && dfYRes != 0.0 )
{
CPLAssert( nPixels == 0 && nLines == 0 );
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY ==
0.0 )
{
dfMinX = adfDstGeoTransform[0];
dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] *
nPixels;
dfMaxY = adfDstGeoTransform[3];
dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] *
nLines;
}
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
}
else if( nForcePixels != 0 && nForceLines != 0 )
{
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY ==
0.0 )
{
dfMinX = adfDstGeoTransform[0];
dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] *
nPixels;
dfMaxY = adfDstGeoTransform[3];
dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] *
nLines;
}
dfXRes = (dfMaxX - dfMinX) / nForcePixels;
dfYRes = (dfMaxY - dfMinY) / nForceLines;
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
adfDstGeoTransform[1] = dfXRes;
adfDstGeoTransform[5] = -dfYRes;
nPixels = nForcePixels;
nLines = nForceLines;
}
else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY
!= 0.0 )
{
dfXRes = adfDstGeoTransform[1];
dfYRes = fabs(adfDstGeoTransform[5]);
nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
adfDstGeoTransform[0] = dfMinX;
adfDstGeoTransform[3] = dfMaxY;
}
/* --------------------------------------------------------------------
*/
/* Create the output file.
*/
/* --------------------------------------------------------------------
*/
printf( "Creating output file is that %dP x %dL.\n", nPixels, nLines
);
hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
GDALGetRasterCount(hSrcDS),
GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)),
papszCreateOptions );
if( hDstDS == NULL )
return NULL;
/* --------------------------------------------------------------------
*/
/* Write out the projection definition.
*/
/* --------------------------------------------------------------------
*/
GDALSetProjection( hDstDS, pszTargetSRS );
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
/* --------------------------------------------------------------------
*/
/* Copy the color table, if required.
*/
/* --------------------------------------------------------------------
*/
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
if( hCT != NULL )
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
return hDstDS;
}
---
Outgoing mail is certified Virus Free.
Checked by AVG anti-virus system (http://www.grisoft.com).
Version: 6.0.507 / Virus Database: 304 - Release Date: 8/4/2003
More information about the Gdal-dev
mailing list