16

I have a raster (USGS DEM actually) and I need to split it up into smaller chunks like the image below shows. That was accomplished in ArcGIS 10.0 using the Split Raster tool. I would like a FOSS method to do this. I've looked at GDAL, thinking surely it would do it (somehow with gdal_translate), but can't find anything. Ultimately, I'd like to be able to take the raster and say how large (4KM by 4KM chunks) I would like it split up into.

enter image description here

  • I have a utility that uses subprocess.Popen to run multiple gdal translates at the same time which I use for extracting a large raster to tiles using a fishnet, particularly useful if the input and/or output is highly compressed (eg LZW or deflate GeoTiff), if neither is highly compressed the process peaks on HDD access and isn't much faster than running one at a time. Unfortunately it's not generic enough to share due to rigid naming conventions but food for thought anyway. The -multi option for GDALWarp often causes trouble and only uses 2 threads (one read, one write) not all available. – Michael Stimson Oct 18 '18 at 4:50
18

gdal_translate will work using the -srcwin or -projwin options.

-srcwin xoff yoff xsize ysize: Selects a subwindow from the source image for copying based on pixel/line location.

-projwin ulx uly lrx lry: Selects a subwindow from the source image for copying (like -srcwin) but with the corners given in georeferenced coordinates.

You would need to come up with the pixel/line locations or corner coordinates and then loop over the values with gdal_translate. Something like the quick and dirty python below will work if using pixel values and -srcwin is suitable for you, will be a bit more work to sort out with coordinates.

import os, gdal
from gdalconst import *

width = 512
height = 512
tilesize = 64

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(tilesize)+", " \
            +str(tilesize)+" utm.tif utm_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)
  • 1
    Hi when I try -projwin option with a geotiff image I get warning saying "Warning: Computed -srcwin -3005000 1879300 50 650 falls completely outside raster extent. Going on however" I am not sure where I am doing wrong seems like it does not using its georeferenced coordinates. – ncelik Oct 26 '16 at 9:06
  • @ncelik that is probably because you're using cell coordinates in your projwin and should be using srcwin instead. If you are experiencing difficulty please post a new question with all the relevant information so we can make suggestions on your specific problem. – Michael Stimson Oct 18 '18 at 4:52
15

My solution, based on the one from @wwnick reads the raster dimensions from the file itself, and covers the whole image by making the edge tiles smaller if needed:

import os, sys
from osgeo import gdal

dset = gdal.Open(sys.argv[1])

width = dset.RasterXSize
height = dset.RasterYSize

print width, 'x', height

tilesize = 5000

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        w = min(i+tilesize, width) - i
        h = min(j+tilesize, height) - j
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(w)+", " \
            +str(h)+" " + sys.argv[1] + " " + sys.argv[2] + "_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)
  • I think it should be sys.argv[1] where it says sys.argv[2], right? – oskarlin Jul 30 '17 at 17:00
  • 3
    sys.argv[2] is used as a prefix for the output files, I believe. Super helpful-- thanks @Ries! – Charlie Hofmann Aug 8 '17 at 14:40
3

You can use r.tile of GRASS GIS. r.tile generates a separate raster map for each tile with numbered map names based on the user defined prefix. Width of tiles (columns) and height of tiles (rows) can be defined.

Using the grass-session Python API only a few lines of Python code are needed to call the r.tile functionality from outside, i.e. to write a standalone script. Using r.external and r.external.out also no data duplication occurs during the GRASS GIS processing step.

Pseudo code:

  1. initialize grass-session
  2. define output format with r.external.out
  3. link input file with r.external
  4. run r.tile which generates the tiles in the format defined above
  5. unlink r.external.out
  6. close grass-session
3

There's a bundled python script specifically for retiling rasters, gdal_retile:

gdal_retile.py [-v] [-co NAME=VALUE]* [-of out_format] [-ps pixelWidth pixelHeight]
               [-overlap val_in_pixel]
               [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
                      CInt16/CInt32/CFloat32/CFloat64}]'
               [ -tileIndex tileIndexName [-tileIndexField tileIndexFieldName]]
               [ -csv fileName [-csvDelim delimiter]]
               [-s_srs srs_def]  [-pyramidOnly]
               [-r {near/bilinear/cubic/cubicspline/lanczos}]
               -levels numberoflevels
               [-useDirForEachRow]
               -targetDir TileDirectory input_files

e.g:

gdal_retile.py -ps 512 512 -targetDir C:\example\dir some_dem.tif

3
+100

For @Aaron who asked:

I am hoping to find a gdalwarp version of @wwnick's answer that utilizes the -multi option for enhanced multicore and multithreaded operations

Slight Disclaimer

This uses gdalwarp, though I'm not totally convinced there will be much performance gain. So far I've been I/O bound - running this script on a large raster cutting it into many smaller parts doesn't seem CPU intensive, so I assume the bottleneck is writing to disk. If you're planning on simultaneously re-projecting the tiles or something similar, then this might change. There are tuning tips here. A brief play didn't yield any improvement for me, and CPU never seemed to be the limiting factor.

Disclaimer aside, here's a script that will use gdalwarp to split up a raster into several smaller tiles. There might be some loss due to the floor division but this can be taken care of by picking the number of tiles you want. It will be n+1 where n is the number you divide by to get the tile_width and tile_height variables.

import subprocess
import gdal
import sys


def gdalwarp(*args):
    return subprocess.check_call(['gdalwarp'] + list(args))


src_path = sys.argv[1]
ds = gdal.Open(src_path)

try:
    out_base = sys.argv[2]
except IndexError:
    out_base = '/tmp/test_'

gt = ds.GetGeoTransform()

width_px = ds.RasterXSize
height_px = ds.RasterYSize

# Get coords for lower left corner
xmin = int(gt[0])
xmax = int(gt[0] + (gt[1] * width_px))

# get coords for upper right corner
if gt[5] > 0:
    ymin = int(gt[3] - (gt[5] * height_px))
else:
    ymin = int(gt[3] + (gt[5] * height_px))

ymax = int(gt[3])

# split height and width into four - i.e. this will produce 25 tiles
tile_width = (xmax - xmin) // 4
tile_height = (ymax - ymin) // 4

for x in range(xmin, xmax, tile_width):
    for y in range(ymin, ymax, tile_height):
        gdalwarp('-te', str(x), str(y), str(x + tile_width),
                 str(y + tile_height), '-multi', '-wo', 'NUM_THREADS=ALL_CPUS',
                 '-wm', '500', src_path, out_base + '{}_{}.tif'.format(x, y))

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.