python - Up-sampling (increasing resolution) raster image using GDAL?

I have a single-band geo-referenced raster image (a DEM) and my goal is to increase the number of pixels in each dimension (x and y) by 3. So given an initial resolution of 1200 * 1200 pixels, I want to upsample it to 3600 * 3600 pixels, keeping the geographic referencing.

This procedure is a part of a large batch processing so simple gui workflow is not an option and I'd like to do it programmatically, using python gdal.

What is the right way to do it and is there any interpolation required to do it?

Answer

The following gdal script is useful to resample an image to a smaller pixel size:

import os
from osgeo import gdal


# Change working directory     
os.chdir("directory with rasters")

# Open raster and get band
in_ds = gdal.Open('raster')
in_band = in_ds.GetRasterBand(1)

# Multiply output size by 3 
out_rows = in_band.YSize * 3
out_columns = in_band.XSize * 3


# Create new data source (raster)
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('band1_resampled.tif', out_columns, out_rows)
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())

# Edit the geotransform so pixels are one-sixth previous size
geotransform[1] /= 3
geotransform[5] /= 3

out_ds.SetGeoTransform(geotransform)

data = in_band.ReadAsArray(buf_xsize=out_columns, buf_ysize=out_rows)  # Specify a larger buffer size when reading data
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)

out_band.FlushCache()
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32, 64])


del out_ds

However, this script does not perform a specific interpolation and the result will look similar to the following picture (with a different resampling size):

Note: image taken from Geoprocessing with Python (by Chris Garrard) book. https://www.manning.com/books/geoprocessing-with-python

Furthermore, you could try using gdal_translate from the gdal command line utilities. (More info here: http://www.gdal.org/gdal_translate.html)

As you need to do a large processing batch it is possible to use python along with this utility as the following example:

import os
import subprocess


os.chdir("directory with the rasters")

result = subprocess.call('gdal_translate -of GTiff -outsize 3600 3600 -r bilinear raster.tif raster_resample.tif')

where:

• - of specifies the output format.


• -outsize specifies the output size in pixels (xsize, ysize)


• -r specifies the resampling algorithm


• raster.tif is the input filename



• raster_resample.tif is the output filename.