Cómo establecer el CRS de un ráster desde la Python Consola de QGIS (PyQGIS) con GDAL

Algunas veces se observa el requerimiento de querer especificar el CRS (Coordinate Reference System) de un ráster que, por alguna extraña razón, no está ya encapsulado en un formato “embebido” como el GTiff. Se me ocurre pensar que podría ser cuando se produce un ráster en un geoproceso (usando algún lenguaje de programación) y “olvidamos” hacer tal asignación mediante las líneas de código correspondiente. En esta fase de creación del ráster, generalmente, es fácil asignar el CRS porque, por defecto, está habilitada la escritura de información al ráster.

Cuando se omite el CRS en un ráster y conocemos precisamente esta información que falta, su incorporación ahora requiere que se habilite expresamente el ráster para escritura. Para probar el procedimiento de asignación del CRS generé, mediante un script de Python, un raster con driver GTiff al cual deliberadamente le omití esa información. Ese raster (out_pyqgis.tif) tiene todos sus values iguales a uno y su extensión coincide con un DEM ubicado en la zona UTM 12N y del cual tomé su metadata (opciones de geotransform para GDAL en Python). Lo pueden bajar de aquí para probar el procedimiento. Si lo abren en QGIS verán como les aparece una ventana que señala que el CRS será tomado del proyecto actual (“CRS was undefined”).

Una vez bajado, lo colocan en el directorio de su preferencia, abren la Python Console de QGIS y ejecutan secuencialmente estas líneas de código (también pueden usar el Editor):

from osgeo import osr, gdal
from gdalconst import *
dataset = gdal.Open('/home/zeito/out_pyqgis.tif', GA_Update)
band = dataset.GetRasterBand(1)
srs = osr.SpatialReference()
srs.ImportFromEPSG(32612)
dataset.SetProjection( srs.ExportToWkt() )
dataset = None

En las líneas 2 y 3 se habilita el raster para escritura (con el path que escogí para él) y en la línea 7 se asigna la proyección cuyo EPSG es 32612 (UTM 12N). Carguen el ráster en un proyecto nuevo de QGIS y observarán que ya no se despliega el señalamiento de “CRS was undefined” porque el CRS original de éste ya ha sido asignado.

Esta entrada fue publicada en GDAL, PyQGIS, QGIS, Software Libre. Guarda el enlace permanente.

2 respuestas a Cómo establecer el CRS de un ráster desde la Python Consola de QGIS (PyQGIS) con GDAL

  1. Observación: En Brighton el código no corría, tuve que cambiar en la linea 2 a:
    from osgeo.gdalconst import *

    Saludos.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s