Objetivo de la practica. Determinar los centros de salud de la provincia del Azuay que se encuentran localizados en areas inundables. Se utilizarán tres capas de información. Ellas son provincias, centros de salud y areas inundables.
Pasos a seguir
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”
library(rgdal) #http://cran.r-project.org/web/packages/rgdal/rgdal.pdf## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/projlibrary(sp)
library(rgeos) #http://cran.r-project.org/web/packages/rgeos/rgeos.pdf## Warning: package 'rgeos' was built under R version 3.0.3## rgeos version: 0.3-6, (SVN revision 450)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Polygon checking: TRUElibrary(maptools)## Warning: package 'maptools' was built under R version 3.0.3## Checking rgeos availability: TRUEsetwd("C:/curso_R/") #definir directorio de trabajo
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFbmNIWnZoaFk0c3M")#descargar datos
#ogrListLayers(".")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 and 10 fields
## Feature type: wkbPolygon with 2 dimensionssummary(provincias)#explorar contenido de datos## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -732143.5 1147852
## y 9445216.4 10189398
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 01 : 1 AZUAY : 1 Min. :0 2012:25 02:5 01:11
## 02 : 1 BOLIVAR : 1 1st Qu.:0 03:6 02: 6
## 03 : 1 CAÃAR : 1 Median :0 04:8 03: 6
## 04 : 1 CARCHI : 1 Mean :0 05:5 04: 1
## 05 : 1 CHIMBORAZO: 1 3rd Qu.:0 90:1 05: 1
## 06 : 1 COTOPAXI : 1 Max. :0
## (Other):19 (Other) :19
## PEE_CODIGO preccount precsum precmean
## 593:25 Min. : 2.00 Min. : 307.6 Min. : 76.91
## 1st Qu.: 5.75 1st Qu.: 1224.6 1st Qu.:175.78
## Median : 9.00 Median : 2339.4 Median :228.27
## Mean :13.67 Mean : 3355.4 Mean :232.67
## 3rd Qu.:20.25 3rd Qu.: 4571.0 3rd Qu.:283.81
## Max. :40.00 Max. :14596.2 Max. :434.52
## NA's :1 NA's :1 NA's :1proj4string(provincias)#explorar sistema de referencia## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"utm17 <- "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs" #Sistema de referencia a utilizar utm zona 17 sur con datum WGS84
provincias <- spTransform(provincias,CRS(utm17))# transformar CRS para que todas las capas tengan la misma definición
#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 and 28 fields
## Feature type: wkbPoint with 2 dimensionssummary(centros_salud)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -80.985997 -75.401706
## coords.x2 -4.945079 1.445067
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 2712
## Data attributes:
## uni_codigo uni_nombre uni_direcc
## 000001 : 1 SAN PABLO : 8 CENTRO DE LA COMUNIDAD: 47
## 000002 : 1 BELLAVISTA : 7 S/D : 44
## 000003 : 1 SAN VICENTE: 7 S/N CENTRO DEL PUEBLO : 42
## 000004 : 1 SANTA ROSA : 7 CALLE PRINCIPAL : 29
## 000005 : 1 ATAHUALPA : 6 N/A : 25
## 000006 : 1 LA VICTORIA: 5 (Other) :2522
## (Other):2706 (Other) :2672 NA's : 3
## are_codigo are_descri
## 06A01 : 56 RIOBAMBA CHAMBO: 56
## 13A01 : 49 PORTOVIEJO : 49
## 03A01 : 45 AZOGUES : 45
## 13A03 : 43 CHONE : 43
## 10A01 : 41 IBARRA : 41
## 13A02 : 40 MANTA : 40
## (Other):2438 (Other) :2438
## par_descri par_urbano can_descri
## QUITO : 114 R:1496 QUITO : 173
## GUAYAQUIL : 109 U:1216 GUAYAQUIL : 121
## PORTOVIEJO : 26 CUENCA : 59
## SANTO DOMINGO DE LOS COLORADOS: 23 RIOBAMBA : 54
## CUENCA : 22 PORTOVIEJO : 49
## ESMERALDAS : 22 SANTO DOMINGO: 46
## (Other) :2396 (Other) :2210
## prv_descri dis_distri cir_codigo
## MANABI : 305 RIOBAMBA-CHAMBO : 57 16D01C03: 16
## GUAYAS : 266 PASTAZA-MERA-SANTA CLARA: 53 04D03C03: 12
## PICHINCHA : 233 PORTOVIEJO : 49 06D04C03: 11
## LOJA : 186 AZOGUES-BIBLIAN-DELEG : 44 03D02C02: 10
## CHIMBORAZO: 162 CHONE-FLAVIO ALFARO : 43 06D05C01: 10
## ESMERALDAS: 154 COLTA-GUAMOTE : 43 02D02C02: 9
## (Other) :1406 (Other) :2423 (Other) :2644
## hau_descri nun_descri tun_descri
## 12 HORAS : 65 PRIMER NIVEL :2522 CENTRO DE SALUD :1478
## 24 HORAS : 227 SEGUNDO NIVEL: 186 HOSPITAL BASICO : 114
## 8 HORAS :1824 TERCER NIVEL : 4 HOSPITAL DE ESPECIALIDADES: 2
## No definido: 596 HOSPITAL ESPECIALIZADO : 16
## HOSPITAL GENERAL : 58
## PUESTO DE SALUD :1044
##
## sun_descri igu_descri FID_1
## No definido: 22 MSP :1993 Min. : 0.0
## PUBLICO :2690 IESS : 686 1st Qu.: 250.0
## ISSFA : 13 Median : 570.0
## SOLCA : 7 Mean : 564.6
## JUNTA DE BENEFICENCIA: 4 3rd Qu.: 845.0
## ONG : 3 Max. :1133.0
## (Other) : 6
## ZONA COD_CIRC_1 DPA_PROVIN DPA_DESPRO
## Min. :0.000 16D01C03: 16 13 : 305 MANABI : 305
## 1st Qu.:3.000 04D03C03: 12 09 : 266 GUAYAS : 266
## Median :5.000 06D04C03: 11 17 : 233 PICHINCHA : 233
## Mean :4.613 03D02C02: 10 11 : 186 LOJA : 186
## 3rd Qu.:7.000 06D05C01: 10 06 : 162 CHIMBORAZO: 162
## Max. :9.000 (Other) :2647 (Other):1554 ESMERALDAS: 154
## NA's : 6 NA's : 6 (Other) :1406
## DPA_CANTON DPA_DESCAN DPA_PARROQ
## 1701 : 173 QUITO : 173 170150 : 114
## 0901 : 121 GUAYAQUIL : 121 090150 : 97
## 0101 : 59 CUENCA : 59 130150 : 26
## 0601 : 54 RIOBAMBA : 54 230150 : 23
## 1301 : 49 PORTOVIEJO : 49 010150 : 22
## 2301 : 46 SANTO DOMINGO: 46 (Other):2425
## (Other):2210 (Other) :2210 NA's : 5
## DES_CIRC POB_CIRC X
## PORTOVIEJO : 26 Min. : 0 Min. :-80.99
## SANTO DOMINGO DE LOS COLORADOS: 23 1st Qu.: 5715 1st Qu.:-79.79
## ESMERALDAS : 22 Median :10245 Median :-79.06
## LOJA : 22 Mean :13410 Mean :-79.08
## MANTA : 21 3rd Qu.:15832 3rd Qu.:-78.54
## MACHALA : 19 Max. :92260 Max. :-75.40
## (Other) :2579
## Y
## Min. :-4.9451
## 1st Qu.:-2.5067
## Median :-1.5084
## Mean :-1.5437
## 3rd Qu.:-0.3563
## Max. : 1.4451
## centros_salud <- spTransform(centros_salud,CRS(utm17))#transformar CRS
#AREAS INUNDABLES
inundacion<-readOGR(".","AREA_INUNDACION") #fuente: editado de INEC## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "AREA_INUNDACION"
## with 1309 features and 2 fields
## Feature type: wkbPolygon with 2 dimensionssummary(inundacion)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 498741.5 1146668
## y 9503323.5 10161011
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## DESCRIPCIO
## ZONAS INUNDADAS PERMANENTEMENTE (MANGLARES Y PANTANOS) :310
## ZONAS INUNDADAS TEMPORALMENTE (CADA EPOCA LLUVIOSA) : 12
## ZONAS PROPENSAS A INUNDACIONES (DESBORDAMIENTO DE RIOS O FUERTES PRESCIPITACIONES):987
##
##
##
## HECTARES
## Min. : 1.0
## 1st Qu.: 196.3
## Median : 418.1
## Mean : 3186.8
## 3rd Qu.: 1278.2
## Max. :430406.9inundacion <- spTransform(inundacion,CRS(utm17))#transformar CRS A.2-Seleccionar área de estudio (Azuay), Intersectar para cortar capas de información
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, sin embargo para capas de poligonos o lineas puede no ser la funcion desada. Esto es porque selecciona el poligono o linea completa, sin cortar con la superfice de interseccion. A continuacion se muestra el resultado utilizando gIntersects (sin cortar poligonos) y posteriormente con gIntersection (cortando poligonos).
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 funcion gIntersects, selecciona los poligonos que intersectan al Azuay, sin embargo no los corta. Por ello, para obtener poligonos cortados con el elemento de interseccion 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="red", 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)
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
#visualizar resultados
polygons1 = list("sp.polygons", buffer_inundaciondf, col="red", fill="red")
polygons2 = list("sp.polygons", inundacion_intersect_azuay, col="blue", fill="blue")
spplot(azuay, "DPA_DESPRO", col="white", fill="grey", sp.layout = list(polygons2,polygons1),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.02796centros_salud_riesgo_centroides <- gCentroid(inundacion_intersect_azuay, byid=TRUE) # centroides de los buffers de centros de salud en riesgos. http://www.inside-r.org/packages/cran/rgeos/docs/gCentroid
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"