jueves, 25 de noviembre de 2010

Análisis estadístico y localización. 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 el análisis estadístico y localización de un parque de bomberos al condado de Wakes (estado de Carolina del Norte).
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

###### PRÁCTICA 5. ANÁLISIS ESTADÍSTICO Y LOCALIZACIÓN #######

#OBJETIVO: añadiremos un parque de bomberos nuevo al condado de Wakes (estado de Carolina del Norte). Emplearemos una serie de criterios para emplazarlo: que maximice la distancia al resto de parques y que esté en un sitio llano y cerca de lagos. Usaremos también herramientas estadísticas básicas.

#PASOS:
# 1) Calcular las zonas cercanas a lagos del condado de Wakes.
# 2) Determinación de zonas llanas en el condado de Wakes.
# 3) Aunar ambos criterios: zonas llanas cerca de lagos.
# 4) Cáculo del mapa de distancias para las zonas seleccionadas .
# 5) Selección del emplazamiento de la nueva estación de bomberos.
# 6) Presentación de la información obtenida: resultado final.

# 1) Calcular las zonas cercanas a lagos del condado de Wakes.
d.mon x1
d.rast map=boundary_county_500m@PERMANENT -o
d.vect map=streets_wake@PERMANENT fcolor=170:170:170 #color=0:0:0 lcolor=0:0: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=firestations@PERMANENT fcolor=255:255:0 size=10 #color=0:0:0 lcolor=0:0:0 display=shape type=point,line,boundary,centroid,area icon=basic/circle layer=1 lsize=8 xref=left yref=center llayer=1
d.vect map=lakes@PERMANENT fcolor=0:0:255 size=3 # color=0:0:0 lcolor=0:0:0 display=shape type=point,line,boundary,centroid,area icon=basic/circle size=3 layer=1 lsize=8 xref=left yref=center llayer=1
g.region vect=streets_wake res=30 ­-a #AJUSTAMOS LA REGIÓN
v.to.rast input=lakes@PERMANENT output=wakes_lakes use=val layer=1 value=1 rows=4096 #convertimos el mapa de lagos vectorial a un mapa de lagos raster para todo el condado de Wakes ocupado por calles y carreteras

r.colors map=wakes_lakes color=rules < "/home/usuario/CursoGIS/actividad5/rule" #cambiamos el color a azul
r.buffer input=wakes_lakes output=cerca_lagos distances=250 units=meters --­­overwrite #delimitar una zona alejada de los lagos. Vamos a considerar, por ejemplo, que 250 metros es una distancia cercana.

# 2) Determinación de zonas llanas en el condado de Wakes.
#Asumimos que una estación de bomberos ocupa unos 8000 m2, que deseamos además en terreno llano (pendiente nula o casi nula). Para un área cuadrada de este tamaño, calculamos el valor del lado: Si A = l2 = 8000 m2 → l = √8000 = 89,4 m y nos quedaremos con el terreno llano (donde hay muy poca diferencia de altura entre sus puntos -0.1 metros, por ejemplo-)
#Usamos r.neighbors para encontrar las zonas que cumplen ese criterio. Especificamos el tamaño de la ventana de cálculo (size) de 3x3 celdas para la ventana (ventana de 90 metros de lado -=3x30m-).
g.region rast=elev_ned_30m res=30 ­-a #restringimos la zona de estudio a la zona del mapa de elevaciones

#calculamos las máximas y mínimas elevaciones para cada celda en el mapa de elevaciones
r.neighbors input=elev_ned_30m@PERMANENT output=sWake_elev_max method=maximum size=3 --overwrite
r.neighbors input=elev_ned_30m@PERMANENT output=sWake_elev_min method=minimum size=3 --overwrite
r.mapcalc "llanas=if(sWake_elev_max­-sWake_elev_min<0.1,1,null())"

# 3) Aunar ambos criterios: zonas llanas cerca de lagos.
r.mapcalc "llanas_y_lagos=if(llanas*cerca_lagos==2,1,null())" #Calculamos las zonas llanas y cerca del lago: la categoría 2 del mapa "cerca_lagos", alberga todos los puntos del mapa que están a menos de 250 metros de un lago y la categoría 1 del mapa "llanas" representa las zonas llanas.

# 4) Cáculo del mapa de distancias para las zonas seleccionadas.
r.mapcalc "coste1=ewres()" #coste unitario, en distancia, es la cantidad de espacio que nos ocupa una celda, el espacio que tenemos que recorrer al atravesarla, que se corresponde con la resolución a la que estamos trabajando.
r.cost -­k in=coste1 out=toFirestations start_rast=bomberos@actividad4 #calculamos la superficie de coste en distancia, utilizando el mapa raster de las estaciones de bomberos de la actividad anterior (no 4), para toda la región del condado que estamos contemplando.

#Queremos saber el coste en metros de atravesar desde zonas llanas y cerca de rios a las estaciones de bomberos, en ese mapa de costes recién calculado. Esas zonas las tenemos ya calculadas en el mapa llanas_y_lagos. Si multiplicamos este mapa por el de costes calculados (toFirestations), los valores nulos del mapa de llanas_y_lagos eliminarán las zonas que no nos interesan del mapa de costes toFirestations, y la categoria 1 -celdas que valen 1- nos respetará los valores del mapa de costes que nos interesan:
r.mapcalc "mtoFirestations=llanas_y_lagos*toFirestations"
d.rast mtoFirestations
d.legend map=mtoFirestations color=black lines=0 thin=1 labelnum=5 at=80,95,2,10

# 5) Selección del emplazamiento de la nueva estación de bomberos.
r.univar map=mtoFirestations # cual o cuáles de los puntos calculados en el mapa mtoFirestations tienen el valor más alto: que vale 4170.
r.stats ­-1lngx input=mtoFirestations fs=space nv=* nsteps=255 #localizamos las coordenadas a las que se refiere la categoria 4170 (metros). Para la categoría 4170, nos encontramos dos puntos. Como son contiguos, nos quedamos con la coordenada espacial de uno de ellos.
echo "637785 218475 1" | v.in.ascii out=nuevaFS fs=space #convertimos en un mapa vectorial de un único punto (como hicimos en la actividad anterior).

# 6) Presentación de la información obtenida: resultado final.
g.region zoom=elev_ned_30m
d.his h_map=elev_ned_30m@PERMANENT i_map=elevation_shade@PERMANENT
d.vect streets_wake
d.vect map=firestations@PERMANENT color=255:255:0 lcolor=0:0:0 fcolor=255:255:0 display=shape type=point,line,boundary,centroid,area icon=basic/diamond size=8 layer=1 lsize=8 xref=left yref=center llayer=1
d.vect map=nuevaFS color=0:0:255 lcolor=0:0:0 fcolor=0:0:153 display=shape type=point,line,boundary,centroid,area
icon=basic/star size=15 layer=1 lsize=8 xref=left yref=center llayer=1

nviz elev_ned_30m vector=streets_wake points=firestations,nuevaFS #en 3D