You could do something like the following: (code modified from the gdal_calculations library - http://code.google.com/p/gdal-calculations)
from osgeo import gdal class Block(object): def __init__(self, dataset, band, xoff, yoff, xsize, ysize): self.xoff = xoff self.yoff = yoff self.xsize = xsize self.ysize = ysize if not band: self.data = dataset.ReadAsArray(xoff, yoff, xsize, ysize) else: self.data = dataset.GetRasterBand(band).ReadAsArray(xoff, yoff, xsize, ysize) def read_blocks_as_array(dataset, band=0, nblocks=1): '''Read GDAL Datasets or individual Bands block by block''' ncols=dataset.RasterXSize() nrows=dataset.RasterYSize() xblock,yblock=dataset.GetRasterBand(1).GetBlockSize() if xblock==ncols: yblock*=nblocks else: xblock*=nblocks for yoff in xrange(0, nrows, yblock): if yoff + yblock < nrows: ysize = yblock else: ysize = nrows - yoff for xoff in xrange(0, ncols, xblock): if xoff + xblock < ncols: xsize = xblock else: xsize = ncols - xoff yield Block(dataset, band, xoff, yoff, xsize, ysize) ds1=gdal.Open(some_raster) ds2=gdal.Open(another_raster) out_driver=gdal.GetDriverByName('GTIFF') out_ds=out_driver.Create (out_raster,ds1.RasterXSize(),ds1.RasterYSize(), 1,ds1.GetRasterBand(1).DataType, ['BIGTIFF=IF_SAFER']) out_ds.SetGeoTransform(ds1.GetGeoTransform()) out_ds.SetProjection(ds1.GetProjection()) out_band=out_ds.GetRasterBand(1) for block in read_blocks_as_array(ds1,band=1,nblocks=10): block2=Block(ds2, 1, block.xoff, block.yoff, block.xsize, block.ysize) ndarray1=block.data ndarray2=block2.data result = ndarray1+ndarray2 out_band.WriteArray(result, block.xoff, block.yoff) -- View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-How-to-implement-tile-read-write-to-gdal-supported-format-file-tp5089581p5089803.html Sent from the GDAL - Dev mailing list archive at Nabble.com. _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev