EVOLUTION-MANAGER
Edit File: warptut.html
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> <html xmlns="http://www.w3.org/1999/xhtml"> <head> <meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/> <meta http-equiv="X-UA-Compatible" content="IE=9"/> <meta name="generator" content="Doxygen 1.8.5"/> <title>My Project: GDAL Warp API Tutorial</title> <link href="tabs.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="jquery.js"></script> <script type="text/javascript" src="dynsections.js"></script> <link href="search/search.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="search/search.js"></script> <script type="text/javascript"> $(document).ready(function() { searchBox.OnSelectItem(0); }); </script> <link href="doxygen.css" rel="stylesheet" type="text/css" /> </head> <body> <div id="top"><!-- do not remove this div, it is closed by doxygen! --> <div id="titlearea"> <table cellspacing="0" cellpadding="0"> <tbody> <tr style="height: 56px;"> <td style="padding-left: 0.5em;"> <div id="projectname">My Project </div> </td> </tr> </tbody> </table> </div> <!-- end header part --> <!-- Generated by Doxygen 1.8.5 --> <script type="text/javascript"> var searchBox = new SearchBox("searchBox", "search",false,'Search'); </script> <div id="navrow1" class="tabs"> <ul class="tablist"> <li><a href="index.html"><span>Main Page</span></a></li> <li class="current"><a href="pages.html"><span>Related Pages</span></a></li> <li> <div id="MSearchBox" class="MSearchBoxInactive"> <span class="left"> <img id="MSearchSelect" src="search/mag_sel.png" onmouseover="return searchBox.OnSearchSelectShow()" onmouseout="return searchBox.OnSearchSelectHide()" alt=""/> <input type="text" id="MSearchField" value="Search" accesskey="S" onfocus="searchBox.OnSearchFieldFocus(true)" onblur="searchBox.OnSearchFieldFocus(false)" onkeyup="searchBox.OnSearchFieldChange(event)"/> </span><span class="right"> <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a> </span> </div> </li> </ul> </div> <!-- window showing the filter options --> <div id="MSearchSelectWindow" onmouseover="return searchBox.OnSearchSelectShow()" onmouseout="return searchBox.OnSearchSelectHide()" onkeydown="return searchBox.OnSearchSelectKey(event)"> <a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark"> </span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark"> </span>Pages</a></div> <!-- iframe showing the search results (closed by default) --> <div id="MSearchResultsWindow"> <iframe src="javascript:void(0)" frameborder="0" name="MSearchResults" id="MSearchResults"> </iframe> </div> </div><!-- top --> <div class="header"> <div class="headertitle"> <div class="title">GDAL Warp API Tutorial </div> </div> </div><!--header--> <div class="contents"> <div class="textblock"><h1><a class="anchor" id="warptut_overview"></a> Overview</h1> <p>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.</p> <p>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 <a href="gdal_datamodel.html">GDAL Data Model</a>, and the general GDAL API.</p> <p>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.</p> <h1><a class="anchor" id="warptut_simple"></a> A Simple Reprojection Case</h1> <p>First we will construct a relatively simple example for reprojecting an image, assuming an appropriate output file already exists, and with minimal error checking.</p> <div class="fragment"><div class="line"><span class="preprocessor">#include "gdalwarper.h"</span></div> <div class="line"></div> <div class="line"><span class="keywordtype">int</span> main()</div> <div class="line">{</div> <div class="line"> GDALDatasetH hSrcDS, hDstDS;</div> <div class="line"></div> <div class="line"> <span class="comment">// Open input and output files. </span></div> <div class="line"></div> <div class="line"> GDALAllRegister();</div> <div class="line"></div> <div class="line"> hSrcDS = GDALOpen( <span class="stringliteral">"in.tif"</span>, GA_ReadOnly );</div> <div class="line"> hDstDS = GDALOpen( <span class="stringliteral">"out.tif"</span>, GA_Update );</div> <div class="line"></div> <div class="line"> <span class="comment">// Setup warp options. </span></div> <div class="line"> </div> <div class="line"> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();</div> <div class="line"></div> <div class="line"> psWarpOptions->hSrcDS = hSrcDS;</div> <div class="line"> psWarpOptions->hDstDS = hDstDS;</div> <div class="line"></div> <div class="line"> psWarpOptions->nBandCount = 1;</div> <div class="line"> psWarpOptions->panSrcBands = </div> <div class="line"> (<span class="keywordtype">int</span> *) CPLMalloc(<span class="keyword">sizeof</span>(<span class="keywordtype">int</span>) * psWarpOptions->nBandCount );</div> <div class="line"> psWarpOptions->panSrcBands[0] = 1;</div> <div class="line"> psWarpOptions->panDstBands = </div> <div class="line"> (<span class="keywordtype">int</span> *) CPLMalloc(<span class="keyword">sizeof</span>(<span class="keywordtype">int</span>) * psWarpOptions->nBandCount );</div> <div class="line"> psWarpOptions->panDstBands[0] = 1;</div> <div class="line"></div> <div class="line"> psWarpOptions->pfnProgress = GDALTermProgress; </div> <div class="line"></div> <div class="line"> <span class="comment">// Establish reprojection transformer. </span></div> <div class="line"></div> <div class="line"> psWarpOptions->pTransformerArg = </div> <div class="line"> GDALCreateGenImgProjTransformer( hSrcDS, </div> <div class="line"> GDALGetProjectionRef(hSrcDS), </div> <div class="line"> hDstDS,</div> <div class="line"> GDALGetProjectionRef(hDstDS), </div> <div class="line"> FALSE, 0.0, 1 );</div> <div class="line"> psWarpOptions->pfnTransformer = GDALGenImgProjTransform;</div> <div class="line"></div> <div class="line"> <span class="comment">// Initialize and execute the warp operation. </span></div> <div class="line"></div> <div class="line"> GDALWarpOperation oOperation;</div> <div class="line"></div> <div class="line"> oOperation.Initialize( psWarpOptions );</div> <div class="line"> oOperation.ChunkAndWarpImage( 0, 0, </div> <div class="line"> GDALGetRasterXSize( hDstDS ), </div> <div class="line"> GDALGetRasterYSize( hDstDS ) );</div> <div class="line"></div> <div class="line"> GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );</div> <div class="line"> GDALDestroyWarpOptions( psWarpOptions );</div> <div class="line"></div> <div class="line"> GDALClose( hDstDS );</div> <div class="line"> GDALClose( hSrcDS );</div> <div class="line"></div> <div class="line"> <span class="keywordflow">return</span> 0;</div> <div class="line">}</div> </div><!-- fragment --><p>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.</p> <p>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.</p> <p>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.</p> <p>Normally error check would be needed after opening files, setting up the reprojection transformer (returns NULL on failure), and initializing the warp.</p> <h1><a class="anchor" id="warptut_options"></a> Other Warping Options</h1> <p>The GDALWarpOptions structures contains a number of items that can be set to control warping behavior. A few of particular interest are:</p> <ol> <li> <p class="startli">GDALWarpOptions::dfWarpMemoryLimit - Set the maximum amount of memory to be used by the GDALWarpOperation when selecting a size of image chunk to operate on. The value is in bytes, and the default is likely to be conservative (small). Increasing the chunk size can help substantially in some situations but care should be taken to ensure that this size, plus the GDAL cache size plus the working set of GDAL, your application and the operating system are less than the size of RAM or else excessive swapping is likely to interfere with performance. On a system with 256MB of RAM, a value of at least 64MB (roughly 64000000 bytes) is reasonable. Note that this value does <b>not</b> include the memory used by GDAL for low level block caching.</p> <p class="endli"></p> </li> <li> <p class="startli">GDALWarpOpations::eResampleAlg - One of GRA_NearestNeighbour (the default, and fastest), GRA_Bilinear (2x2 bilinear resampling) or GRA_Cubic. The GRA_NearestNeighbour type should generally be used for thematic or colormapped images. The other resampling types may give better results for thematic images, especially when substantially changing resolution.</p> <p class="endli"></p> </li> <li> <p class="startli">GDALWarpOptions::padfSrcNoDataReal - This array (one entry per band being processed) may be setup with a "nodata" value for each band if you wish to avoid having pixels of some background value copied to the destination image.</p> <p class="endli"></p> </li> <li> <p class="startli"><a class="anchor" id="#warpoptions"></a> GDALWarpOptions::papszWarpOptions - This is a string list of NAME=VALUE options passed to the warper. See the GDALWarpOptions::papszWarpOptions docs for all options. Supported values include:</p> <ul> <li> <p class="startli">INIT_DEST=[value] or INIT_DEST=NO_DATA: This option forces the destination image to be initialized to the indicated value (for all bands) or indicates that it should be initialized to the NO_DATA value in padfDstNoDataReal/padfDstNoDataImag. If this value isn't set the destination image will be read and the source warp overlayed on it.</p> <p class="endli"></p> </li> <li> WRITE_FLUSH=YES/NO: This option forces a flush to disk of data after each chunk is processed. In some cases this helps ensure a serial writing of the output data otherwise a block of data may be written to disk each time a block of data is read for the input buffer resulting in a lot of extra seeking around the disk, and reduced IO throughput. The default at this time is NO. </li> </ul> <p class="endli"></p> </li> </ol> <h1><a class="anchor" id="warptut_output"></a> Creating the Output File</h1> <p>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.</p> <div class="fragment"><div class="line"><span class="preprocessor">#include "gdalwarper.h"</span></div> <div class="line"><span class="preprocessor">#include "ogr_spatialref.h"</span></div> <div class="line"></div> <div class="line">...</div> <div class="line"></div> <div class="line"> GDALDriverH hDriver;</div> <div class="line"> GDALDataType eDT;</div> <div class="line"> GDALDatasetH hDstDS;</div> <div class="line"> GDALDatasetH hSrcDS;</div> <div class="line"></div> <div class="line"> <span class="comment">// Open the source file. </span></div> <div class="line"></div> <div class="line"> hSrcDS = GDALOpen( <span class="stringliteral">"in.tif"</span>, GA_ReadOnly );</div> <div class="line"> CPLAssert( hSrcDS != NULL );</div> <div class="line"> </div> <div class="line"> <span class="comment">// Create output with same datatype as first input band. </span></div> <div class="line"></div> <div class="line"> eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));</div> <div class="line"></div> <div class="line"> <span class="comment">// Get output driver (GeoTIFF format)</span></div> <div class="line"></div> <div class="line"> hDriver = GDALGetDriverByName( <span class="stringliteral">"GTiff"</span> );</div> <div class="line"> CPLAssert( hDriver != NULL );</div> <div class="line"></div> <div class="line"> <span class="comment">// Get Source coordinate system. </span></div> <div class="line"></div> <div class="line"> <span class="keyword">const</span> <span class="keywordtype">char</span> *pszSrcWKT, *pszDstWKT = NULL;</div> <div class="line"></div> <div class="line"> pszSrcWKT = GDALGetProjectionRef( hSrcDS );</div> <div class="line"> CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );</div> <div class="line"></div> <div class="line"> <span class="comment">// Setup output coordinate system that is UTM 11 WGS84. </span></div> <div class="line"></div> <div class="line"> OGRSpatialReference oSRS;</div> <div class="line"></div> <div class="line"> oSRS.SetUTM( 11, TRUE );</div> <div class="line"> oSRS.SetWellKnownGeogCS( <span class="stringliteral">"WGS84"</span> );</div> <div class="line"></div> <div class="line"> oSRS.exportToWkt( &pszDstWKT );</div> <div class="line"></div> <div class="line"> <span class="comment">// Create a transformer that maps from source pixel/line coordinates</span></div> <div class="line"> <span class="comment">// to destination georeferenced coordinates (not destination </span></div> <div class="line"> <span class="comment">// pixel line). We do that by omitting the destination dataset</span></div> <div class="line"> <span class="comment">// handle (setting it to NULL). </span></div> <div class="line"></div> <div class="line"> <span class="keywordtype">void</span> *hTransformArg;</div> <div class="line"></div> <div class="line"> hTransformArg = </div> <div class="line"> GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, </div> <div class="line"> FALSE, 0, 1 );</div> <div class="line"> CPLAssert( hTransformArg != NULL );</div> <div class="line"></div> <div class="line"> <span class="comment">// Get approximate output georeferenced bounds and resolution for file. </span></div> <div class="line"></div> <div class="line"> <span class="keywordtype">double</span> adfDstGeoTransform[6];</div> <div class="line"> <span class="keywordtype">int</span> nPixels=0, nLines=0;</div> <div class="line"> CPLErr eErr;</div> <div class="line"></div> <div class="line"> eErr = GDALSuggestedWarpOutput( hSrcDS, </div> <div class="line"> GDALGenImgProjTransform, hTransformArg, </div> <div class="line"> adfDstGeoTransform, &nPixels, &nLines );</div> <div class="line"> CPLAssert( eErr == CE_None );</div> <div class="line"></div> <div class="line"> GDALDestroyGenImgProjTransformer( hTransformArg );</div> <div class="line"></div> <div class="line"> <span class="comment">// Create the output file. </span></div> <div class="line"></div> <div class="line"> hDstDS = GDALCreate( hDriver, <span class="stringliteral">"out.tif"</span>, nPixels, nLines, </div> <div class="line"> GDALGetRasterCount(hSrcDS), eDT, NULL );</div> <div class="line"> </div> <div class="line"> CPLAssert( hDstDS != NULL );</div> <div class="line"></div> <div class="line"> <span class="comment">// Write out the projection definition. </span></div> <div class="line"></div> <div class="line"> GDALSetProjection( hDstDS, pszDstWKT );</div> <div class="line"> GDALSetGeoTransform( hDstDS, adfDstGeoTransform );</div> <div class="line"></div> <div class="line"> <span class="comment">// Copy the color table, if required.</span></div> <div class="line"></div> <div class="line"> GDALColorTableH hCT;</div> <div class="line"></div> <div class="line"> hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );</div> <div class="line"> <span class="keywordflow">if</span>( hCT != NULL )</div> <div class="line"> GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );</div> <div class="line"></div> <div class="line"> ... proceed with warp as before ...</div> </div><!-- fragment --><p>Some notes on this logic:</p> <ul> <li> <p class="startli">We need to create the transformer to output coordinates such that the output of the transformer is georeferenced, not pixel line coordinates since we use the transformer to map pixels around the source image into destination georeferenced coordinates.</p> <p class="endli"></p> </li> <li> <p class="startli">The GDALSuggestedWarpOutput() function will return an adfDstGeoTransform, nPixels and nLines that describes an output image size and georeferenced extents that should hold all pixels from the source image. The resolution is intended to be comparable to the source, but the output pixels are always square regardless of the shape of input pixels.</p> <p class="endli"></p> </li> <li> <p class="startli">The warper requires an output file in a format that can be "randomly" written to. This generally limits things to uncompressed formats that have an implementation of the Create() method (as opposed to CreateCopy()). To warp to compressed formats, or CreateCopy() style formats it is necessary to produce a full temporary copy of the image in a better behaved format, and then CreateCopy() it to the desired final format.</p> <p class="endli"></p> </li> <li> <p class="startli">The Warp API copies only pixels. All colormaps, georeferencing and other metadata must be copied to the destination by the application.</p> <p class="endli"></p> </li> </ul> <h1><a class="anchor" id="warptut_perfomance"></a> Performance Optimization</h1> <p>There are a number of things that can be done to optimize the performance of the warp API.</p> <ol> <li> <p class="startli">Increase the amount of memory available for the Warp API chunking so that larger chunks can be operated on at a time. This is the GDALWarpOptions::dfWarpMemoryLimit parameter. In theory the larger the chunk size operated on the more efficient the I/O strategy, and the more efficient the approximated transformation will be. However, the sum of the warp memory and the GDAL cache should be less than RAM size, likely around 2/3 of RAM size.</p> <p class="endli"></p> </li> <li> <p class="startli">Increase the amount of memory for GDAL caching. This is especially important when working with very large input and output images that are scanline oriented. If all the input or output scanlines have to be re-read for each chunk they intersect performance may degrade greatly. Use GDALSetCacheMax() to control the amount of memory available for caching within GDAL.</p> <p class="endli"></p> </li> <li> <p class="startli">Use an approximated transformation instead of exact reprojection for each pixel to be transformed. This code illustrates how an approximated transformation could be created based on a reprojection transformation, but with a given error threshold (dfErrorThreshold in output pixels).</p> <div class="fragment"><div class="line">hTransformArg = </div> <div class="line"> GDALCreateApproxTransformer( GDALGenImgProjTransform, </div> <div class="line"> hGenImgProjArg, dfErrorThreshold );</div> <div class="line">pfnTransformer = GDALApproxTransform;</div> </div><!-- fragment --><p class="endli"></p> </li> <li> <p class="startli">When writing to a blank output file, use the INIT_DEST option in the GDALWarpOptions::papszWarpOptions to cause the output chunks to be initialized to a fixed value, instead of being read from the output. This can substantially reduce unnecessary IO work.</p> <p class="endli"></p> </li> <li> <p class="startli">Use tiled input and output formats. Tiled formats allow a given chunk of source and destination imagery to be accessed without having to touch a great deal of extra image data. Large scanline oriented files can result in a great deal of wasted extra IO.</p> <p class="endli"></p> </li> <li> <p class="startli">Process all bands in one call. This ensures the transformation calculations don't have to be performed for each band.</p> <p class="endli"></p> </li> <li> <p class="startli">Use the GDALWarpOperation::ChunkAndWarpMulti() method instead of GDALWarpOperation::ChunkAndWarpImage(). It uses a separate thread for the IO and the actual image warp operation allowing more effective use of CPU and IO bandwidth. For this to work GDAL needs to have been built with multi-threading support (default on Win32, default on Unix since GDAL 1.8.0, for previous versions –with-threads was required in configure).</p> <p class="endli"></p> </li> <li> <p class="startli">The resampling kernels vary is work required from nearest neighbour being least, then bilinear then cubic. Don't use a more complex resampling kernel than needed.</p> <p class="endli"></p> </li> <li> <p class="startli">Avoid use of esoteric masking options so that special simplified logic case be used for common special cases. For instance, nearest neighbour resampling with no masking on 8bit data is highly optimized compared to the general case.</p> <p class="endli"></p> </li> </ol> <h1><a class="anchor" id="warptut_other_opts"></a> Other Masking Options</h1> <p>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.</p> <p> $Id: warptut.dox 21106 2010-11-09 19:19:21Z rouault $ </p> </div></div><!-- contents --> <!-- start footer part --> <hr class="footer"/><address class="footer"><small> Generated by  <a href="http://www.doxygen.org/index.html"> <img class="footer" src="doxygen.png" alt="doxygen"/> </a> 1.8.5 </small></address> </body> </html>