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.
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
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)
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)
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:
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))
import gdal
import os
def create_tiles(tile_size, input_filename, in_path='/media/Data/avinash/input/',
out_path='/home/nesac/PycharmProjects/GIS_Data_Automation/data/tiles/'):
ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
output_filename = 'tile_'
xsize, ysize = (band.XSize, band.YSize)
tile_size_x, tile_size_y = tile_size
tile_list = {}
complete_x = xsize // tile_size_x
complete_y = ysize // tile_size_y
residue_x = xsize % tile_size_x
residue_y = ysize % tile_size_y
# for part A
for i in range(complete_x):
for j in range(complete_y):
Xmin = i * tile_size_x
Xmax = i * tile_size_x + tile_size_x - 1
Ymin = j * tile_size_y
Ymax = j * tile_size_y + tile_size_y - 1
# do patch creation here
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
tile_size_x) + ", " + str(tile_size_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
x_residue_count = 1
y_residue_count = 1
# for part B
for j in range(complete_y):
Xmin = tile_size_x * complete_x
Xmax = tile_size_x * complete_x + residue_x - 1
Ymin = j * tile_size_y
Ymax = j * tile_size_y + tile_size_y - 1
# do patch creation here
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
residue_x) + ", " + str(tile_size_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
# for part C
for i in range(complete_x):
Xmin = i * tile_size_x
Xmax = i * tile_size_x + tile_size_x - 1
Ymin = tile_size_y * complete_y
Ymax = tile_size_y * complete_y + residue_y - 1
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
tile_size_x) + ", " + str(residue_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
# for part D
Xmin = complete_x * tile_size_x
Ymin = complete_y * tile_size_y
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
residue_x) + ", " + str(residue_y) + " " + \
str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)
