8

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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Basile
  • 3,543
  • 3
  • 21
  • 49
  • 3
    In ArcGIS the tool is called Resample http://resources.arcgis.com/en/help/main/10.2/index.html#//00170000009t000000 I suggest you'll want to use Bilinear interpolation, cubic gouges valleys and enhances peaks.. – Michael Stimson Feb 13 '18 at 00:17
  • If you also need an ArcPy answer after trying the above then please ask a new question with your code attempt. – PolyGeo Feb 13 '18 at 04:42

2 Answers2

15

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):

enter image description here

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.
Marcelo Villa
  • 5,928
  • 2
  • 19
  • 38
1

If you can use rasterio, these modules might help:

https://rasterio.readthedocs.io/en/stable/topics/resampling.html

Leo
  • 952
  • 6
  • 23