The GDAL Warp API (declared in gdalwarper.h) provides services for high performance image warping using application provided geometric transformation functions (GDALTransformerFunc), a variety of resampling kernels, and various masking options. Files much larger than can be held in memory can be warped.
This tutorial demonstrates how to implement an application using the Warp API. It assumes implementation in C++ as C and Python bindings are incomplete for the Warp API. It also assumes familiarity with the GDAL Data Model, and the general GDAL API.
Applications normally perform a warp by initializing a GDALWarpOptions structure with the options to be utilized, instantiating a GDALWarpOperation based on these options, and then invoking the GDALWarpOperation::ChunkAndWarpImage() method to perform the warp options internally using the GDALWarpKernel class.
First we will construct a relatively simple example for reprojecting an image, assuming an appropriate output file already exists, and with minimal error checking.
#include "gdalwarper.h" int main() { GDALDatasetH hSrcDS, hDstDS; // Open input and output files. GDALAllRegister(); hSrcDS = GDALOpen( "in.tif", GA_ReadOnly ); hDstDS = GDALOpen( "out.tif", GA_Update ); // Setup warp options. GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->hSrcDS = hSrcDS; psWarpOptions->hDstDS = hDstDS; psWarpOptions->nBandCount = 1; psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); psWarpOptions->panSrcBands[0] = 1; psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount ); psWarpOptions->panDstBands[0] = 1; psWarpOptions->pfnProgress = GDALTermProgress; // Establish reprojection transformer. psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 0.0, 1 ); psWarpOptions->pfnTransformer = GDALGenImgProjTransform; // Initialize and execute the warp operation. GDALWarpOperation oOperation; oOperation.Initialize( psWarpOptions ); oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ); GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg ); GDALDestroyWarpOptions( psWarpOptions ); GDALClose( hDstDS ); GDALClose( hSrcDS ); return 0; }
This example opens the existing input and output files (in.tif and out.tif). A GDALWarpOptions structure is allocated (GDALCreateWarpOptions() sets lots of sensible defaults for stuff, always use it for defaulting things), and the input and output file handles, and band lists are set. The panSrcBands and panDstBands lists are dynamically allocated here and will be free automatically by GDALDestroyWarpOptions(). The simple terminal output progress monitor (GDALTermProgress) is installed for reporting completion progress to the user.
GDALCreateGenImgProjTransformer() is used to initialize the reprojection transformation between the source and destination images. We assume that they already have reasonable bounds and coordinate systems set. Use of GCPs is disabled.
Once the options structure is ready, a GDALWarpOperation is instantiated using them, and the warp actually performed with GDALWarpOperation::ChunkAndWarpImage(). Then the transformer, warp options and datasets are cleaned up.
Normally error check would be needed after opening files, setting up the reprojection transformer (returns NULL on failure), and initializing the warp.
The GDALWarpOptions structures contains a number of items that can be set to control warping behaviour. A few of particular interest are:
In the previous case an appropriate output file was already assumed to exist. Now we will go through a case where a new file with appropriate bounds in a new coordinate system is created. This operation doesn't relate specifically to the warp API. It is just using the transformation API.
#include "gdalwarper.h" #include "ogr_spatialref.h" ... GDALDriverH hDriver; GDALDataType eDT; GDALDatasetH hDstDS; GDALDatasetH hSrcDS; // Open the source file. hSrcDS = GDALOpen( "in.tif", GA_ReadOnly ); CPLAssert( hSrcDS != NULL ); // Create output with same datatype as first input band. eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)); // Get output driver (GeoTIFF format) hDriver = GDALGetDriverByName( "GTiff" ); CPLAssert( hDriver != NULL ); // Get Source coordinate system. const char *pszSrcWKT, *pszDstWKT = NULL; pszSrcWKT = GDALGetProjectionRef( hSrcDS ); CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 ); // Setup output coordinate system that is UTM 11 WGS84. OGRSpatialReference oSRS; oSRS.SetUTM( 11, TRUE ); oSRS.SetWellKnownGeogCS( "WGS84" ); oSRS.exportToWkt( &pszDstWKT ); // Create a transformer that maps from source pixel/line coordinates // to destination georeferenced coordinates (not destination // pixel line). We do that by omitting the destination dataset // handle (setting it to NULL). void *hTransformArg; hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 ); CPLAssert( hTransformArg != NULL ); // Get approximate output georeferenced bounds and resolution for file. double adfDstGeoTransform[6]; int nPixels=0, nLines=0; CPLErr eErr; eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines ); CPLAssert( eErr == CE_None ); GDALDestroyGenImgProjTransformer( hTransformArg ); // Create the output file. hDstDS = GDALCreate( hDriver, "out.tif", nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, NULL ); CPLAssert( hDstDS != NULL ); // Write out the projection definition. GDALSetProjection( hDstDS, pszDstWKT ); GDALSetGeoTransform( hDstDS, adfDstGeoTransform ); // Copy the color table, if required. GDALColorTableH hCT; hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) ); if( hCT != NULL ) GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT ); ... proceed with warp as before ...
Some notes on this logic:
There are a number of things that can be done to optimize the performance of the warp API.
hTransformArg = GDALCreateApproxTransformer( GDALGenImgProjTransform, hGenImgProjArg, dfErrorThreshold ); pfnTransformer = GDALApproxTransform;
The GDALWarpOptions include a bunch of esoteric masking capabilities, for validity masks, and density masks on input and output. Some of these are not yet implemented and others are implemented but poorly tested. Other than per-band validity masks it is advised that these features be used with caution at this time.