Script de python para producir un ráster (*.asc) con la suma de los píxeles vecinos

En un foro sobre SIG un usuario tenía la necesidad de, a partir de los centroides de los píxeles de un ráster, obtener el valor de la sumatoria de los píxeles adyacentes a una distancia dada (tipo buffer) y que cada píxel adquiriera ese valor. Entre la asistencia recibida se sugería el uso de la herramienta de ‘Focal Statistics’ que se encuentra en el ‘Spatial Analyst’ de ArcGis. No obstante, en ausencia de cualquier otra información y asumiendo que los píxeles adyacentes están entre 3 (vértices), 5 (lados) u ocho (centro) en el ráster, es posible codificar un sencillo script de python para procesar la información a partir de un ráster en formato *.asc.

El script usa la matriz de píxeles del *.asc para generar otra que tiene un “anillo” de ceros y así garantizar que en el recorrido de todos los píxeles de ésta, especialmente en vértices y lados de la imagen, la evaluación de la suma se realice siempre para los 8 píxeles adyacentes. El script propuesto fue el siguiente:

#! /usr/bin/env python

#Librerias
import numpy as np
from StringIO import StringIO 
from os import system

system("clear")

pfile=open('raster.asc','r')

data=pfile.read()

pfile.close()

data=np.genfromtxt(StringIO(data)) #Se sobre entiende que los 
                                   #delimitadores son espacios

rows=len(data)
columns=len(data[0])

print data, "\n"

data_extend=np.zeros((rows+2,columns+2))

suma_vec=np.zeros((rows,columns))

for i in range(1,rows+1):
	for j in range(1,columns+1):
		data_extend[i][j]=data[i-1][j-1]

print data_extend

print

for i in range(rows):
	for j in range(columns):
		suma_vec[i][j] = data_extend[i:i+3,j:j+3].sum() - data[i][j]

print suma_vec

Su ejecución produce la siguiente salida:

[[ 2.  2.  1.  1.  5.  5.  5.]
 [ 2.  2.  8.  8.  5.  2.  1.]
 [ 7.  1.  1.  8.  2.  2.  2.]
 [ 8.  7.  8.  8.  8.  8.  5.]
 [ 8.  8.  1.  1.  5.  3.  9.]
 [ 8.  1.  1.  2.  5.  3.  9.]] 

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  2.  2.  1.  1.  5.  5.  5.  0.]
 [ 0.  2.  2.  8.  8.  5.  2.  1.  0.]
 [ 0.  7.  1.  1.  8.  2.  2.  2.  0.]
 [ 0.  8.  7.  8.  8.  8.  8.  5.  0.]
 [ 0.  8.  8.  1.  1.  5.  3.  9.  0.]
 [ 0.  8.  1.  1.  2.  5.  3.  9.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]

[[  6.  15.  21.  27.  21.  18.   8.]
 [ 14.  24.  24.  31.  33.  27.  16.]
 [ 20.  43.  50.  48.  49.  33.  18.]
 [ 31.  42.  35.  34.  37.  36.  24.]
 [ 32.  42.  36.  38.  38.  52.  28.]
 [ 17.  26.  13.  13.  14.  31.  15.]]

En ella se verifica, matriz intermedia, que a partir de la original se genera la del “anillo” de ceros y la última contiene la suma de los 8 píxeles adyacentes; tomando en cuenta sólo los diferentes de cero. Sin embargo, lo anterior sólo fue para probar el algoritmo. Corriendo el script con un ráster de verdad; como el DEM de la imagen siguiente:

se seleccionó, al azar, una matriz de 9 píxeles y se verificó que en ésta el píxel central de la matriz 3×3 contiene, en el ráster sum_vec.asc, el valor de la suma de sus vecinos en el DEM; tal como se ilustra en la imagen anterior. No obstante, el código se modificó ligeramente para evitar la impresión de matrices intermedias por pantalla y sólo imprimir el sum_vec.asc a un archivo. El código fue el siguiente:

#! /usr/bin/env python

#Librerias
import numpy as np
from StringIO import StringIO 
from os import system

system("clear")

pfile=open('tiznados_canoa2.asc','r')

print "Leyendo..."

data=pfile.read()

pfile.close()

print "Terminado..."

data=np.genfromtxt(StringIO(data)) #Se sobre entiende que los 
                                   #delimitadores son espacios

rows=len(data)
columns=len(data[0])

data_extend=np.zeros((rows+2,columns+2))

suma_vec=np.zeros((rows,columns))

print "Obteniendo matriz con anillo de ceros..."

for i in range(1,rows+1):
	for j in range(1,columns+1):
		data_extend[i][j]=data[i-1][j-1]

print "Terminado..."

print "Calculando suma y escribiendo a archivo..."

pfile=open('suma_vec.asc','w')

for i in range(rows):
	for j in range(columns):
		suma_vec[i][j] = data_extend[i:i+3,j:j+3].sum() - data[i][j]
		pfile.write(str(suma_vec[i][j],)+" ")
	pfile.write("\n") 	

pfile.close()

print "Listo..."

y produjo esta salida en poco más de 3 minutos y medio:

Leyendo...
Terminado...
Obteniendo matriz con anillo de ceros...
Terminado...
Calculando suma y escribiendo a archivo...
Listo...

real	3m31.640s
user	3m20.129s
sys	0m2.540s
zeito@debian:~$ 

About these ads