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.

Anuncios
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