[gdal-dev] How to implement tile read / write to gdal supported format file?

[Luke]

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.