Crear un ráster *.tif copiando parte de la metadata de un ráster *.asc mediante librería GDAL en Python

En el post anterior se produjo un ráster *.asc que reflejaba la distancia media a los 5 puntos más cercanos de una línea. Aunque el resultado fue satisfactorio, la manera como manipulé su escritura no me convenció del todo porque “mezclé procedimientos” (no emplee exclusivamente la librería de gdal) para copiar lo que me interesaba al ráster de salida. Además, como el driver AAIGrid (*.asc) no permite la opción de Create, sólo copia, entonces se incluye íntegramente toda la metadata del ráster padre y no me fue posible añadir sus propias estadísticas. Esto trae como consecuencia que sólo admite la paleta con tonos de grises para el ensanchamiento de contraste y para visualizarlo adecuadamente tuve que exportarlo al ambiente de GRASS con la consecuente pérdida de tiempo.

Indagando en estas referencias:

Creating a multispectral image from scratch

gdal.h File Reference

pude encontrar la manera de escribir directamente al array de datos de un *.tif sin preocuparme por el encabezado del ráster (algo que posiblemente también es factible de hacer en el *.asc pero habría que averiguarlo). La razón es porque el driver GTiff, a diferencia del driver AAIGrid, si admite la opción Create. Por tanto, posteriormente, sólo copio lo que me interesa: la extensión y la proyección del ráster (además de los valores de los píxeles).

El código completo es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8

# Import libs
import numpy as np
from osgeo import osr, gdal
from os import system
from StringIO import StringIO
from math import sqrt

system ("clear")

#Abre el ráster original
dataset = gdal.Open('tiznados_canoa.asc')

band = dataset.GetRasterBand(1)

print "columnas = %d filas = %d" % (band.XSize, band.YSize)

#Se obtiene la metadata del ráster

#geotransform[0] = top left x
#geotransform[1] = w-e pixel resolution
#geotransform[3] = top left y
#geotransform[5] = n-s pixel resolution

geotransform = dataset.GetGeoTransform()

#Se obtienen los valores mínimo y máximo del ráster padre
min = band.GetMinimum()
max = band.GetMaximum()

#Determina los valores mínimo y máximo si no existen
if min is None or max is None:
	(min,max) = band.ComputeRasterMinMax(1)

print "Ráster padre:\nMin=%.3f, Max=%.3f" % (min,max) #Imprime valores mínimo y máximo

# Establece el nombre del ráster de salida
output_file = "out2.tif"

# Crea el gtif con las filas y columnas del ráster padre y tipo float 
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, band.XSize, band.YSize, 1, gdal.GDT_Float32)

#Declara el arreglo del raster y le asigna 1 a todos los píxeles
raster = np.ones((band.YSize,band.XSize))

#Escribe al ráster de salida
dst_ds.GetRasterBand(1).WriteArray( raster )

#Cierra el dataset del ráster padre
dataset = None

#Establece la extensión del ráster de salida
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( [ geotransform[0], geotransform[1], 0, geotransform[3], 0, geotransform[5] ] )
 
# Establece la referencia espacial del ráster de salida 
srs = osr.SpatialReference()
srs.SetUTM( 19, 1 ) # Zona UTM 19: El uno es true para el hemisferio norte
srs.SetWellKnownGeogCS( 'WGS84' )
dst_ds.SetProjection( srs.ExportToWkt() )

band = dst_ds.GetRasterBand(1)

#Se obtienen los valores mínimo y máximo del ráster de salida
min = band.GetMinimum()
max = band.GetMaximum()

#Determina los valores mínimo y máximo si no existen
if min is None or max is None:
	(min,max) = band.ComputeRasterMinMax(1)

print "Ráster de salida:\nMin=%.3f, Max=%.3f" % (min,max) #Imprime valores mínimo y máximo

#Cierra el dataset del ráster de salida
dst_ds = None

Su ejecución en cónsola produce esta salida:

columnas = 3603 filas = 3630
Ráster padre:
Min=0.000, Max=1876.000
Ráster de salida:
Min=1.000, Max=1.000

real	0m0.916s
user	0m0.296s
sys	0m0.304s
zeito@debian:~$

y su visualización con QGIS, con transparencia de 30 % sobre el ráster original, es la siguiente:

gdal

donde se corrobora que tiene unos (1) en todos sus píxeles.

Esta entrada fue publicada en Código Python, Teledetección. Guarda el enlace permanente.

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