5

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.

enter image description here

  • Rather than adapt the image to the tile system how about editing the tile index to match image? (Old question, I know) – matt wilkie Mar 12 '14 at 17:39
  • Hey @mattwilkie, thanks for giving it a look. In this case, the tile index/grid is the TMS structure, so each tile corresponds to the TMS bounds. So for those tiles along the top, I'd want to render tiles that had transparency not only where the image is already NODATA, but also where the tile bounds exceeded the original raster. I probably should have mentioned the TMS tiling constraint in the original question. :/ – elrobis Mar 12 '14 at 18:43
  • 1
    you might be able to do something with gdalbuildvrt -te, set target extent, and -hidenodata, "...Useful when you want to control the background color of the dataset. By using along with the -addalpha option, you can prepare a dataset which doesn't report nodata value but is transparent in areas with no data." – matt wilkie Mar 12 '14 at 18:47
3

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']
  • Hey thanks for adding this. I put this question aside quite some time ago but every now and then think about revisiting it. I'm sure yours is a useful contribution. If and when I return to this I'll certainly investigate some of your code there, as I have a hunch it may reveal the clues I need to modify gdal_retile.py to handle my use case (i.e. tile a county-wide, high resolution image to fit the TMS grid). – elrobis Oct 18 '17 at 20:03
  • 1
    My pleasure. After reading your problem more closely, I don't think this will help you at all (unless you decide to loop through each tile and pad it, but even then you'd probably want a different approach). If I were you, I would go with your first option: make a raster of global extents filled with 0's (or whatever nodata value you want). Then you can create a virtual raster (gdalbuildvrt) with this global one and all your tiles, setting your nodata value appropriately. Finally, you can just ReadAsArray() from that vrt. – Jon Oct 18 '17 at 20:18
  • Ok cool, interesting points here. You're planning on leaving your answer here though, right? :) – elrobis Oct 18 '17 at 20:21
  • Probably I shouldn't as it doesn't directly answer your question, but I will because I assume anyone trying to solve the problem I was trying to solve will find this question, just as I did :) And it is at least tangentially relevant. – Jon Oct 18 '17 at 20:22
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.)

  • Hey there--yes I'm definitely still interested in this problem, although I've got it paused just now. If I'm understanding you correctly, are you looping over your tile grid and passing the tile extents into gdal_translate to get each tile? – elrobis Feb 9 '16 at 15:40
  • At the time, I was trying to modify the gdal_retile utility's source code, because that utility runs ultra fast. Essentially I just wanted to add an option to snap it to the TMS tile grid, but I ran into trouble somewhere at band.ReadRaster(), because the TMS grid fell outside of the actual raster extents. – elrobis Feb 9 '16 at 15:45
  • 1
    Yes I am passing the tile grid co-ordinates in a loop to chunk a raster image. So, starting from bottom left (in my case 0,0) apply 'gdal_translate -projwin $((x*20000)) $(((y+1)*20000)) $(((x+1)*20000)) $((y*20000))' to each subsequent tile. Would this work for you? – ProudGIS Feb 10 '16 at 11:34
  • It might--it's certainly an intriguing alternative. I've up-voted your answer for now, but I'm going to wait to accept until I can spend a little more time researching a comparable technique using the GDAL python bindings (like the rest of gdal_retile does). However I'll make a point to try this too--if it's fast enough I may just accept it. I appreciate the help and insight. – elrobis Feb 11 '16 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.