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.