Source: https://grasswiki.osgeo.org/w/index.php?title=Unleash_the_power_of_GRASS_GIS_at_US-IALE_2017/es&direction=prev&oldid=24660
Timestamp: 2020-08-12 04:24:42+00:00

Document:
Unleash the power of GRASS GIS at US-IALE 2017/es - GRASS-Wiki
Unleash the power of GRASS GIS at US-IALE 2017/es
< Unleash the power of GRASS GIS at US-IALE 2017
Revision as of 14:53, 22 June 2017 by Dat (talk | contribs)
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.
1.1 Extensiones de GRASS GIS (addons)
1.2 Paquetes de R
2 Datos para este taller
3 Introducción a GRASS GIS
3.1 Estructura de la base de datos espacial de GRASS GIS
3.2 Creando una base de datos de GRASS para el tutorial
3.3 Mostrar y explorar los datos
3.4 Módulos de GRASS GIS
3.4.1 Encontrando y corriendo un módulo
3.4.2 Correr un módulo como un comando
3.4.3 Línea de comandos vs. Interfaz GUI
3.4.4 Parámetros de módulo
3.5 Región computacional
3.6 Correr módulos
3.7 Vista 3D
4 Análisis de rásters y vectoriales
4.1 Distancia desde el borde del bosque
4.2 Etadísticas de putos (Point statistics)
4.2.1 Importar un Archivo Shapefile
4.2.2 Generando una malla hexagonal
4.2.3 Calculando estadísticas de puntos en polígonos
5 Análisis de estructura de paisaje
6 Realizando scripts con Python
6.1 Llamando módulos GRASS GIS
6.2 Llamado módulos de GRASS GIS con entrada o salida de texto
6.3 Funciones simplificadas (wrapper functions)
6.4 Ejercicio
7 Creando una Localización GRASS GIS nueva e importar datos
8 Realizando scripts con R
9 Manipulación y visualización de datos Espacio-temporales
9.1 Ejemplo de modelación de radiación solar
10 Procesamiento de datos Lidar
10.1 Análisis Lidar usando binning
10.2 Interpolación
Se usará GRASS GIS 7.2 y R (>=3.1).
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:
Instale GRASS GIS desde el administrador de paquetes:
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:
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:
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.
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
Encontrando y corriendo un módulo
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.
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:
De manera alterna, puede rellenar parámetro por parámetro en la GUI, cuando ya tiene el comando.
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
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.
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í:
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:
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:
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.
Podemos explorar nuestra área de estudio en vista 3D.
Añada elev_lid792_1m y quite de la selección (deseleccione) o remueva cualquier otra capa.
Defina la región computacional a este ráster. Cambie a la vista 3D (en la esquina derecha del Visualizador de mapas).
Ajuste la vista (perspectiva, altura, exageración vertical)
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.
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 putos (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.
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.
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:
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:
# 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:
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.
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:
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().
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.
Como comparación, abajo el mismo ejemplo, pero usando la función read_commad():
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:
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).
# área de celdas en unidades de mapa (en Localizaciones proyectadas)
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):
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.
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).
Inicie GRASS GIS
Seleccione Nueva en la parte izquierda de la pantalla de bienvenida para iniciar el Ayudante de Localizaciones
En el ayudante, ingrese el nombre de la nueva Localización, por ejemplo PRISM, presione Siguiente
Seleccione Seleccionar código EPSG del sistema de referencia espacial, presione Siguiente
Ingrese el número 4269 en el campo de Código EPSG, presione Siguiente
Seleccione 1 en el diálogo de transformación del Datum
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:
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")
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")
plot(samples@data)
Calculamos el modelo lineal multivariado:
linmodel <- lm(temp ~ elevation + latitude, samples)
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:
Crear un conjunto de datos vacío: strds, str3ds o stvds (t.create),
registrar los mapas GRASS GIS (t.register),
revisar el conjunto de datos espacio temporales creado (t.list), t.info, t.rast.list)
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.
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:
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
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.
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.
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
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.
R and the rgrass7 package
Processing raster time series with GRASS GIS and R
Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=Unleash_the_power_of_GRASS_GIS_at_US-IALE_2017/es&oldid=24286"

References: resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 resolución 
 resolución 
 resolución