Álgebra 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 realiza álgebra de mapas R.
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 calcular la zona inundable del estado de Carolina del Norte, generando mapas de cotas, pendientes, curvas de nivel y relieves. Queremos visualizar en R los siguientes imágen mediante la interfaz GRASS-R.

mapa final

cota690_state

orientacion_state

inundaciones_state


pendientes025_state


pendientes_state

state_shaded

#!/bin/sh

########### PRÁCTICA 3. ÁLGEBRA DE MAPAS ##################
 ##Herramientas de GRASS que utilizan álgebra de mapas para realizar cáculos complejos.
 ###Objetivo general: Calcularemos la zona inundable del estado de Carolina del Norte, generando en nuestro MAPSET los mapas correspondientes, así como cotas, pendientes, curvas de nivel y relieves.
 #Partiendo del mapa de elevaciones de todo el estado   (elev_state_500m@PERMANENT) vamos a calcular unas cuantas características, utilizando álgebra de mapas tanto directamente (como hicímos en la práctica anterior) como indirectamente (a través de herramientas de GRASS que utilizan el algebra).
 ###Objetivos específicos:
   #  1) Obtener el mapa de relieve del estado y visualizarlo.
   #  2) Obtener el mapa de curvas de nivel.
   #  3) Obtener el mapa de pendientes y el mapa de orientación del estado.
   #  4) Calcular el mapa de zonas inundables (necesitaremos calcular la cota y las pendientes que afectan a estas zonas previamente)

#### 0) Seleccionar la región activa (depende de los mapas con que queremos trabajar)
 d.mon x0
 g.region rast=elev_state_500m@PERMANENT
 d.rast map=elev_state_500m@PERMANENT

#### 1) Calcular el mapa de relieve del estado y visualizarlo.
 #Utilizamos la orden r.shaded.relief para generar asi un nuevo mapa, state_shaded, suponiendo que el sol esta en la posicion dada por ese acimut y altura.
 r.shaded.relief map=elev_state_500m@PERMANENT shadedmap=state_shaded altitude=30 azimuth=270 zmult=1 scale=1
 #Visualizamos el mapa utilizando la opacidad de la capa raster
 d.his h_map=elev_state_500m@PERMANENT i_map=state_shaded

#### 2) Obtener el mapa vectorial de curvas de nivel para el estado de Carolina del Norte.
 #Aunque vamos a obtener un mapa vectorial, los caculos necesarios son caculos raster. Vamos a utilizar la orden r.contour y le pediremos que nos extraiga las curvas de nivel cada 50m (step=50).
 r.contour input=elev_state_500m@PERMANENT output=curvas_nivel_state step=50 cut=0

#### 3) Obtener el mapa de pendientes y el mapa de orientacion del estado.
#Utilizamos la orden r.slope.aspect. Construiremos dos mapas: pendientes_state y orientacion_state. Tenemos con ello las pendientes del terreno y hacia donde estan orientadas.
 r.slope.aspect elevation=elev_state_500m@PERMANENT slope=pendientes_state format=degrees prec=float aspect=orientacion_state zfactor=1.0 min_slp_allowed=0.0

#### 4) Calcular el mapa de zonas inundables:  Consideraremos la hipótesis de que las zonas inundables del estado de Carolina del norte están por debajo de la cota 690 y se presentan en terrenos con pendientes inferiores al 2.5%.
 ### a) Calculamos el mapa de cota 690 (puntos del estado de Carolina del Norte que alcanzan y están por debajo de esa altura -690m-) que afectan la zona de estudio: Tendremos dos categorías en nuestro mapa raster (los puntos que pertenecen al conjunto de puntos que cumplen esa característica y los que no): 1 y 0
 r.mapcalculator amap=elev_state_500m@PERMANENT formula="if(A<691,1,0)" outfile=cota690_state_mapcalc

 ### b) Calculamos el mapa de pendientes que afectan la zona de estudio: Hacemos lo mismo que con el mapa de cotas
 r.mapcalculator amap=pendientes_state formula="if(A<2.5,1, 0)" outfile=pendientes025_state_mapcalc

 ### c) Ahora, para ver las zonas inundables por debajo de esa cota para esas pendientes, sólo tenemos que multiplicar los dos mapas obtenidos, cota690_state_mapcalc y pendientes025_state_mapcalc, con álgebra de mapas:
  r.mapcalculator amap=cota690_state_mapcalc bmap=pendientes025_state_mapcalc formula="A*B" outfile=inundaciones_state
  r.colors map=inundaciones_state color=byg #Le cambiamos el color
  d.rast  inundaciones_state #Donde el valor 1 (la zona verde) corresponde a las zonas inundables que buscábamos. 
  d.his h_map=inundaciones_state i_map=state_shaded  #Lo superponemos al mapa de relieves, para comprobarlo visualmente
  d.legend map=inundaciones_state color=black  lines=0 thin=1  labelnum=5 at=75,95,5,10 #agregamos una leyenda