10

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?

  • 1
    gdal_translate and gdalwarp and other gdal utilities are not python packages/modules, they are standalone executables. You can use os.system or (preferably) subprocess.Popen to call them. – user2856 Oct 9 '14 at 0:55
  • @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. – djq Oct 9 '14 at 3:15
  • You can use the gdal python API to do nearly everything the gdal_translate and gdalwarp utilities can, it's just a bit more complicated. Check out gdal.AutoCreateWarpedVRT and gdal Driver.CreateCopy. – user2856 Oct 9 '14 at 7:01
16

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
  • 1
    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. – CaMa Mar 27 '17 at 15:05

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.