¿Ráster de georreferenciación con GDAL y Python?

10

Quiero georreferenciar un ráster usando pythony GDAL. Mi enfoque actual es llamar gdal_translatey gdalwarpusar os.systemy una lista fea de puntos de control en tierra. Realmente me gustaría una forma de hacer esto de forma nativa python.

Este es el proceso actual que estoy usando:

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

Hay una pregunta y respuesta anterior de 2012 que establece que gdal_translatese puede acceder después de la importación gdal. No estoy seguro de si está obsoleto o si está mal, pero cuando lo ejecuto from osgeo import gdalno lo veo gdal.gdal_translatecomo una opción.

No sé si existe, pero me encantaría poder traducir y reproyectar rásteres de manera pitónica. Por ejemplo:

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

¿Es posible este enfoque?

djq
fuente
1
gdal_translate y gdalwarp y otras utilidades de gdal no son paquetes / módulos de Python, son ejecutables independientes. Puede usar os.system o (preferiblemente) subprocess.Popen para llamarlos.
user2856
@Luke, entonces, ¿no hay forma de interactuar con estas utilidades de gdal? Me encantaría explorar cómo escribir un contenedor de python para estos si no existe.
DJ
Puede usar la API gdal python para hacer casi todo lo que pueden hacer las utilidades gdal_translate y gdalwarp, es un poco más complicado. Echa un vistazo a gdal.AutoCreateWarpedVRT y gdal Driver.CreateCopy.
usuario2856

Respuestas:

17

Aquí hay un ejemplo que hace aproximadamente lo que pides. Los parámetros principales son la geotransformmatriz que usa gdal para describir una ubicación ráster (posición, escala de píxeles y sesgo) y el epsgcódigo de la proyección. Con eso, el siguiente código debería georreferenciar correctamente el ráster y especificar su proyección.

No probé tanto, pero parecía funcionar para mí. Tendría que asegurarse de que las coordenadas que puse sean correctas y probablemente cambien la proyección y la escala / tamaño de píxeles. Espero eso ayude.

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
tapa amarilla
fuente
2
Es correcto, pero si desea ser más preciso, en caso de rotación, debe usar la transformación Affine y calcular los parámetros correctos con el paquete Affine ( descargue desde aquí ). Uso: 'Affine.scale (PARAMETER) * Affine.rotation (ANGLE)'. PARÁMETRO son de la relación de píxeles de las dos imágenes, puede usar 'gdalinfo' o PIL para conocer el tamaño de píxel, ANGLE es el ángulo.
CaMa