Georeferencing raster using GDAL and Python?

Question

I want to georeference a raster using python and GDAL. My current approach is to call gdal_translate and gdalwarp using os.system and an ugly list of ground control points. I'd really like a way to do this natively within python.

This is the current process I am using:

import os
os.system('gdal_translate -of GTiff -gcp 1251.92 414.538 -7.9164e+06 5.21094e+06 -gcp 865.827 107.699 -7.91651e+06 5.21104e+06 "inraster.tif" "outraster1.tif"')
os.system('gdalwarp -r bilinear -tps -co COMPRESS=NONE "outraster2.tif" "outraster3.tif"')

There is a previous question and answer from 2012 which states that gdal_translate can be accessed after importing gdal. I'm not sure if is obsolete, or whether it is wrong but when I run from osgeo import gdal I do not see gdal.gdal_translate as an option.

I don't know if it exists but I would love if I could translate and reproject rasters in a pythonic way. For example:

# translate
gcp_points = [(1251.92, 414.538), (-7.9164e+06, 5.21094e+06)]
gdal.gdal_translate(in_raster, gcp_points, out_raster1)
# warp
gdal.gdalwarp(out_raster1, out_raster2, 'bilinear', args*)

Is such an approach possible?

• djq: @Luke so is there no way of interacting with these gdal utilities? I'd be happy to explore how to write a python wrapper for these if one does not exist.

Answer 1

Here is an example that does roughly what you ask for. The main parameters are the geotransform array that gdal uses to describe a raster location (position, pixel scale, and skew) and the epsg code of the projection. With that, the following code should properly georeference the raster and specify its projection.

I did not test this much, but it seemed to work for me. You would have to make sure the coordinates I put in are correct and probably change the projection and the pixel scale/size. Hope it helps.

from osgeo import gdal, osr
src_filename ='/path/to/source.tif'
dst_filename = '/path/to/destination.tif'
# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)
# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)
# Specify raster location through geotransform array
# (uperleftx, scalex, skewx, uperlefty, skewy, scaley)
# Scale = size of one pixel in units of raster projection
# this example below assumes 100x100
gt = [-7916400, 100, 0, 5210940, 0, -100]
# Set location
dst_ds.SetGeoTransform(gt)
# Get raster projection
epsg = 3857
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()
# Set projection
dst_ds.SetProjection(dest_wkt)
# Close files
dst_ds = None
src_ds = None

• CaMa: It is right but if you want to be more accurate, in case of rotation, you have to use the Affine Transformation and calculate the right parameters with Affine package (download from here). Usage: 'Affine.scale(PARAMETER)*Affine.rotation(ANGLE)'. PARAMETER are from the pixel ratio of the two image, you can use 'gdalinfo' or PIL to know the pixel size, ANGLE is the angle.
• g123456k: In addition, here is how I found "gt" for the solution: reference = r'D:\....tif' src = gdal.Open(reference) ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform(). And simply add gt = [ulx, xres, xskew, uly, yskew, yres ]
• Tarun Maganti: How do I do this with multiple GCP points? 2D array maintaining [(uperleftx, scalex, skewx, uperlefty, skewy, scaley),(...)]?