[gdal-dev] gdalwarp rather than gdal_translate

Frank Warmerdam warmerdam at pobox.com
Fri Feb 15 13:38:02 EST 2008


Diana Esch-Mosher wrote:
> Frank,
> 
> So I have I done something incorrectly?
> 
>                    ds=gdal.Open(image) # image is a BIG nitf
>                    band=ds.GetRasterBand(1)
>                    sizeX= ds.RasterXSize
>                    sizeY= ds.RasterYSize
>                    driver=gdal.GetDriverByName('GTiff')
>                    sizeXNew = int(math.ceil(factor*sizeX))
>                    sizeYNew = int(math.ceil(factor*sizeY))
>                    resampledFile=driver.Create(outFilename, sizeXNew, 
> sizeYNew, bands=1, datatype=band.DataType)
>                    resampledFile.SetProjection(ds.GetProjection())
>                    resampledFile.SetGeoTransform(ds.GetGeoTransform())
>                    resampledFile.SetMetadata(ds.GetMetadata())
>                    cmd = string.join(["gdalwarp -rc" ,image,outFilename])
>                    p = subprocess.Popen(cmd, 
> shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
>                    (stdout, stderr) = p.communicate(input=None)

Diana,

You need to adjust the pixel size in the geotransform to the same degree
that the image size in pixels and lines is being adjusted.


This code in gdal_translate.cpp accomplishes it in that context, though
it is more complicated than you will need since it also accounts for
subwindowing.

     else if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None
         && nGCPCount == 0 )
     {
         adfGeoTransform[0] += anSrcWin[0] * adfGeoTransform[1]
             + anSrcWin[1] * adfGeoTransform[2];
         adfGeoTransform[3] += anSrcWin[0] * adfGeoTransform[4]
             + anSrcWin[1] * adfGeoTransform[5];

         adfGeoTransform[1] *= anSrcWin[2] / (double) nOXSize;
         adfGeoTransform[2] *= anSrcWin[3] / (double) nOYSize;
         adfGeoTransform[4] *= anSrcWin[2] / (double) nOXSize;
         adfGeoTransform[5] *= anSrcWin[3] / (double) nOYSize;

         poVDS->SetGeoTransform( adfGeoTransform );
     }

Keep in mind that the geotransform is essentially an origin and pixel
size - not extents.

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



More information about the gdal-dev mailing list