Géoréférencement raster utilisant GDAL et Python?

Je souhaite géoréférencer un raster à l'aide de pythonet GDAL. Mon approche actuelle est d'appeler gdal_translateet d' gdalwarputiliser os.systemet une vilaine liste de points de contrôle au sol. J'aimerais vraiment un moyen de le faire nativement à l'intérieur python.

Voici le processus que j'utilise actuellement:

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"')

Il y a une précédente question et réponse de 2012 qui indique que l' gdal_translateon peut y accéder après l'importation gdal. Je ne sais pas si est obsolète, ou si c'est faux, mais quand je cours, from osgeo import gdalje ne vois pas gdal.gdal_translated'option.

Je ne sais pas si elle existe mais j'aimerais beaucoup pouvoir traduire et reprojeter des rasters de manière pythonique. Par exemple:

# 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*)

Une telle approche est-elle possible?

Voici un exemple qui fait à peu près ce que vous demandez. Les principaux paramètres sont le geotransformtableau que gdal utilise pour décrire un emplacement de trame (position, échelle de pixels et asymétrie) et le epsgcode de la projection. Avec cela, le code suivant doit géoréférencer correctement le raster et spécifier sa projection.

Je n'ai pas beaucoup testé, mais cela a semblé fonctionner pour moi. Vous devez vous assurer que les coordonnées que j'ai saisies sont correctes et probablement changer la projection et l'échelle / taille des pixels. J'espère que cela aide.

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