Superficie de coste. 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 en R para calcular una superficie de coste.
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 el camino más corto al parque de bomberos más próximo (mediante una superficie de coste) y obtener la siguiente imágen mediante la interfaz GRASS-R.

#!/bin/sh

########### PRÁCTICA 4. SUPERFICIES DE COSTE ##################
### OBJETIVO: sobre el mapa vectorial de carreteras y calles del condado de Wake (streets_wakes) y sabiendo dónde están las estaciones de bomberos en dicho condado (gracias al mapa vectorial firestations), vamos a poder saber qué estación de bomberos llegará antes a un incendio en un punto determinado.

### METODOLOGÍA:
# •  Visualizar las estaciones de bomberos y las calles y carreteras sobre el condado de Wakes.
# •  Preparación de mapas para el cálculo: paso de vectorial a raster del mapa de calles y carreteras y del mapa de estaciones de bomberos.
# •  Conversión del mapa de calles y carreteras de millas por hora a kilómetros por hora.
# •  Cáculo del tiempo en cada celda del mapa de calles y carreteras .
# •  Calcular el mapa raster de la superficie de coste en tiempo.
# •  Cálculo del camino más rápido a un punto concreto.
# •  Conversión del camino mínimo a formato vectorial.
# •  Presentación de la información obtenida: resultado final.


## 1)• Condado de Wake: estaciones de bomberos y calles y carreteras
 d.mon x0
 g.region vect=firestations@PERMANENT
 d.rast boundary_county_500m@PERMANENT
 d.vect map=firestations@PERMANENT color=0:0:0 lcolor=0:0:0 fcolor=255:255:0 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1
 d.vect map=streets_wake@PERMANENT
 v.info -c map=streets_wake@PERMANENT layer=1 #deberemos ver cómo están conectadas y la rapidez de las vías que las unen. Mirar en el atributo SPEED del que disponemos.

## 2) • Preparación de los mapas de calles y carreteras y de estaciones de bomberos para el cálculo. Necesitamos convertir estos dos mapas vectoriales a mapas raster, pero atendiendo a los datos que nos interesan para cada uno de ellos.
 g.region vect=streets_wake res=30  #cambiamos la región y resolución de trabajo
 v.to.rast input=firestations@PERMANENT output=bomberos use=cat layer=1 value=1 rows=4096 # convertimos el mapa vectorial de estaciones a raster con una diferenciación entre las categorías de las estaciones de bomberos
 v.to.rast input=streets_wake@PERMANENT output=vias_mph use=attr column=SPEED layer=1 value=1 rows=4096 # convertimos el mapa vectorial de vías a raster con una diferenciación entre la velocidad (en mph) de las calles y carreteras

### • Conversión de unidades en el mapa vectorial de calles y carreteras. Una milla son 1'609 kilómetros
 r.mapcalculator amap=vias_mph formula=A*1.609 outfile=vias

### • Cálculo del tiempo por celda en el mapa ráster de calles y carreteras. Necesitamos calcular cuánto cuesta pasar (en tiempo) de una celda a otra para posteriormente calcular el coste global de ir de un punto cualquiera a otro del mapa.
 r.mapcalc "tceldaMinutos= if(isnull(vias), ewres()/50, ewres()/(1000*vias/60) )" # tiempo en min que cuesta atravesar una celda del mapa raster de vías (km/h): división entre la resolución del mapa y la velocidad 
 r.info vias
 ##Consideraciones:
 #• Para las celdas del mapa vias que no sean nulas, es decir, que pertenezcan a una calle o carretera, antes de dividir por las velocidades del mapa vias, las pasamos de km/h a m/min : (1000*vias/60) .
 #• Para las celdas nulas, las que no pertenezcan a una calle o carretera, vamos a dividir por una velocidad extremadamente lenta, para que atravesar esas casillas para alcanzar un punto del mapa no merezca la pena. Nos aseguramos, de todas formas, que esa velocidad esté fuera del rango de velocidades que representa nuestro mapa de vias (haciendo un r.info vias, como puedes ver en la siguiente imagen) y elegimos, por ejemplo, 3km/h, que son 50 m/min.

###• Calcular el mapa raster de la superficie de coste en tiempo
#como esta operación toma más de 6 horas trabajando, cambiamos la región activa al mapa raster elev_ned_30m#PERMANENT, que es una parte del condado de Wakes. En ella podremos calcular el camino propuesto por el ejercicio.
 g.region rast=elev_ned_30m res=30 -­a
 r.cost input=tceldaMinutos output=distanciaBomberos -­k start_rast=bomberos #obtenemos el mapa de distancia hasta las estaciones de bomberos

### • Cálculo del camino más rápido a un punto concreto.
# Vamos a suponer que el incendio se ha dado en la coordenada(637743,217989). Para calcular desde qué estación de bomberos se llegará antes (el camino mínimo), utilizamos la herramienta r.drain. 
 r.drain input=distanciaBomberos output=camino coordinate=637743,217989

### • Conversión del camino mínimo a formato vectorial.
 r.to.vect input=camino output=camino feature=line

### • Presentación de la información obtenida: resultado final.
 g.region vect=streets_wake@PERMANENT
 d.rast boundary_county_500m@PERMANENT
 d.vect map=streets_wake@PERMANENT 
 d.vect map=camino type=point,line,boundary,centroid,area,face display=shape icon=basic/x size=8 layer=1 width=0 wscale=1 color=white fcolor=200:200:200 rgb_column=GRASSRGB llayer=1 lcolor=red bgcolor=none bcolor=none lsize=8 font=romans xref=left yref=center #Dibujamos el camino encontrado sobre la red de calles y carreteras. 
 d.vect map=firestations@PERMANENT color=0:0:0 lcolor=0:0:0 fcolor=255:255:0 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=5 layer=1 lsize=8 xref=left yref=center llayer=1 #visualizamos el condado de Wakes y las estaciones de bomberos. 

 echo "637743 217989 1" | v.in.ascii out=fuego fs=space #Crearemos un mapa vectorial de puntos para indicar y visualizar el incendio. 
 d.vect map=fuego color=255:0:0 lcolor=0:0:0 fcolor=255:0:0 display=shape type=point,line,boundary,centroid,area icon=basic/diamond size=10 layer=1 lsize=8 xref=left yref=center llayer=1

 g.region rast=elev_ned_30m@PERMANENT
 d.his h_map=elev_ned_30m@PERMANENT i_map=elevation_shade@PERMANENT
 d.vect map=streets_wake@PERMANENT 
 d.vect map=camino col=white
 d.vect firestations col=blue icon=basic/circle size=10
 d.vect fuego col=red icon=basic/marker size=20

exi nviz elev_ned_30m vector=streets_wake,camino points=firestations,fuego #ahora en 3D