Consulta y transformación de mapas. Interfaz R-GRASS

Esta es una práctica dictada en el curso "Análisis espacial con sistemas GIS (1ª Edición)" de la Universidad de Granada (UGR), en el cual se realizan los pasos básicos de consulta y transformación de mapas.

Para este ejercicio debemos instalar GRASS (versiones 6.2, 6.3 o 6.4) y R (la versión que deseen), así como el paquete spgrass6 de R.

Ejemplo en R: queremos visualizar en R la siguiente imágen mediante la interfaz GRASS-R.

#!/bin/sh

#ACTIVIDAD 2: Queremos obtener el mapa de elevaciones de un condado concreto del estado de North Carolina, simplificarlo y que sólo muestre cuatro categorías: valles, estribaciones, baja montaña y alta montaña. Los pasos a seguir son:

#• Consultar el mapa vectorial del estado de North Carolina: ver qué información se guarda de cada condado y localizar el condado cuyo nombre es Wilkes.
 #Cargar y visualizar los mapas: Mapa vectorial de condados del estado de North Carolina: boundary_county@PERMANENT; Mapa raster de elevaciones del estado de North Carolina: elev_state_500m@PERMANENT
 d.mon x1
 d.rast map=elev_state_500m@PERMANENT -o
 d.vect map=boundary_county@PERMANENT -c color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1
 #Consultar el mapa vectorial: mostrar columnas de atributos (v.info)
 v.info -c map=boundary_county@PERMANENT layer=1
 #Consultar interactivamente un condado al azar (v.what) y para esa coordenada obtener el valor de altura (r.what)
 v.what -a map=boundary_county@PERMANENT east_north=632337.807606,221543.624161 distance=18232.6621924
 ##Salida
 #East: 632337.807606
 #North: 221543.624161
 #Map: boundary_county 
 #Mapset: PERMANENT
 #Centroid
 #sigue...

 g.region rast=elev_state_500m@PERMANENT
 r.what input=elev_state_500m@PERMANENT cache=500 null=* east_north=632337.807606,221543.624161
 ##Salida: 632337.807606|221543.624161||126.13894254

 #Buscar el condado de Wilkes utilizando el atributo de área que contiene su nombre: (“consulta SQL” en Administrador)
 d.vect map=boundary_county@PERMANENT -c color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1 where="NAME='WILKES'"

#• Crear un mapa vectorial del condado de Wilkes, extrayéndolo del mapa vectorial del estado de North Carolina.
 #Crear mapa vectorial (“guardar objetos visualizados en un nuevo vectorial” en Administrador y darle un nombre “wilkes_country” que se guardará en el MAPSET; visualizar y cargar el nuevo mapa)
 v.extract input=boundary_county@PERMANENT output=WILKES_country type=point,line,boundary,centroid,area layer=1 where="NAME='WILKES'"

#• Convertir ese mapa vectorial del condado de Wilkes en un mapa raster.
 #Zoom al mapa seleccionado
 d.vect map=WILKES_country color=0:0:0 lcolor=0:0:0 fcolor=170:170:170 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1
 #Cambiar la región activa a la región del mapa vectorial que hemos creado (display: set current región to match display)
 g.region vect=WILKES_country
 #Convertir vector en raster (v.to.rast uss=”atr” type=”point, line, area” layer=1 column=”cat” value=1 rows=4096 labelcolumn=”NAME”)
 v.to.rast input=WILKES_country output=WILKES_country use=attr column=cat layer=1 value=1 rows=4096
 #verificar con : r.info map=WILKES_country
 d.rast map=WILKES_country

#• Calcular la zona del mapa de elevaciones raster del estado de North Carolina que se corresponde con el condado de Wilkes y crear el mapa raster de elevaciones del condado de Wilkes.
 #Algebra de datos: 
 r.mapcalculator amap=WILKES_country bmap=elev_state_500m@PERMANENT formula="if(isnull(A),A,B)" outfile=elev_WILKES_country
#• Añadir una leyenda a este último mapa y visualizarlo.
 d.rast map=elev_WILKES_country
 d.legend map=elev_WILKES_country color=black  lines=0 thin=1  labelnum=5 at=5,95,5,10
#• Reclasificar ese último mapa en las cuatro categorias anteriores, según la altura: (r.reclass) p10
 #• Valles, si el punto (pixel) tiene un valor entre 271 y 400 metros.
 #• Estribaciones, si el pixel tiene un valor entre 400 y 700 metros.
 #• Baja montaña, si el pixel tiene un valor entre 700 y 1000 metros.
 #• Alta montaña, si el pixel tiene un valor entre 1000 y 1168 metros.
 r.reclass input=elev_WILKES_country output=elev_WILKES_country_reclass < "/home/usuario/CursoGIS/actividad2/rules"

# Visualizar lo anterior y añadirle una leyenda con los nombres de las categorias.
 d.rast map=elev_WILKES_country_reclass
 d.legend map=elev_WILKES_country_reclass color=black  lines=0 thin=1  labelnum=5 at=80,95,2,10