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

En un post anterior se obtuvo un ráster con la distancia media a los 5 puntos más cercanos de una línea con 8 puntos. No obstante, el requerimiento consistía en hacerlo a una línea de 300 puntos. Después de sustituir el ineficiente método de la burbuja por el método sort de python, escribir al array de un *.tif en lugar de hacerlo en un *.asc, todo mediante la librería GDAL de Python, generé trescientos puntos aleatorios para la zona del ráster, le añadí sus coordenadas en la tabla atributiva y luego las copié en el archivo de entrada de datos desde donde se leen los puntos de la recta. Al correr el programa éste tardó 7 horas y cuarto en producir el ráster. El resultado, por cónsola, fue el siguiente:

columnas = 3603 filas = 3630
Min=0.000, Max=1876.000
band.XSize =  3603 band.YSize =  3630
Leyendo datos de py_linea.dat ...
Listo...
Ejecutando...
Min=3060.508, Max=19854.041
Listo...

real	435m31.506s
user	432m9.637s
sys	0m10.217s
zeito@debian:~$ 

y la visualización del ráster, en QGIS, se da a continuación:

media5

El código python fue el siguiente:

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

#Definición de funciones

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

def promedio_5_Lista(lista):

	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 

# 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")

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 = dataset.GetGeoTransform()

#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

# Set file vars
output_file = "out.tif"

# Create gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, band.XSize, band.YSize, 1, gdal.GDT_Float32)

print "band.XSize = ", band.XSize, "band.YSize = ", band.YSize

pfile2 = open('py_linea.dat', 'r')

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

pfile2.close()

print "Listo..."

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] 

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

print "Ejecutando..."

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

		A = Distancia(X, Y, x, y, larray)
		A.sort()
		B = promedio_5_Lista(A)
		raster[i][j] = B
		x += geotransform[1] #tamaño de celda del ráster 
	x = geotransform[0]
	y -= -geotransform[5] #tamaño de celda del ráster

dst_ds.GetRasterBand(1).WriteArray( raster )

dataset = None

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

# 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] ] )
 
# set the reference info 
srs = osr.SpatialReference()
srs.SetUTM( 19, 1 )
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
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

dst_ds = None

print "Listo..."

Considerando que el ordenamiento de 8 valores para cada píxel, en 13.078.890 de ellos, se realizó en alrededor de 14 minutos, una extrapolación grosera a 300 valores da una estimación de 8 horas con 3/4 de computo. Comparado con el tiempo real que fue 7 horas con 1/4 permite señalar que el ordenamiento parece ser el que determina el tiempo de ejecución del programa.

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

Una respuesta a 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

  1. david arturomendez matute dijo:

    muchas gracias me gustaria que me enviara un com hacer ejecutar un raster de 4 puntos en mi computador david mendez mi correo :davidarturomendezmatute@yahoo.com

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