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
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:
donde se corrobora que tiene unos (1) en todos sus píxeles.
