GDAL and Python: How to "pad the extent" or otherwise handle band.ReadRaster() for area falling slightly outside the image bounds?

Question

In a nutshell, I'm trying to create a bare-metal tiling utility using GDAL and Python. As you can see in the graphic below, which shows the perimeter of a raster and the OSM-structured tile index, some tiles overlap the raster's native extent. This prevents me from calling band.ReadRaster() on the image, as doing so throws an exception.

So I'm looking for a way to handle this. I have a couple ideas, but I'm not sure how to implement them; I'll go ahead and list them:

1) Somehow "Pad" the raster extent with alpha/NODATA values before running my script. Usually people want to crop collars, so I'm not real sure how to add a NODATA collar. And I'm not fond of this approach for two big reasons. First, I would need to remember to pad every image before tiling it this way. Second, each tile index will have a different footprint, meaning each iteration would need its own padded source raster. If I choose this route, though, should I just make a huge NODATA VRT for the full extent of the tile index, then gdal_merge it with my source raster? That sounds somewhat efficient, as I could just output another VRT and avoid waiting on a big render. But how would I actually make the NODATA VRT?? Calling ogrinfo would provide the extents, but then what?

2) Somehow pad on-the-fly in the script/runtime. This approach seems more elegant, but how can I handle band.ReadRaster() on an area that falls outside the image extent? Any tips/tricks?

Alternatively, is there something even better I haven't thought of?

Here's a graphic illustrating the problem. The tile index is in red, and you can see where the perimeter tiles extend beyond the source raster.

Answer 1

I'm not sure that this answers the OP's question, but I found this after searching for how to pad a geotiff with gdal. I am writing this answer because I tried the method of @ProudGIS and I was getting strange behavior. The output tif or vrt using -projwin was 8 pixels off in the y-direction and 12 pixels off in the x-direction after padding. All the metadata was correct for the padded tiff (i.e. origin, resolution, and nrows/ncols). I couldn't figure out the problem, so I wrote some code that pads a numpy array and writes it to a new tiff/vrt.

import numpy as np
import gdal
import os

def pad_geotiff(pathtopad, npad, outpath=0, padval=0):
    """
    Pads a geotiff image by adding npad pixels to each edge.
    """
    base, folder, file, ext = gh.parse_path(pathtopad)
    if outpath == 0:
        outpath = os.path.join(base,folder,file) + '_pad' + ext

    topad = gdal.Open(pathtopad)
    gt = topad.GetGeoTransform()
    colortable = topad.GetRasterBand(1).GetColorTable()
    data_type = topad.GetRasterBand(1).DataType
    Itopad = topad.ReadAsArray()

    ulx = gt[0] - gt[1] * npad
    uly = gt[3] - gt[5] * npad
#    lrx = gt[0] + gt[1] * (topad.RasterXSize + npad)
#    lry = gt[3] + gt[5] * (topad.RasterYSize + npad)

    # Make new geotransform
    gt_new = (ulx, gt[1], gt[2], uly, gt[4], gt[5])

    # Make padded raster (pad with zeros)
    raster = np.pad(Itopad, npad, mode='constant', constant_values=padval)

    write_tile(raster, gt_new, topad, outpath, dtype=data_type, color_table=colortable)

    return outpath


def write_tile(raster, gt, data_obj, outputpath, dtype=gdal.GDT_UInt16, options=0, color_table=0, nbands=1, nodata=False):

    width = np.shape(raster)[1]
    height = np.shape(raster)[0]

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(raster)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    wkt = data_obj.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())


    # Close output raster dataset 
    dest = None

def parse_path(path):

    """
    Parses a file or folderpath into: base, folder (where folder is the 
    outermost subdirectory), filename, and extention. Filename and extension
    are empty if a directory is passed.
    """

    if path[0] != os.sep:
        path = os.sep + path

    filename = ''
    extension = ''

    split_for_ext = path.split('.')
    if len(split_for_ext) > 1:
        extension = '.' + split_for_ext[-1]
    else:
        extension = ''

    # Remove trailOing '/'
    path = os.path.normpath(split_for_ext[0])

    if len(extension) > 0:
        filename = path.split(os.sep)[-1]
        path = os.path.join(*path.split(os.sep)[:-1])

    path = list(filter(None,path.split(os.sep)))

    folder = path[-1]
    base = os.sep + os.path.join(*path[:-1])

    return base, folder, filename, extension

I think all the options are self-explanatory. Provide the path to pad_geotiff and a pad width (npad). Optionally provide output filepath (outpath) and the value to use to pad (padval). It will write a padded geotiff and return the path of this geotiff (outpath).

You can pass in custom creation options to write_tile if you want to save with compression, set blocksize, etc. For example:

    options = ['COMPRESS=PACKBITS',
#               'NUM_THREADS=ALL_CPUS',
               'PREDICTOR=2',
               'BLOCKXSIZE=256',
               'BLOCKYSIZE=256',
               'TILED=YES']

Answer 2

Although quite old, I stumbled upon this page with the same question- How do I fill the extents of a raster to a defined tile grid? The problem I was having is that all grid squares at the edge of my raster that were partially filled by the raster extents were being missed by a gdal_translate loop. Seeing as it took me a while to figure out the issues I was having I thought I would share a solution in case anyone else finds their way here.

gdal_translate -projwin will process tiled areas and processing in a loop can give you chunked rasters to a grid.

The problem occurs when the xy's create a square outside the raster extents. My raster extents on the X axis runs to 650000, yet I need a 20k tile chunk that starts from 640000. I run the below:

gdal_translate -projwin 640000, 260000, 660000, 240000 fullraster.tif clippedraster.tif

If you recieve an error such as 'Computed -srcwin outside raster extents' then (double check your co-ordinates first) it is likely to be down to your GDAL versioning. If you are running a GDAL version below 1.10 this command will fail and no tif will be generated even though there is still unclipped valid raster data.

From GDAL version 1.10 gdal_translate will fill your raster extents automatically with no data values to your defined grid square. It will warn you that the computed -srcwin falls partially outside raster extents but will continue to process tif.

(You can use swithches -epo: (Error when Partially Outside) and -eco: (Error when Completely Outside) to manage the automatic functionality if required.)