Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/

Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications

Curso completo en: http://rpubs.com/daniballari/analisis_espacial_indice



Objetivo de la práctica: Identificar los centros de salud de la provincia del Azuay que se encuentran localizados en areas inundables. Se utilizarán tres capas de información y operaciones vectoriales como son buffers, intersecciones, etc.. Ellas son provincias, centros de salud y areas inundables.

Temario

A-Preparación de datos
A.1-Cargar capas y editar SRS
A.2-Seleccionar área de estudio (Azuay).Intersección para cortar capas de información
B-Análisis espacial
B.1-Extender área de centros de salud buffer 100m
B.2-Consulta espacial de buffers que se superponen con áreas inundables
C-Otros procesos complementarios
C.1-Crear nuevos campos y atributos
C.2-Combinar capas resultado

¿Cómo hacer esta práctica con Qgis? Aqui

¿Teoría sobre analisis espacial vectorial?
Olaya Victor (Ed.) (2011) Sistemas de Información Geográfica, Capítulo 18 “Operaciones geométricas con datos vectoriales”

Cargar librerias y datos

library(rgdal) 
library(sp)
library(rgeos) 
library(maptools)

setwd("C:/R_ecohidrologia/Vectorial") #definir directorio de trabajo

A-Preparación de datos

A.1-Cargar capas

#PROVINCIAS
provincias<-readOGR(".","nxprovincias") #fuente editado de INEC archivo shp
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "nxprovincias"
## with 25 features
## It has 10 fields
# summary(provincias)#explorar contenido de datos

#CENTROS DE SALUD
centros_salud<-readOGR(".","AL_CENTOS_SALUD_2012") #fuente: editado de INEC, archivo shp
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "AL_CENTOS_SALUD_2012"
## with 2712 features
## It has 28 fields
# summary(centros_salud)

#AREAS INUNDABLES
inundacion<-readOGR(".","AREA_INUNDACION") #fuente: editado de INEC
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "AREA_INUNDACION"
## with 1309 features
## It has 2 fields
# summary(inundacion)
#Explorar sistemas de referencias
proj4string(provincias)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(centros_salud)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(inundacion)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
#Es necesario transformar el sistema de referencia de los centros de salud

utm17 <-"+init=epsg:32717" #Sistema de referencia a utilizar utm zona 17 sur con datum WGS84
centros_salud <- spTransform(centros_salud,CRS(utm17))#transformar CRS

# Transformaremos los sistemas de referenicas de las demás capas para que tengan la misma definición
provincias <- spTransform(provincias,CRS(utm17))
inundacion <- spTransform(inundacion,CRS(utm17))

A.2-Seleccionar área de estudio (Azuay), Intersectar para cortar capas de información

plot(provincias)

names(provincias)#explorar nombres de campos de variable provincias
##  [1] "DPA_PROVIN" "DPA_DESPRO" "DPA_VALOR"  "DPA_ANIO"   "REI_CODIGO"
##  [6] "REN_CODIGO" "PEE_CODIGO" "preccount"  "precsum"    "precmean"
azuay <- subset(provincias, DPA_DESPRO=='AZUAY')# Se selecciona solo aquellos registros que tengan el valor de 'AZUAY' en el campo 'DPA_DESPRO'  
plot(azuay)  

centros_salud$azuay <- as.vector(gIntersects(centros_salud, azuay, byid=TRUE))#intersectar los centros de salud con el polígono de Azuay. Devuelve valores TRUE para los puntos que están contenidos en Azuay. #La funcion gIntersects hace que los datos resultantes pierdan sus atributos, por ello se realiza el siguiente artilugio para recuperar sus atributos: El resultado TRUE o FALSE es asignado a un nuevo campo llamado azuay (centros_salud$azuay) y posteriormente se crea una nueva capa con la seleccion (subset) de los valores TRUE del campo azuay.     


centros_salud_azuay <- subset(centros_salud, azuay==TRUE) 

#La funcion gIntersects es util para interseccion de capas de tipo puntual(como lo es la capa de centros de salud), sin embargo para capas de poligonos o lineas (como lo es la capa de inundación)  no es la función deseada. Esto es porque selecciona el polígono o línea completa, sin cortar con la superfice de interseccion. A continuación se muestra el resultado utilizando gIntersects (sin cortar polígonos) y posteriormente con gIntersection (cortando polígonos).  

inundacion$azuay <- as.vector(gIntersects(inundacion, azuay, byid=TRUE))
inundacion_azuay <- subset(inundacion, azuay==TRUE)

#visualizar resultados
points = list("sp.points", centros_salud_azuay, col="red", pch=19)
polygons = list("sp.polygons", inundacion_azuay, col="blue", fill="blue")
spplot(azuay, "DPA_DESPRO", col="white", fill="grey", sp.layout = list(polygons,points),colorkey=FALSE)

# En la figura se observa que la función gIntersects, selecciona los polígonos que intersectan al Azuay, sin embargo no los corta. Por ello, para obtener polígonos cortados con el elemento de intersección se realizan los siguientes pasos:  
inundacion_intersect <- gIntersection(inundacion, azuay, byid=c(TRUE,FALSE))# byid=c(TRUE,FALSE) permite almacenar el ID de los elementos a conservar. 

#En este caso para recuperar los atributos originales de los elementos intersectados se debe recurrir a la funcion 'spCbind' de la libreria 'maptools'  

#Libreria que contiene la funcion spCbind para unir (acoplar = bind) columnas a un SpatialPolygonsDataFrame ya existente. http://cran.r-project.org/web/packages/maptools/maptools.pdf

inundacion_intersect_azuay <- spCbind(as(inundacion_intersect, "SpatialPolygonsDataFrame"),inundacion@data[row.names(inundacion_azuay),] ) # se crea una nueva capa 'inundacion_intersect_azuay' que es la union de 'inundacion_intersect'(previamente transformado a SpatialPolygonsDataFrame) con la tabla de atributos (data.frame) de inundacion. Para evitar unir data.frames con diferentes número de elementos se seleccionan los elementos que pertenezcan al azuay (capa inundacion_azuay), mediante al función row.names. Dado que row.names (id de cada fila) es un elemento común en todas las capas de inundaciones, la unión se realiza a través de este id.  

#visualizar resultados
points = list("sp.points", centros_salud_azuay, col="green", pch=19)
polygons = list("sp.polygons", inundacion_intersect_azuay, col="blue", fill="blue")
spplot(azuay, "DPA_DESPRO", col="white", fill="grey", sp.layout = list(polygons,points),colorkey=FALSE)

l1 = list("sp.points", centros_salud_azuay, pch = 19, col="green")
polygons = list("sp.polygons", azuay, col="white", fill="grey")
spplot(inundacion_intersect_azuay, "azuay", col="blue", fill="blue", sp.layout = list(l1, polygons),colorkey=FALSE)

B-Análisis espacial

B.1-Extender área de centros de salud buffer 100m

buffer <- gBuffer(centros_salud_azuay, byid=TRUE, width=100) #cálculo de buffers de 100m
class(buffer)#tipo de dato devuelto
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

B.2-Selección de buffers que se superponen con áreas inundables

buffer_inundacion <-  gIntersection(buffer,inundacion_intersect_azuay,byid=c(TRUE, FALSE))
buffer_inundaciondf<- spCbind (as(buffer_inundacion, "SpatialPolygonsDataFrame"),buffer@data[row.names(buffer_inundacion),] )#recuparar atributos

centros_salud_riesgo_centroides <- gCentroid(buffer_inundaciondf, byid=TRUE) # centroides de los buffers de centros de salud en riesgos. http://www.inside-r.org/packages/cran/rgeos/docs/gCentroid  

#visualizar resultados
l1 = list("sp.points", centros_salud_azuay, pch = 19, col="green")
l2= list("sp.points", centros_salud_riesgo_centroides, pch = 19, col="red")
polygons = list("sp.polygons", azuay, col="white", fill="grey")
spplot(inundacion_intersect_azuay, "azuay", col="blue", fill="blue", sp.layout = list(l1, l2, polygons),colorkey=FALSE)

C-Otros procesos complementarios

C.1-Crear nuevos campos y atributos

buffer_inundaciondf$riesgo <- "inundacion" #cear un campo 'riesgo' cuyo valor es constante para todos los elementos
buffer_inundaciondf$id <- 1:length(buffer_inundaciondf) #cear un campo 'id' cuyo valor se asigna de forma secuencial desde 1 hasta el numero todal de elementos contenidos en el'SpatialPolygonsDataFrame'.  
inundacion_intersect_azuay$area_m <- gArea(inundacion_intersect_azuay, byid=TRUE) # En unidades del sistema de referencia, en este caso metros. http://www.inside-r.org/packages/cran/rgeos/docs/gArea

inundacion_intersect_azuay$area_km <- 1.0e-6 * gArea(inundacion_intersect_azuay, byid=TRUE) # en km2
sum(inundacion_intersect_azuay$area_km)#area total inundable en el azuay 75km2
## [1] 75.02796
centros_salud_riesgo<- centros_salud[row.names(centros_salud_riesgo_centroides),] # subset la variable centros de salud original para obtener coordenadas y atributos de los centros de salud en riesgo.
names(centros_salud_riesgo)
##  [1] "uni_codigo" "uni_nombre" "uni_direcc" "are_codigo" "are_descri"
##  [6] "par_descri" "par_urbano" "can_descri" "prv_descri" "dis_distri"
## [11] "cir_codigo" "hau_descri" "nun_descri" "tun_descri" "sun_descri"
## [16] "igu_descri" "FID_1"      "ZONA"       "COD_CIRC_1" "DPA_PROVIN"
## [21] "DPA_DESPRO" "DPA_CANTON" "DPA_DESCAN" "DPA_PARROQ" "DES_CIRC"  
## [26] "POB_CIRC"   "X"          "Y"          "azuay"

Has llegado al final de este material!