Proyecciones GeoTIFF con Python GDAL

Question

Tengo un GeoTIFF que importé en el entorno de Python con GDAL. La imagen tiene datos de referencia espacial, que son los parámetros que describen el modelo geográfico:

(-104.31292762744124, 0.00030593847623094916, 0.0, 32.77290675637778, 0.0, -0.0002583980094641447)

No hay metadatos de proyección ( raster.GetProjection() ). ¿Cómo puedo utilizar estas coordenadas para crear una proyección (quiero coordenadas del mundo real para cada ubicación de los píxeles)?

He mirado en osr.SpatialReference y pyproj pero no he podido encontrar ninguna explicación sencilla o rápida. Lo siento si esto es básico, pero la documentación en línea simplemente no fue útil para mí.

Answer 1

Esta es la matriz de transformación afín que describe (en este orden en su caso): [x coordinate of top left corner, pixel width, rotation on x axis, y coordinate of top left corner, rotation on y axis, pixel height] . Afortunadamente, esto es suficiente para hacer lo que quieres.

Puede que te resulte un poco más fácil de usar rasterio que es una interfaz para GDAL mucho más fácil de usar que los enlaces de GDAL en Python.

Así es como funcionaría:

# you import rasterio like this
import rasterio

# you can open a dataset like this
with rasterio.open("myfile.tif") as dataset:

   # you can see the affine transformation matrix (that you had above) like this:
   print(dataset.transform)

   # you can convert from pixel coordinates to geographical like this:
   # (this example uses the middle of the dataset)
   x, y = dataset.xy(dataset.height // 2, dataset.width // 2)

   # and vice-versa like this:
   # (replace x and y with some coordinates that fall within your dataset)
   row, col = dataset.index(x, y)

Tenga en cuenta que Espacio de la imagen (píxel) están en el orden y,x (fila, col), mientras que espacio de coordenadas utiliza x, y (lng lat / easting northing). También conviene saber que Espacio de la imagen tiene el origen en la parte superior izquierda, donde Espacio de coordenadas tiene el origen en la parte inferior izquierda, (de ahí que la altura del píxel sea negativa en su conjunto de datos).

Si quieres añadir el CRS también (no es necesario para lo anterior) - podrías hacerlo:

   # assign CRS as WGS84
   dataset.crs = rasterio.crs.CRS.from_string('EPSG:4326')

He adivinado su Sistema de Referencia de Coordenadas (CRS) allí - si es uno diferente a WGS84, entonces usted puede obtener el código EPSG correcto de aquí y cambiarlo por 4236 en el fragmento anterior.

---

Alternativamente - si no quieres cambiar a usar rasterio, podrías hacer la conversión desde ubicaciones de píxeles ( px_x , px_y ) a coordenadas geográficas utilizando su matriz de transformación afín ( t ) así:

# convert px_x and px_y rto geographical coordinates
x, y = (topleft_x + (px_x * resolution_x), topleft_y - (px_y * resolution_x))

# so in your case, if t is the affine transformation object you had:
x, y = (t[0] + (px_x * t[1]), t[3] - (px_y * t[5]))

Answer 2

Suponiendo que conoces la referencia espacial de tus datos (supongo que es WGS84 echando un vistazo a las coordenadas de origen que proporcionaste y al tamaño de los píxeles), puedes establecer la referencia espacial de esta manera:

from osgeo import gdal
from osgeo import osr

# create WGS84 Spatial Reference
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)

# open raster in writing mode
fn = r'C:\path\to\my_raster.tif'
ds = gdal.Open(fn, 1)

# set spatial reference
ds.SetProjection(sr.ExportToWkt())

# save raster
ds = None