Producir un ráster con la distancia media a los 5 puntos más cercanos de una línea

Siguiendo con el problema de encontrar en un ráster la distancia media a los 5 puntos más cercanos de una línea con 300 puntos y luego volcarlos, para cada píxel, en un nuevo ráster se produjo el siguiente código:

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

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

from osgeo import gdal,ogr

system("clear")

print "Posibles rásteres a procesar"

system("ls *.asc")

#Definición de funciones

#Función con ordenamiento tipo burbuja

def ordenamientoBurbuja(lista,tam):
    for i in range(1,tam):
        for j in range(0,tam-i):
            if(lista[j] > lista[j+1]):
                k = lista[j+1]
                lista[j+1] = lista[j]
                lista[j] = k;

#Función que promedia los 5 primeros valores del arreglo ordenado

def promedio_5_Lista(lista,tam):

	sum = 0
	
	for i in range(0, 5):
		sum += lista[i]	
	
	return sum/5

#Función que determina la distancia a n puntos con coordenadas x,y

def Distancia(x_l, y_l, x_p, y_p, l_array):
	d = np.zeros(l_array)
	for i in range(l_array):	
		d[i] = sqrt((x_l[i] - x_p)*(x_l[i] - x_p) + (y_l[i] - y_p)*(y_l[i] - y_p))

	return d 
 
#Comienzo del programa

#Se introduce el nombre del ráster

nameraster = raw_input(str("\nIntroduzca el nombre del ráster = ? "))

dataset = gdal.Open(nameraster)

band = dataset.GetRasterBand(1)

#Se obtienen los valores mínimo y máximo del ráster
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 "Min=%.3f, Max=%.3f" % (min,max) #Imprime valores mínimo y máximo

print 'Driver: ', dataset.GetDriver().ShortName

format = dataset.GetDriver().ShortName
driver = gdal.GetDriverByName(format)

print 'Band Type=', gdal.GetDataTypeName(band.DataType)

#Se obtiene la metadata del ráster
geotransform = dataset.GetGeoTransform()

if not geotransform is None:
	print 'Origen = (',geotransform[0], ',',geotransform[3],')'
	print 'Tamaño de Pixel = (',geotransform[1], ',',geotransform[5],')'

pfile1=open(nameraster,'r')

print "Leyendo datos de %s ..." % nameraster
data = pfile1.read()
 
pfile1.close()

data = data.split("\n")

for i in range(5):
	print data[i]

print "columnas = %d filas = %d" % (band.XSize, band.YSize)
 
pfile2 = open('py_linea.dat', 'r')

print "Leyendo datos de py_linea.dat ..." 
data2 = pfile2.read()

pfile2.close()

data2=np.genfromtxt(StringIO(data2))

#Determina longitud del arreglo de puntos de la recta
larray = len(data2)

X = np.zeros(larray)
Y = np.zeros(larray)

for i in range(larray):
	X[i] = data2[i][0]
	Y[i] = data2[i][1]

#Coordenada x del vértice superior izquierdo del ráster
x = geotransform[0]
#Coordenada y del vértice superior izquierdo del ráster
y = geotransform[3] 

print "Espere..."

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

dst_ds = driver.CreateCopy('media_5min.asc', dataset, 0 )

d_5min = np.zeros((band.XSize,band.YSize))

#Imprime el encabezado del ráster de salida
for i in range(5):
	pfile.write(data[i]+"\n")

for i in range(band.YSize):
	for j in range(band.XSize):

		A = Distancia(X, Y, x, y, larray)
		ordenamientoBurbuja(A,len(A))
		B = promedio_5_Lista(A,len(A))
		d_5min[i][j] = B
		pfile.write(str(d_5min[i][j],)+" ")
		x += geotransform[1] #tamaño de celda del ráster 
	pfile.write("\n")
	x = geotransform[0]
	y -= -geotransform[5] #tamaño de celda del ráster

pfile.close()

# Una vez hecho, se cierra adecuadamente el dataset
dst_ds = None
dataset = None

print "Listo!"

El programa, en términos generales, opera de la siguiente manera:

1. Se usa la librería del sistema (os) para borrar la pantalla y producir un listado de los posibles archivos ráster, en formato *.asc, a procesar.

2. Con la librería de gdal se tienen las herramientas para abrir el dataset y obtener la banda, el driver y la metadata necesaria para automatizar el proceso de barrido del ráster. La obtención del driver permitió conocer que los ráster AAIGrid admiten copia en modo de escritura por lo que pueden ser modificados. Es una forma rápida de producir un ráster de salida con la proyección original pero tiene la desventaja de que, por ahora, no pude introducirle las estadísticas propias (preserva las del original).

3. El arreglo de puntos de la recta se lee a partir de un archivo externo (en este caso ocho) y se utiliza en dos bucles anidados que barren el ráster para determinar la distancia de todos los puntos a cada píxel, se ordenan de menor a mayor y se determina el promedio de los 5 menores antes de ser imprimidos al ráster de salida; al cual previamente se le colocaron las 5 líneas del encabezado del ráster original.

La ejecución del código produce está salida por cónsola:

zeito@debian:~$ time python py_ord5min.py

Posibles rásteres a procesar
media_5min2.asc  media_5min.asc  prueba.asc  suma_vec2.asc  tiznados_canoa2.asc
media_5min3.asc  prueba2.asc	 raster.asc  suma_vec.asc   tiznados_canoa.asc

Introduzca el nombre del ráster = ? tiznados_canoa.asc
Min=0.000, Max=1876.000
Driver:  AAIGrid
Band Type= Int32
Origen = ( 609585.502177 , 1106092.58697 )
Tamaño de Pixel = ( 30.6105447446 , -30.6105447446 )
Leyendo datos de tiznados_canoa.asc ...
Leyendo datos de py_linea.dat ...
Espere...
Listo!

real	21m38.110s
user	21m17.716s
sys	0m1.824s
zeito@debian:~$ 

y la visualización es la siguiente:

raster4

Es como una especie de buffer con 5 valores de referencia. Aunque funciona (pero no de la forma deseada), seguiré buscando la manera de introducir la proyección sin producir la copia y aligerar el proceso de escritura en otro tipo de ráster.

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

3 respuestas a Producir un ráster con la distancia media a los 5 puntos más cercanos de una línea

  1. david arturomendez matute dijo:

    BUEN DIA MUCHAS GRACIAS POR LA INFORMACION ME GUSTARIA QUE ME MANDARAS MAS EJERCICIOS DE UN RASTER DE 4 PUNTOS A DISTANCIA EN UNA IMAGEN O MAPA PARA YO PRACTICAR EN MI COMPUTADOR

  2. Pingback: Crear un ráster *.tif copiando parte de la metadata de un ráster *.asc |

  3. Pingback: Ráster con la distancia media a los 5 puntos más cercanos de una línea de 300 puntos mediante librería GDAL en Python |

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