Unleash the power of GRASS GIS at US-IALE 2017/es

Este es un material preparado para un taller de la US-IALE 2017, que lleva el nombre de Unleash the power of GRASS GIS y se realizó el 11 de Abril de 2017 en Baltimore. Este taller es introductorio a las capacidades de procesamiento de GRASS GIS relevantes para la ecología de paisaje.

Software

Se usará GRASS GIS 7.2 y R (>=3.1).

MS Windows

Si no tiene R, por favor instálelo primero. Luego descargue los binarios de GRASS GIS 7.2 (version de 64-bit, o versión de 32-bit) de grass.osgeo.org. Durante la instalación puede descargar el conjunto de datos de muestra de Carolina del Norte, por favor seleccione esta opción. También puede descargar los datos después (vea la siguiente sección).

Mac OSX Instale GRASS GIS 7.2 usando homebrew osgeo4mac:

brew tap osgeo/osgeo4mac
brew install grass7

Ubuntu Linux

Instale GRASS GIS desde el administrador de paquetes:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install grass

Para otras distribuciones de Linux, por favor encuentre GRASS GIS en su administrador de paquetes.

Extensiones de GRASS GIS (addons)

Se usarán las siguientes extensiones:

• r.forestfrag
• r.sun.hourly
• r.diversity

Las puede instalar desde la GUI, o usando la línea de comandos, por ejemplo:

 g.extension r.diversity

Paquetes de R

install.packages(c("rgdal", "rgrass7"))

Datos para este taller

Para este taller, descargue los siguientes conjuntos de datos:

• Localización GRASS de muestra de North Carolina
• points_of_interest.zip

Además, para la parte de scripts de R, necesita:

• US Census Bureau: Límites cartográficos de los estados de EEUU
• PRISM 30-year temperatura media anual normal, malla de 4 km
• datos de elevación SRTM para EEUU

Para la parte de lidar:

• mosaico LAS para Raleigh
• versión txt del mosaico (todos los puntos) y y solamente los primeros puntos de retorno

Introducción a GRASS GIS

Aquí se provee una reseña del proyecto GRASS GIS grass.osgeo.org que puede ser útil para revisar si ud. es un usuario primerizo. Para este ejercicio no es necesario entender completamente cómo funciona GRASS GIS. Sin embargo, puede necesitar saber cómo colocar sus datos en el directorio de bases de datos correcto de GRASS GIS, así como algunas ideas básicas sobre las funciones de GRASS. Aquí introducimos los conceptos esenciales ecesarios para correr este tutorial:

Estructura de la base de datos espacial de GRASS GIS

GRASS usa una terminología y estructura de la base de datos única (base de datos GRASS) que es importante entender para correr este tutorial, ya que es posible que necesite colocar los datos (ej. Localización en una base de datos específica de GRASS. En lo que sigue revisaremos la terminología y se dan las direcciones paso a paso de cómo descargar y colocar sus datos en el lugar correcto.

• Una Base de datos espacial de GRASS GIS (GRASS database) consiste de un directorio con Localizaciones (proyectos) específicos donde los datos (datos capas/mapas) se guardan.
• Localización es un directorio con datos relacionados a una localidad geográfica o a un proyecto. Todos los datos dentro de una Localización tienen el mismos sistema de coordenadas de referencia.
• Directorio de mapas es una colección de mapas dentro de una Localización, contiene los datos relacionados a una tarea específica, usuario o subproyecto.

Creando una base de datos de GRASS para el tutorial

Por favor descargue la Localización de muestra de GRASS de North Carolina, ponga atención de la ubicación en la cual se guardan los archivos dentro de su computadora.

Ahora, cree (a menos que ya lo haya hecho) una carpeta llamada grassdata (base de datos de GRASS) en su carpeta Home (o Documentos). Descomprima los datos descargados en esta carpeta. Ahora debe teer la Localización nc_spm_08_grass7 en grassdata.

• 

  GRASS GIS 7.2 startup dialog with North Carolina sample dataset

Mostrar y explorar los datos

Ahora que tenemos los datos en la base de datos de GRASS correcta, lanzamos la Interfaz Gráfica de Usuario (GUI) en el Directorio de mapas user1.

La interfaz GUI permite mostrar los datos vectoriales y ráster, así como acercarse y alejarse. Otras opciones de exploración y visualización también son posibles usando consultas o añadiendo leyendas, por ejemplo. La captura de pantalla de abajo muestra cómo se pueden añadir diferentes capas de mapas (izquierda) y mostrar los metadatos de las capas de datos.

• 

  Añadir capa de mapa ráster.

• 

  Añadir leyeda de ráster.

• 

  Administrador de capas y Display de mapas. Las anotaciones muestran cómo añadir una capa ráster, hacer consultas, añadir leyenda.

• 

  Mostrar metadatos de mapa ráster al hacer click derecho en la capa.

Módulos de GRASS GIS

Una de las ventajas d eGRASS es la diversidad y cantidad de módulos que le permiten analizar los diferentes aspectos espaciales y temporales. GRASS GIS tiene más de 500 módulos distintos en la distribución núcleo y más de 230 módulos de extensión addon que pueden ser usados para preparar y analizar capas de datos.

La funcionalidad de GRASS está disponible a través de módulos (herramientas, funciones). Existe en GRASS una convención de nombrar a los módulos de acuerdo a su función, de la siguiente manera:

Prefijo	Función	Ejemplo
r.*	procesamiento ráster	r.mapcalc: álgebra de mapas
v.*	procesamiento vectorial	v.clean: limpieza topológica
i.*	procesamiento de imágenes	i.segment: reconocimiento de objetos
db.*	administración de bases de datos	db.select: seleccionar valores de una tabla
r3.*	procesamiento ráster 3D	r3.stats: estadísticas de ráster 3D
t.*	procesamiento de datos temporales	t.rast.aggregate: agregación temporal
g.*	administración general de datos	g.rename: renombrar mapa
d.*	pantalla	d.rast: mostrar mapa ráster

Estos son los grupos de módulos principales. Hay algunos otros que son específicos. Note también que algunos módulos tienen varios puntos dentro del nombre. Esto sugiere otro agrupamiento. Por ejemplo, los módulos que inician con v.net. tienen que ver con análisis de redes vectoriales.

El nombre del módulo ayuda a entender su función, por ejemplo v.in.lidar inicia con v así que tiene que ver con mapas vectoriales, el nombre continúa con in lo que indica que el módulo es para importar datos a la Base de datos espacial de GRASS GIS y finalmente lidar indica que tiene que ver con nubes de puntos lidar.

Encontrando y corriendo un módulo

Para encontrar un módulo para su análisis, escriba el término en la caja de búsqueda en la pestaña Módulos (Modules) en el Administrador de capas (Layer manager), luego presione Enter hasta que encuentre el módulo.

De manera alternativa, puede buscar a través de el árbol de módulos en la pestaña Módulos (Modules). Puede también buscar en el menú principal. Por ejemplo, para encontrar información sobre un mapa ráster, use:Ráster → Reportes y estadísticas → Metadatos básicos del ráster (Raster → Reports and statistics → Basic raster metadata).

• 

  Búsqueda de un módulo en el árbol de módulos (busca en nombre, descripción, palabras clave)

• 

  Los módulos también están en el menú principal

Correr un módulo como un comando

Si ya sabe el nombre del módulo, puede simplemente usar la línea de comandos. La GUI ofrece una pestaña de consola de Comandos, con una línea de comandos específicamente construida para correr los módulos de GRASS GIS. Si escribe el nombre ahí, obtendrá sugerencias para completar automáticamete el nombre. Luego de presionar Enter, obtendrá una GUI para el módulo.

• 

  Sugerencia de autocompletado cuando se escribe el nombre del módulo: Al escribir el prefijo r. tenemos una lista de módulos que inician con el prefijo.

Puede usar la línea de comandos para correr todo el comando cuando ya lo tiene, esto es, módulo y lista de parámetros, en las instrucciones.

Línea de comandos vs. Interfaz GUI

Los módulos de GRASS pueden ser ejecutados tanto a través de la GUI o por la interfaz de línea de comandos. La GUI ofrece una aproximación más amigable para ejecutar los módulos donde el usuario puede navegar en las capas de datos cuando quiera analizar y modificar las opciones de procesamiento simplemente haciendo click en las cajas de la GUI. La GUI también ofrece un manual fácilmente accesible de cómo ejecutar el módulo. LA interfaz de línea de comandos permite ejecutar un módulo usando las opciones del comando específicas para el módulo. Esto es conveniente cuando se están corriendo análisis similares con pequeñas modificaciones o cuando se está familiarizado con los comandos de los módulos para un procesamiento eficiente. En este taller se dan los comandos, que pueden ser copiados y pegados en la línea de comando, aunque también se puede usar de manera paralela la GUI y la línea de comandos, dependiendo de la preferencia personal. Eche un vistazo cómo la GUI y la línea de comandos representan la misma herramienta.
Ejercicio: calcular el aspecto (orientación) dado un modelo digital de elevación usando el módulo r.slope.aspect usando tanto la GUI del módulo y la línea de comandos.

• ¿Cómo encontrar módulos? Los módulos están organizados por su funcionalidad en el menú de la wxGUI, o puede buscarlos en la pestaña de Buscar módulos. Si ya sabe qué modulo usar, puede simplemente escribirlo en la consola de comandos de la wxGUI.

Parámetros de módulo

• 

  Diálogo del módulo

El mismo análisis se puede hacer usando el siguiente comando:

r.neighbors -c input=elevation output=elev_smooth size=5

De manera alterna, puede rellenar parámetro por parámetro en la GUI, cuando ya tiene el comando.

Región computacional

Antes de usar el módulo para calcular un nuevo mapa ráster, debe definir correctamente una región computacional. Todos los cálculos que se hacen sobre rásters se realizan dentro de una extensión especificada y con una resolución dada.

La región computacional es un concepto de rásters importante dentro de GRASS GIS. En GRASS una región computacional puede ser definida dentro de una región ocupada por datos de mayor extensión para un análisis de prueba más rápido o por ejemplo dentro de regiones basadas en unidades administrativas. Proveemos algunos puntos para tener en cuenta cuando se use la función de región computacional:

• definida por extensión de la región y la resolución ráster
• aplica a todas las operaciones ráster
• persiste entre sesiones de GRASS, puede ser diferente para diferentes Directorios de mapas (mapsets)
• ventajas: mantiene la consistencia, evita el recorte, para tareas computacionalmente demandantes definir una región más pequeña, revisar si el resultado es bueno y luego definir la región computacional para toda el área de estudio y volver a correr el análisis.
• corra g.region -p, o en el menú Configuraciones - Región - Mostrar región para ver la configuración actual de la región.

• 

  Concepto de región computacional: un ráster con una extensión grade (azul) se muestra junto con otro más pequeño (verde). La región computacional (rojo) se define para que sea igual al ráster pequeño, así los cálculos estarán limitados a la región que ocupa el ráster menor aún cuando la entrada sea el ráster mayor. (No se muestra en la imagen: también la resolución, no solo la extensión, es igual a la resolución del ráster menor.)

• 

  Modos sencillos para definir la región computacional desde la GUI. A la izquierda, definir la región para que sea igual a un mapa ráster. A la derecha, seleccionar la opción señalada y luego defina la región dibujando un rectángulo.

• 

  Definir la región computacional (extensión y resolución) para que sea igual a un ráster (pestaña de Capas en el Administrador de capas)

Los valores numéricos de la región computacional pueden ser revisados usando:

g.region -p

Luego de ejecutar el comando, obtendrá algo como esto:

norte:      220750
sur:      220000
oeste:       638300
este:       639000
nsres:      1
eores:      1
rows:       750
columnas:       700
celdas:      525000

La región computacional también puede ser definida usando un mapa vectorial. En este caso, solamente se define la extensión (ya que el mapa vectorial no tiene resolución - al menos no en el sentido del mapa ráster). En la GUI, esto se puede hacer en la misma forma que con el mapa ráster. En la línea de comandos se ve así:

g.region vector=lakes

La resolución se puede definir de manera separada usando el parámetro res del módulo g.region. Las unidades son las mismas de la Localización actual, en nuestro caso metros. Esto se puede hacer en la pestaña Resolución del díalogo de g.region o en la línea de comandos del siguiente modo (usando también la bandera -p para mostrar los nuevos valores:

g.region res=3 -p

La nueva resolución puede modificarse ligeramente en este caso para ajustar a la región, que no estamos cambiando. Sin embargo, a menudo queremos que la resolución sea igual a los valores que proveemos y estamos de acuerdo con una pequeña modificación a la extensión. Para eso está la bandera -a.

El siguiente comando de ejemplo usará la extensión del vectorial llamado lakes, con una resolución 10, modificando la extensión para alinearla a que la resolución sea de 10 metros, y muestra los valores de la nueva región computacional:

g.region vector=lakes res=10 -a -p

Correr módulos

Encuentre el módulo para calcular pendiente y aspecto en el menú o en el árbol de módulo en Ráster → Análisis de terreno → Pendiente y aspecto o simplemente corra r.slope.aspect.

• 

  Seleccione el mapa ráster de entrada de elevación.

• 

  Ingrese los nombres de los mapas ráster de salida. También note que el comando correspondiente aparece en la parte inferior del diálogo GUI.

• 

  Maual de uso incluido en el diálogo GUI para ver más detalles.

• 

  Presione Ejecutar (1). Cuando se haya calculado, el resultado se añade al Administrador de capas y al Visualizador de mapas. Use Cerrar (2) para cerrar la ventana.

Vista 3D

Podemos explorar nuestra área de estudio en vista 3D.

1. Añada elev_lid792_1m y quite de la selección (deseleccione) o remueva cualquier otra capa.
2. Defina la región computacional a este ráster. Cambie a la vista 3D (en la esquina derecha del Visualizador de mapas).
3. Ajuste la vista (perspectiva, altura, exageración vertical)
4. En la pestaña Datos, defina

In Data tab, set Resolución modo Fino a 1 y defina slope (calculado en el paso anterior) como el color de la superficie.

1. Cuando termine, cambie de nuevo a la vista 2D.

Vista 3D de el MDE elev_lid792_1m con slope como color

Análisis de rásters y vectoriales

Distancia desde el borde del bosque

Use álgebra de mapas para extraer una clase de bosque dada (en este caso la 5) del ráster de clasificación de uso de suelo:

r.mapcalc "forest = if(landclass96 == 5, 1, null())"

La función if() que se usó tiene tres parámetros, con la siguiente sintáxis:

if(condición, valor si verdadero, valor si falso)

Después usamos el operador == que evalúa como verdadero solo cuando ambos lados son iguales. Finalmente usamos la función null() que representa valor nulo NULL (sin dato).

Ahora podemos obtener la distancia al borde del bosque usando el módulo r.grow.distance, que calcula distancias a áreas con valores en áreas sin valores (con NULL) o al contrario. De manera predeterminada nos da la distancia al borde del bosque desde fuera del bosque, pero ahora usaremos la bandera -n para obtener la distancia al borde desde dentro del mismo bosque:

r.grow.distance -n input=forest distance=distance

• 

  Definiendo el color verde para el mapa ráster de bosque (Click derecho en la capa de mapa ráster, seleccione Establecer paleta de colores interactivamente).

• 

  Diálogo de r.grow.distance y las distancias calculadas como fondo.

Etadísticas de puntos (Point statistics)

Importar un Archivo Shapefile

Descargue el archivo Shape de muestra points_of_interest.zip y descomprímalo.

Importe el archivo usando el módulo v.in.ogr. Note que necesita especificar la ruta de búsqueda completa.

v.in.ogr input=/ruta/a/points_of_interest.shp output=points_of_interest

Generando una malla hexagonal

Para calcular la densidad de puntos en una malla hexagonal a partir del mapa vectorial de points_of_interest use el mapa vectorial mismo para definir la extesión de la región computacional. La resolución está basada en el tamaño deseado de los hexágonos.

g.region vector=points_of_interest res=2000 -pa

Aunque la Región computacional generalmente no es usada para el procesamiento de vectoriales, la malla hexagonal es creada como un mapa vectorial basándose en la extensión previamente seleccionada y en el tamaño de la malla.

v.mkgrid map=hexagons -h

Calculando estadísticas de puntos en polígonos

Lo que sigue es el conteo del número de puntos por hexágonos usando el módulo v.vect.stats.

v.vect.stats points=points_of_interest areas=hexagons count_column=count

Por último, el siguiente comando define la paleta de colores a viridis basada en la columna del conteo. Use la paleta de colores ryb si tiene GRASS GIS 7.0.

v.colors map=hexagons use=attr column=count color=viridis

• 

  Hexágonos coloreados y módulos usados para crearlos.

Análisis de estructura de paisaje

Para el siguiente ejemplo usaremos el ráster landuse96_28m de nuestro conjunto de datos de North Carolina, donde los parches representan diferentes coverturas de suelo (usos). Usamos el módulo Vamos a usa el módulo r.neighbors y las extensiones (addons) r.diversity y r.forestfrag. Primero instalamos las extensiones:

 g.extension r.diversity
 g.extension r.forestfrag

Cerramos la GUI de GRASS y la abrimos de nuevo.

Riqueza Primero calculamos la riqueza (cantidad de clases únicas o distintas), usando r.neighbors. Al usar una ventana móvil, creamos un ráster que representa la variabilidad espacial de la riqueza. El tamaño de la ventana está en celdas.

 g.region raster=landuse96_28m
 r.neighbors input=landuse96_28m output=richness method=diversity size=15

Índices de paisaje

La extensión r.diversity calcula diversos índices de paisaje usando una ventana móvil. Está basado en los módulos r.li para análisis de estructura de paisaje. En este ejemplo calculamos los índices de diversidad de Simpson, Shannon y Renyi con un rango de ventanas móviles:

 r.diversity input=landuse96_28m prefix=index alpha=0.8 size=9-21 method=simpson,shannon,renyi

Esto genera 9 rásters con nombres como index_simpson_size_9. Podemos añadirlos al Visualizador de mapas usando Añadir varios ráster o vectoriales' en el Administrador de capas, en el menú Archivo.

Para verlos, necesitamos definir la misma paleta de colores para mapas ráster del mismo índice.

 r.colors map=index_shannon_size_21,index_shannon_size_15,index_shannon_size_9 color=viridis
 r.colors map=index_renyi_size_21_alpha_0.8,index_renyi_size_15_alpha_0.8,index_renyi_size_15_alpha_0.8 color=viridis
 # usamos la paleta de colores grey1.0 por que Simpson es de 0 a 1
 r.colors map=index_simpson_size_21,index_simpson_size_15,index_simpson_size_9 color=grey1.0

Fragmentación del bosque

Fragmentación del bosque calculado con la extensión r.forestfrag

El módulo r.forestfrag calcula la fragmentación del bosque siguiendo una metodología propuesta por Riitters et al. (2000).

Primero se marcan todas las celdas que tienen bosque como 1 y todo lo demás como 0:

 # Definir la región
 g.region raster=landclass96
 # listar las clases:
 r.category map=landclass96
 # seleccionar clases
 r.mapcalc "forest = if(landclass96 == 5, 1, 0)"

Use el nuevo mapa ráster de presencia de bosque para calcular el índice de fragmentación del bosque con una ventana de tamaño 15:

 r.forestfrag input=forest output=fragmentation window=15

Reporte de la distribución de categorías de fragmentación:

 r.report map=fragmentation units=k,p

Realizando scripts con Python

El modo más sencillo de ejecutar el código de Python que use paquetes de GRASS GIS es simplemente usar el Editor Simple de Python integrado a GRASS GIS accesible desde la barra de herramientas o la pestaña Python en el Administrador de capas. Otra opción es escribir código de Python en su editor de texto favorito, como Notepad++ (note que los editores de Python son editores de texto simple). Luego corra el script en GRASS GIS usando el menú principal Archivo -> Lanzar script.

• 

  Editor simple de Python integrado a GRASS GIS (desde la versión 7.2) con la pestaña Python en el fondo y que contiene una consola de Python interactiva.

• 

  Pestaña Python con una consola de Python interactiva.

La GRASS GIS 7 Pytho Scripting Library provee funciones para llamar los módulos de GRASS dentro de scripts como subprocesos. Las funciones más usadas incluyen:

• run_command: usualmente usado con módulos cuya salida son datos ráster/vectorial y no se espera un texto como salida.
• read_command: usado cuando se tiene interés en el texto de salida
• parse_command: usado con módulos que producen texto como salida en pares clave=valor.
• write_command: usado con módulos que se espera un texto de entrada, ya sea de la entrada estándar o un archivo.

Además de esto, esta librería provee varias funciones simplificadas (wrapper functions) para módulos que se llaman frecuentemente.

Llamando módulos GRASS GIS

Usaremos la Consola de Python de GRASS GIS para correr los comandos. Para scripts más largos, puede crear un archivo de texto, guardarlo en el directorio de trabajo actual y correrlo con python myscript.py a partir de la consola de comandos GUI o de una terminal.

Tip: Cuando copie fragmentos de código a la consola GUI de Python, de click derecho en la posición y seleccione Paste Plus en el menú de contexto. De otro modo no va a funcionar para fragmetos de código de varias líneas.

Empezamos por importar la librería GRASS GIS Python Scriptig:

import grass.script as gscript

Antes de correr cualquier módulo de GRASS, necesita definir la región computacional usando g.region. En este ejemplo, definimos la extensión de la región computacional y la resolución a la capa ráster elevation.

gscript.run_command('g.region', raster='elevation')

La función run_command() es la más frecuentemente usada. Aquí aplicamos la operación de promedio focal (r.neighbors) para suavizar la capa ráster de elevación. Note que la sintáxis es similar al usado en la consola, solo que las banderas se especifican como parámetros.

gscript.run_command('r.neighbors', input='elevation', output='elev_smoothed', method='average', flags='c')

Si corremos los comandos de Python desde la consola GUI de Python, podemos usar AddLayer para añadir las capas recién creadas:

AddLayer('elev_smoothed')

Llamado módulos de GRASS GIS con entrada o salida de texto

La salida de texto de los módulos puede ser capturada usando la función read_command().

gscript.read_command('g.region', flags='p')

gscript.read_command('r.univar', map='elev_smoothed', flags='g')

Algunos módulos pueden producir una salida en formato de valores con clave (key-value), esto se realiza con la bandera -g. La función parse_command() automaticamente analiza esta salida y regresa un diccionario. En este ejemplo llamamos g.proj para mostrar los parámetros de la proyección de la Localización actual.

gscript.parse_command('g.proj', flags='g')

Como comparación, abajo el mismo ejemplo, pero usando la función read_commad():

gscript.read_command('g.proj', flags='g')

Algunos módulos requieren que el texto de entrada esté en un archivo o sea proveído en la entrada estándar. Usando la función write_command() podemos pasar el texto al módulo. Aquí, estamos creando un nuevo vectorial con un punto con v.in.ascii. Note que el parámetro stdin no es usado como un parámetro de módulo, pero su contenido es pasado al subproceso como la entrada estándar.

gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (635818, 221342), output='point')

Si corremos los comandos Python desde la consola GUI de Python, podemos usar AddLayer para añadir la capa recién creada:

AddLayer('point')

Funciones simplificadas (wrapper functions)

Algunos módulos tienen funciones simplificadas (wrapper functions) para simplificar las tareas de uso común. Por ejemplo podemos obtener la información sobre una capa ráster con raster_info que es la función simplificada de r.info, o una capa vectorial con vector_info.

gscript.raster_info('elevation')
gscript.vector_info('point')

Otro ejemplo es usando la simplificada de r.mapcalc para álgebra ráster:

gscript.mapcalc("elev_strip = if(elevation > 100 && elevation < 125, elevation, null())")
gscript.read_command('r.univar', map='elev_strip', flags='g')

La función region es un modo conveniente de obtener la configuración de la región actual. Regresa un diccionario con valores convertido a los tipos apropiados (flotantes y enteros).

region = gscript.region()
print region
# área de celdas en unidades de mapa (en Localizaciones proyectadas)
region['nsres'] * region['ewres']

Podemos enlistar los datos guardados en una Localización de GRASS GIS con la simplificada de g.list. Con list_grouped, las capas de mapas son agrupadas por Directorios de mapas (capas ráster en este ejemplo):

gscript.list_grouped(type=['raster'])
gscript.list_grouped(type=['raster'], pattern="landuse*")

Aquí hay un ejemplo de un wrapper diferente para g.list list_pairs que estrucutra la salida como listas de pares (nombre, Directorio de mapas). Obtenemos el Directorio de mapas actual con la wrapper de g.gisev.

current_mapset = gscript.gisenv()['MAPSET']
gscript.list_pairs('raster', mapset=current_mapset)

Ejercicio

Usando "elev_*" como prefijo en el nombre, exportar todas las capas ráster desde su Directorio de mapas como GeoTiff (vea r.out.gdal. No olvide definir la región actual g.region) para cada mapa para que sea igual a la extensión y resolución de cada uno, ya que pueden variar.

Creando una Localización GRASS GIS nueva e importar datos

Para el siguiente ejemplo con R, necesitamos crear primero una nueva Localización GRASS e importar datos. Vamos a usar los datos de temperatura PRISM y los límites vectoriales de los estados de EEUU. Vamos a crear una nueva Localización basada en el código EPSG 4269 (datum NAD83).

1. Inicie GRASS GIS
2. Seleccione Nueva en la parte izquierda de la pantalla de bienvenida para iniciar el Ayudante de Localizaciones
3. En el ayudante, ingrese el nombre de la nueva Localización, por ejemplo PRISM, presione Siguiente
4. Seleccione Seleccionar código EPSG del sistema de referencia espacial, presione Siguiente
5. Ingrese el número 4269 en el campo de Código EPSG, presione Siguiente
6. Seleccione 1 en el diálogo de transformación del Datum
7. Revise las definiciones PROJ.4 y presione "Finalizar"

Se crea automáticamente el Directorio de mapas PERMANENT, así que puede iniciar su sesión GRASS.

Descargue y descomprima el siguiente conjunto de datos en una carpeta:

• US Census Bureau: Límites cartográficos de los estados de EEUU
• temperatura media anual normal de 30 años PRISM, malla de 4km
• datos de elevación SRTM para los EEUU

Definir el Directorio de trabajo actual en GRASS GIS para hacer más fácil el encontrar los archivos extraidos. En el menú Configuraciones - Entorno de trabajo de GRASS - Cambiar de directorio de trabajo defina el directorio en el cual tenga los datos. De manera alternativa en la consola de comandos de la GUI escriba cd y seleccione el directorio.

Ahora podemos importarlo a GRASS GIS:

 r.import input=PRISM_tmean_30yr_normal_4kmM2_annual_bil.bil output=temp_mean
 v.import input=cb_2015_us_state_20m.shp output=boundaries
 r.import input=US_elevation.tif output=elevation

No especificamos la ruta completa a los archivos gracias a que definimos más arriba el directorio de trabajo. Si usamos el diálogo, podemos explorar los archivos usando el botón Explorar.

Ahora visualizamos los datos y definimos la rampa de colores apropiada para la temperatura:

 r.colors map=temp_mean@PERMANENT color=celsius

Realizando scripts con R

Usar R con GRASS GIS se puede hacer de dos maneras:

• Usar R dentro de una sesión de GRASS GIS, i.e. puede iniciar R (o RStudio) desde la línea de comandos de GRASS GIS.
  • Trabaja con datos en la Base de datos espaciales GRASS GIS usando GRASS GIS
  • No inicia la función initGRASS() (GRASS GIS ya está corriendo).
• Usando GRASS GIS dentro de una sesión de R, i.e. se conecta a la Base de datos espaciales de GRASS GIS desde R (o RStudio).
  • Pone datos en la Base de datos GRASS GIS solamente para realizar las operaciones de GRASS GIS.
  • Use la función initGRASS() para iniciar GRASS GIS dentro de R.

Vamos a correr R dentro de una sesión de GRASS GIS (la primer manera). Lanzar R dentro de GRASS GIS e instalar los paquetes rgrass7 y rgdal

install.packages("rgrass7")
install.packages("rgdal")
library("rgrass7")
library("rgdal")

Podemos ejecutar módulos GRASS usando la función execGRASS:

execGRASS("g.region", raster="temp_mean", flags="p")

Vamos a analizar la relación entre la temperatura, la elevación y la latitud.

Primero vamos a generar un ráster con los valores de latitud:

execGRASS("r.mapcalc", expression="latitude = y()")

Nota: en un sistema de coordenadas proyectada puede usar r.latlong.

Después, vamos a generar puntos aleatorios y muestear el conjunto de datos.

execGRASS("v.random", output="samples", npoints=1000)
# esto va a restrigir el muestreo dentro de los límites de EEUU
# estamos sobreescribiendo muestras vectoriales, así que necesitamos usar la bandera overwrite (sobreescribir)
execGRASS("v.random", output="samples", npoints=1000, restrict="boundaries", flags=c("overwrite"))
# creamos la tabla de atributos
execGRASS("v.db.addtable", map="samples", columns=c("elevation double precision", "latitude double precision", "temp double precision"))
# muestreamos rásters individuales
execGRASS("v.what.rast", map="samples", raster="temp_mean", column="temp")
execGRASS("v.what.rast", map="samples", raster="latitude", column="latitude")
execGRASS("v.what.rast", map="samples", raster="elevation", column="elevation")

Nota: hay un módulo de extensión (addon) r.sample.category que simplifica este proceso y muestrea múltiples rásters.

Ahora abrimos el administrador de tablas de atributos de GRASS GIS para inspeccionar los valores muestreados o usar v.db.select para enlistar los valores. Exploramos el conjunto de datos en R:

samples <- readVECT("samples")
summary(samples)
plot(samples@data)

Calculamos el modelo lineal multivariado:

linmodel <- lm(temp ~ elevation + latitude, samples)
summary(linmodel)

Predecimos la temperatura usando este modelo:

maps <- readRAST(c("elevation", "latitude"))
maps$temp_model <- predict(linmodel, newdata=maps)
spplot(maps, "temp_model")
# write modeled temperature to GRASS raster and set color ramp
writeRAST(maps, "temp_model", zcol="temp_model")
execGRASS("r.colors", map="temp_model", color="celsius")

Comparamos el modelo lineal simple con los datos reales:

execGRASS("r.mapcalc", expression="diff = temp_mean - temp_model")
execGRASS("r.colors", map="diff", color="differences")

En la GUI de GRASS, añadimos las capas temp_mean y temp_model, las seleccionamos y vamos a Archivo - Ventana de visualización móvil para comparar visualmente el modelo y la temperatura real.

• 

  Comparasión de la temperatura media anual PRISM y la temperatura modelada basada en latitud y elevación. Izquierda: diferencias entre el modelado y el real, en grados Celsius. Derecha: Usando el Visualizador de ventana móvil para evaluar visualmente el modelo (temperatura modelada del lado derecho).

Manipulación y visualización de datos Espacio-temporales

GRASS GIS7 tiene una nueva librería y una serie de módulos para el manejo y análisis de datos espacio-temporales. El sistema (framework) temporal de GRASS no trabaja con capas de mapas individuales, sino con conjuntos de datos espacio-temporales. Un conjunto de datos es una colección de capas de mapas individuales con una marca de tiempo asignada. Todos los datos espaciales son guardados en la base de datos GRASS estándar, el sistema temporal administra el metadato temporal en una base de datos temporales separada.

Información general en el manual: https://grass.osgeo.org/grass72/manuals/temporalintro.html

Conceptos y terminología importante:

• Conjunto de datos ráster espacio temporales (strds) son diseñados para manejar series de tiempos de mapas ráster. Los módulos que procesan strds tienen el prefijo t.rast.
• Conjunto de datos ráster espacio temporales 3D (str3ds) son diseñados para manejar series de tiempo de mapas ráster 3D. Los módulos que procesan str3ds tienen el prefijo t.rast3d.
• Conjuntos de datos vectoriales espacio temporales (stvds) son diseñados para manejar series de tiempo de mapas vectoriales. Los módulos que procesan stvds tienen prefijo t.vect.
• Los Conjuntos de datos pueden usar intervals (intervalos) (definidos por marca de tiempo de inicio y término) o instances (solo tiempo de inicio)
• Los conjuntos de datos pueden usar tiempos absolute (absolutos) (marcas de tiempo como 2017-04-06 22:39:49) o relative (relativos) (por ej. (e.g., 4 years, -90 days).
• La granularity (granularidad) es el divisor común menor de la extensión temporal (y posibles vacíos (gaps)) de todos los mapas del conjunto de datos.
• La topology (topología) temporal analiza las relaciones temporales entre intervalos de tiempo.
• El sampling (muestreo) temporal es usado para determinar el estado de un proceso durante un segundo proceso.

Vista general del flujo de trabajo:

1. Crear un conjunto de datos vacío: strds, str3ds o stvds (t.create),
2. registrar los mapas GRASS GIS (t.register),
3. revisar el conjunto de datos espacio temporales creado (t.list), t.info, t.rast.list)
4. hacer el análisis t.rast.aggregate, t.info, t.rast.univar, t.vect.univar, ... ¡y muchos más!

Ejemplo de modelación de radiación solar

Vamos a generar datos espacio-temporales y luego explorarlos usando los módulos temporales y visualizarlos como una animación.

Vamos a usar r.sun.hourly para calcular una serie de irradiancias solares durante un día en un modelo digital de superficie (incluye copas de los árboles) derivado previamente en la parte de lidar de este taller.

Descargue la extensión (addon) si no lo ha hecho previamente:

g.extension r.sun.hourly

Convierta la fecha del fía de hoy (o cualquier fecha) a día del año corriendo este comando en la pestaña de la consola de Python en la GUI.

 from datetime import datetime
 datetime.now().timetuple().tm_yday
 # or for an arbitrary day:
 datetime.datetime(2017, 6, 21).timetuple().tm_yday

Use este número para la opción día. Calcular la series ráster de la irradiancia del haz (sea paciente) con el siguiente comando. La serie de tiempo es registrada automáticamente en un conjunto de datos ráster espacio temporales.

 # first set region
 g.region raster=el_D783_6m -p
 r.sun.hourly -t elevation=el_D783_6m start_time=6 end_time=20 day=101 year=2017 time_step=0.5 beam_rad_basename=solar nprocs=4

Ahora cerciorémonos que el conjunto de datos temporales haya sido creado:

 t.list type=strds

Enliste los metadatos y la información analítica sobre el conjunto de datos:

 t.info input=solar type=strds

Muestre cuales capas ráster están registradas:

 t.rast.list input=solar

Defina la paleta de colores para el conjunto de datos creado. Use la herramienta interactiva en el díalogo de t.rast.colors (pestaña Definir, bajo reglas) o cree un archivo de texto con la siguiente regla de colores y páselo a la opción reglas:

 0% 60:60:60
 70% yellow
 100% 255:70:0

 t.rast.colors input=solar rules=rules.txt

Ahora visualizaremos los resultados como una animación.

Inicie la herramienta de animaciones de GRASS GIS desde Archivo - Herramienta de animación y cree una nueva animación. Seleccione conjunto de datos ráster espacio temporales, vea la siguiente figura para más detalles.

• 

  Usando la Herramienta de animación de GRASS GIS para animar la irradiancia solar durante un día. Calculado con r.sun.hourly.

Procesamiento de datos Lidar

Vamos a importar una nube de puntos lidar como un archivo LAS y usar binning (use binning) para analizar nuestros datos, y luego interpolarlos para crear un modelo digital de superficie. Vamos a usar r.in.lidar y v.in.lidar para trabajar con archivos LAS. En caso que no esté disponible en su instalación GRASS (GRASS necesita estar compilado con libLAS), use r.in.xyz y v.in.ascii con los archivos de texto correspondientes, el flujo de trabajo es muy similar.

Datos necesarios:

• mosaico LAS para Raleigh
• versión ttxt para el mosaico (todos los puntos) and y solamente los primeros puntos de retorno

Análisis Lidar usando binning

Calcular la densidad de la nube de puntos:

 g.region  n=224034 s=223266 e=640086 w=639318 res=1
 r.in.lidar -o input=tile_0793_016_spm.las method=n output=tile_density_1m resolution=1
 r.in.lidar -o input=tile_0793_016_spm.las method=n output=tile_density_6m resolution=6

Podemos usar el binning para calcular el rango de los valores de elevación y otras estadísticas:

 r.in.lidar -o input=tile_0793_016_spm.las method=range output=tile_range_6m resolution=6
 # el mapa y las siguientes estadísticas muestran que hay un punto discrepante alrededor de los 600 m
 r.report map=tile_range_6m units=c
 # cambiar la paleta de colores al histograma ecualizado para ver el constraste
 r.colors map=tile_range_6m color=viridis -e

Podemos flitrar los datos por clase o retorno:

 r.in.lidar -o input=tile_0793_016_spm.las method=mean output=tile_ground_6m resolution=6 class_filter=2
 r.colors map=tile_ground_6m color=elevation

Interpolación

Primero importamos los archivos LAS como puntos vectoriales, nos quedamos solamente con los primeros puntos de retorno y limitamos la importación verticalmente para impedir usar los puntos discrepantes encontrados en el paso anterior.

 v.in.lidar -obt input=tile_0793_016_spm.las  output=tile_points return_filter=first zrange=0,300

Después interpolamos:

 g.region -a -p vector=tile_points res=2
 v.surf.rst input=tile_points elevation=tile_dsm tension=25 smooth=1 npmin=100

Puede visualizar el DSM usando la vista 3D (primero deseleccione todas las demás capas y luego cambie a 3D) o calculando el sombreado de relieve con r.relief.

Ver también

• GRASS and Python
• GRASS Python Scripting Library
• R and the rgrass7 package
• Processing raster time series with GRASS GIS and R