[Gdal-dev] Extracting subimage in GDAL

[Frank Warmerdam]

Ozy Sjahputera wrote:
> Greetings,
> 
> I am new to GDAL so I apologize if this topic has been discussed 
> before.  I used Google to search for several key words -- like 
> "subimage" -- in the list archives.  The only one I found to be close is 
> this: http://lists.maptools.org/pipermail/gdal-dev/2004-May/002776.html 
> where Frank pointed out the use of GDAL utility gdal_translate() to 
> extract sub-images.
> 
> I need to extract subimage by specifying the new UL and LR coordinates 
> (in UTM).  I understand that gdal_translate with -projwin option can do 
> this.  However, I need to do this as a function call within a CPP 
> program.  Is there a function in GDAL library that is equivalent to the 
> utility gdal_translate?  

Ozy,

No.  But if you know the source window in pixel/line coordinates you
can load it into a buffer with the GDALRasterBand::RasterIO() method
or the GDALDataset::RasterIO() method.

> Or where can I find the source code for 
> gdal_translate if I want to learn more about how gdal_translate is 
> implemented?

The gdal_translate source can be found in gdal/apps/gdal_translate.c
or on the web:
   http://www.gdal.org/srctree/apps/gdal_translate.c

The pertinent portion for -projwin is:

/* -------------------------------------------------------------------- */
/*      Compute the source window from the projected source window      */
/*      if the projected coordinates were provided.  Note that the      */
/*      projected coordinates are in ulx, uly, lrx, lry format,         */
/*      while the anSrcWin is xoff, yoff, xsize, ysize with the         */
/*      xoff,yoff being the ulx, uly in pixel/line.                     */
/* -------------------------------------------------------------------- */
     if( dfULX != 0.0 || dfULY != 0.0
         || dfLRX != 0.0 || dfLRY != 0.0 )
     {
         double	adfGeoTransform[6];

         GDALGetGeoTransform( hDataset, adfGeoTransform );

         if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 )
         {
             fprintf( stderr,
                      "The -projwin option was used, but the geotransform is\n"
                      "rotated.  This configuration is not supported.\n" );
             GDALClose( hDataset );
             CPLFree( panBandList );
             GDALDestroyDriverManager();
             exit( 1 );
         }

         anSrcWin[0] = (int)
             ((dfULX - adfGeoTransform[0]) / adfGeoTransform[1] + 0.001);
         anSrcWin[1] = (int)
             ((dfULY - adfGeoTransform[3]) / adfGeoTransform[5] + 0.001);

         anSrcWin[2] = (int) ((dfLRX - dfULX) / adfGeoTransform[1] + 0.5);
         anSrcWin[3] = (int) ((dfLRY - dfULY) / adfGeoTransform[5] + 0.5);

         if( !bQuiet )
             fprintf( stdout,
                      "Computed -srcwin %d %d %d %d from projected window.\n",
                      anSrcWin[0],
                      anSrcWin[1],
                      anSrcWin[2],
                      anSrcWin[3] );

         if( anSrcWin[0] < 0 || anSrcWin[1] < 0
             || anSrcWin[0] + anSrcWin[2] > GDALGetRasterXSize(hDataset)
             || anSrcWin[1] + anSrcWin[3] > GDALGetRasterYSize(hDataset) )
         {
             fprintf( stderr,
                      "Computed -srcwin falls outside raster size of %dx%d.\n",
                      GDALGetRasterXSize(hDataset),
                      GDALGetRasterYSize(hDataset) );
             exit( 1 );
         }
     }


Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org