Script de python para la extracción por máscara de una imagen ráster preservando la resolución original (GRASS-QGIS)

Hasta ahora sólo se han considerado procedimientos para archivos vectoriales en el código que se ha propuesto usando el módulo grass.script de las librerías de python. Hoy le tocará el turno al tratamiento de los archivos ráster donde se hará una extracción por máscara de un DEM, mediante un archivo vectorial, preservando la resolución original del DEM. Para ello, primero se alineará el vectorial con el ráster a través de g.region, se rasterizará el vectorial por un valor constante y se volcarán a su interior los valores del DEM usando r.mapcalc para, finalmente, eliminar el ráster temporal y restituir la g.region original.

El código (py_raster.py) es el siguiente:

#! /usr/bin/env python

import os
import grass.script as grass

os.system ("clear")

c = grass.region()

print "Verificando la resolucion inicial de la region"
print "ewres: %f" % c['ewres']
print "nsres: %f" % c['nsres']

#Obteniendo el nombre de la dataset
d=grass.list_pairs('vect')

nom_dataset = d[0][1]

print "\nVectoriales en la dataset:\n"

data_vect=grass.mlist_grouped('vect')

data_vect=data_vect[nom_dataset]

h=3
for i in range(len(data_vect)):
	print data_vect[i],
	if i == h:
		print 
		h +=4

print "\n"

vectorial= raw_input("vectorial a rasterizar = ? ")

print "\nRasteres en la dataset:\n"

data_rast=grass.mlist_grouped('rast')

data_rast=data_rast[nom_dataset]

h=3
for i in range(len(data_rast)):
	print data_rast[i],
	if i == h:
		print 
		h +=4

print "\n"

raster = raw_input("raster de alineacion = ? ")

grass.use_temp_region() #Al finalizar el script recupera la region original de grass

#Estableciendo la nueva g.region alineada con el vectorial
grass.run_command("g.region", vect=vectorial, align=raster)

salida=vectorial + "_raster"

#Rasterizando el vectorial
grass.run_command("v.to.rast", input=vectorial, value="1", output=salida, use="val", overwrite="True", quiet="True")

#Haciendo la extraccion por mascara
grass.mapcalc("poligono_mascara=if($salida, $raster)", salida=salida, raster=raster)

#Eliminando archivo temporal salida
grass.run_command("g.remove", rast=salida, quiet="True")

print "Verificando la resolucion final de la region"
print "ewres: %f" % c['ewres']
print "nsres: %f" % c['nsres']

que ejecutado en cónsola de GRASS-QGIS produce la siguiente salida:

A continuación se observa que el script produjo el resultado deseado:

Esta entrada fue publicada en Código Python, GRASS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Script de python para la extracción por máscara de una imagen ráster preservando la resolución original (GRASS-QGIS)

  1. Pingback: Interpolación de mapas de precipitación usando QGIS |

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