Análisis vectorial con sf
Daniela Ballari
Julio, 2018
Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ
Objetivo de la práctica: Explorar herramientas de análisis espaical vectorial a través de la identificación de 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. Las capas son provincias, centros de salud y areas inundables. Además se profundizará en la lectura de archivos de geoinformación, manejo de sistemas de referencia espaciales y visualización con la libreria sf.
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ó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(sf)## Warning: package 'sf' was built under R version 3.4.4setwd("C:/R_ecohidrologia/Vectorial") #definir directorio de trabajo A-Preparación de datos
A.1-Cargar capas.
Los datos a utilizar en esta práctica están disponibles Aqui
#PROVINCIAS
provincias<- st_read("nxprovincias.shp")#fuente editado de INEC archivo shp## Reading layer `nxprovincias' from data source `C:\R_ecohidrologia\Vectorial\nxprovincias.shp' using driver `ESRI Shapefile'
## Simple feature collection with 25 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -732143.5 ymin: 9445216 xmax: 1147852 ymax: 10189400
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs# summary(provincias)#explorar contenido de datos
#CENTROS DE SALUD
centros_salud<-st_read("AL_CENTOS_SALUD_2012.shp") #fuente: editado de INEC## Reading layer `AL_CENTOS_SALUD_2012' from data source `C:\R_ecohidrologia\Vectorial\AL_CENTOS_SALUD_2012.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2712 features and 28 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -80.986 ymin: -4.945079 xmax: -75.40171 ymax: 1.445067
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs# summary(centros_salud)
#AREAS INUNDABLES
inundacion<-st_read("AREA_INUNDACION.shp") #fuente: editado de INEC## Reading layer `AREA_INUNDACION' from data source `C:\R_ecohidrologia\Vectorial\AREA_INUNDACION.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1309 features and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 498741.5 ymin: 9503324 xmax: 1146668 ymax: 10161010
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs# summary(inundacion)Visualización de mapas básico
plot(provincias)
plot(centros_salud)## Warning: plotting the first 9 out of 28 attributes; use max.plot = 28 to
## plot all
plot(inundacion)
Explorando objetos sf
class(provincias)## [1] "sf" "data.frame"names(provincias)## [1] "DPA_PROVIN" "DPA_DESPRO" "DPA_VALOR" "DPA_ANIO" "REI_CODIGO"
## [6] "REN_CODIGO" "PEE_CODIGO" "preccount" "precsum" "precmean"
## [11] "geometry"provincias$DPA_DESPRO## [1] AZUAY BOLIVAR
## [3] CAÑAR CARCHI
## [5] COTOPAXI CHIMBORAZO
## [7] EL ORO ESMERALDAS
## [9] GUAYAS IMBABURA
## [11] LOJA LOS RIOS
## [13] MANABI MORONA SANTIAGO
## [15] NAPO PASTAZA
## [17] PICHINCHA TUNGURAHUA
## [19] ZAMORA CHINCHIPE GALAPAGOS
## [21] SUCUMBIOS ORELLANA
## [23] SANTO DOMINGO DE LOS TSACHILAS SANTA ELENA
## [25] ZONA NO DELIMITADA
## 25 Levels: AZUAY BOLIVAR CAÑAR CARCHI CHIMBORAZO COTOPAXI ... ZONA NO DELIMITADAprovincias$geometry## Geometry set for 25 features
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -732143.5 ymin: 9445216 xmax: 1147852 ymax: 10189400
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## First 5 geometries:## MULTIPOLYGON (((770262.3 9716934, 770300.8 9716...## MULTIPOLYGON (((737900.3 9868222, 737892.7 9868...## MULTIPOLYGON (((770262.3 9716934, 770204.1 9716...## MULTIPOLYGON (((778609 10132516, 778691.2 10132...## MULTIPOLYGON (((786958.2 9920679, 786958.7 9920...head(st_coordinates(provincias))## X Y L1 L2 L3
## [1,] 770262.3 9716934 1 1 1
## [2,] 770300.8 9716975 1 1 1
## [3,] 770358.2 9717035 1 1 1
## [4,] 770416.0 9717096 1 1 1
## [5,] 770422.6 9717101 1 1 1
## [6,] 770483.7 9717143 1 1 1A.2-Seleccionar área de estudio (Azuay)
azuay = provincias[provincias$DPA_DESPRO == "AZUAY", ]
provincias2 = st_union(provincias)
plot(st_geometry(provincias2), col="grey",axes = TRUE, reset=TRUE)
plot(st_geometry(st_centroid(provincias)), pch = 16, col = 'black', add = TRUE)## Warning in st_centroid.sf(provincias): st_centroid assumes attributes are
## constant over geometries of xplot(azuay["DPA_PROVIN"], col="red",add=TRUE)
Explorar sistemas de referencias
st_crs(provincias)## Coordinate Reference System:
## EPSG: 32717
## proj4string: "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs"st_crs(centros_salud)## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"st_crs(inundacion)## Coordinate Reference System:
## EPSG: 32717
## proj4string: "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs"# Si un SRS fue olvidado, se asigna con
nuevo_crs = st_set_crs(provincias, 32717) # set CRS
#Es necesario transformar el sistema de referencia de los centros de salud
centros_salud<- st_transform(centros_salud, 32717)plot(provincias2, col="lightgrey")
plot(centros_salud[1], pch=20,cex = 0.5, col="red",add=TRUE)
Calculo de superficie
st_area(azuay)## 8325667109 m^2attributes(st_area(azuay))## $units
## $numerator
## [1] "m" "m"
##
## $denominator
## character(0)
##
## attr(,"class")
## [1] "symbolic_units"
##
## $class
## [1] "units"st_area(azuay) / 1000000## 8325.667 m^2units::set_units(st_area(azuay), km^2)## 8325.667 km^2Transformar a data.frame
provincias.df = st_set_geometry(provincias, NULL)
class(provincias.df)## [1] "data.frame"head(provincias.df)## DPA_PROVIN DPA_DESPRO DPA_VALOR DPA_ANIO REI_CODIGO REN_CODIGO
## 1 01 AZUAY 0 2012 05 01
## 2 02 BOLIVAR 0 2012 02 01
## 3 03 CAÑAR 0 2012 05 01
## 4 04 CARCHI 0 2012 04 01
## 5 05 COTOPAXI 0 2012 02 01
## 6 06 CHIMBORAZO 0 2012 02 01
## PEE_CODIGO preccount precsum precmean
## 1 593 9 1307.4945 145.2772
## 2 593 5 1238.6102 247.7220
## 3 593 4 660.2622 165.0656
## 4 593 7 2050.4350 292.9193
## 5 593 9 2428.4772 269.8308
## 6 593 9 1529.6052 169.9561Extraer coordenadas
provincias$geometry## Geometry set for 25 features
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -732143.5 ymin: 9445216 xmax: 1147852 ymax: 10189400
## epsg (SRID): 32717
## proj4string: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## First 5 geometries:## MULTIPOLYGON (((770262.3 9716934, 770300.8 9716...## MULTIPOLYGON (((737900.3 9868222, 737892.7 9868...## MULTIPOLYGON (((770262.3 9716934, 770204.1 9716...## MULTIPOLYGON (((778609 10132516, 778691.2 10132...## MULTIPOLYGON (((786958.2 9920679, 786958.7 9920...geom.prov<-st_coordinates(provincias)
class(geom.prov)## [1] "matrix"head(geom.prov)## X Y L1 L2 L3
## [1,] 770262.3 9716934 1 1 1
## [2,] 770300.8 9716975 1 1 1
## [3,] 770358.2 9717035 1 1 1
## [4,] 770416.0 9717096 1 1 1
## [5,] 770422.6 9717101 1 1 1
## [6,] 770483.7 9717143 1 1 1Crear objeto sf
provincias.sf<- st_sf(data.frame(n = provincias$DPA_DESPRO), geom = provincias$geometry)
class(provincias.sf)## [1] "sf" "data.frame"provincias.sf2<- st_sf(data.frame(n = provincias$DPA_DESPRO), geom = st_sfc(st_multipoint(geom.prov[,c(1,2)])))
class(provincias.sf2)## [1] "sf" "data.frame"Seleccionar centros de salud localizados en el Azuay
centros_salud.azuay<- st_intersection(centros_salud, azuay)## Warning: attribute variables are assumed to be spatially constant
## throughout all geometriesplot(provincias2, col="lightgrey")
plot(centros_salud.azuay[1], pch=20,cex = 0.5, col="red",add=TRUE)
plot(st_geometry(azuay), col="lightgrey")
plot(st_geometry(centros_salud.azuay), pch=20,cex = 0.5, col="red",add=TRUE)
inundacion.azuay <- st_intersection(inundacion, azuay)## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries#visualizar resultados
plot(st_geometry(azuay), col="lightgrey",axes = TRUE)
plot(st_geometry(inundacion.azuay), pch=20,cex = 0.5, col="blue",border="blue",add=TRUE)
plot(st_geometry(centros_salud.azuay), pch=20,cex = 0.5, col="red",add=TRUE)
Area de zonas inundables
st_area(inundacion.azuay)## Units: m^2
## [1] 728576.2 1954415.8 5051608.6 7088599.7 1107148.6 15423539.8
## [7] 7642272.1 18232166.0 17799634.5units::set_units(st_area(inundacion.azuay), km^2)## Units: km^2
## [1] 0.7285762 1.9544158 5.0516086 7.0885997 1.1071486 15.4235398
## [7] 7.6422721 18.2321660 17.7996345B-Análisis espacial
B.1-Extender área de centros de salud buffer 100m
buffer<- st_buffer(centros_salud.azuay, 100)#cálculo de buffers de 100m
class(buffer)#tipo de dato devuelto## [1] "sf" "data.frame"B.2-Selección de buffers que se superponen con áreas inundables
buffer_inundacion <- st_intersection(buffer,inundacion.azuay) ## Warning: attribute variables are assumed to be spatially constant
## throughout all geometriescentros_salud_riesgo_centroides <- st_centroid(buffer_inundacion) ## Warning in st_centroid.sf(buffer_inundacion): st_centroid assumes
## attributes are constant over geometries of x#visualizar resultados
plot(st_geometry(azuay), col="lightgrey",axes = TRUE)
plot(st_geometry(inundacion.azuay), pch=20,cex = 0.5, col="blue",border="blue",add=TRUE)
plot(st_geometry(centros_salud_riesgo_centroides), pch=20,cex = 1, col="yellow",add=TRUE)
C-Otros procesos complementarios
C.1-Crear nuevos campos y atributos
buffer_inundacion$riesgo <- "inundacion" #cear un campo 'riesgo' cuyo valor es constante para todos los elementos
buffer_inundacion$id <- 1:nrow(buffer_inundacion) #cear un campo 'id' cuyo valor se asigna de forma secuencial desde 1 hasta el numero todal de elementos contenidos en el'SpatialPolygonsDataFrame'.
inundacion.azuay$area_m <- st_area(inundacion.azuay) # En m2
inundacion.azuay$area_m <- st_area(inundacion.azuay) # En m2
inundacion.azuay$area_km <- units::set_units(st_area(inundacion.azuay), km^2) # en km2
sum(inundacion.azuay$area_km)#area total inundable en el azuay 75km2## 75.02796 km^2head(centros_salud.azuay$uni_nombre)#centros de salud en riesgo## [1] CHILCAPAMBA ZHUCAY GULLANZHAPA
## [4] REMIGIO CRESPO SIMON BOLIVAR CENTRO DE SALUD NO 1
## 2495 Levels: 1-E 10 DE AGOSTO 10 DE NOVIEMBRE ... ZURMITrabajo Calificado 6p
Localiza datos de tipo vectorial (puntos, lineas y polígonos) para tu zona de estudio. Por ejemplo estaciones metereológicas o pluviómetros (puntos), lagos (polígonos), usos de suelo (polígonos), rios (líneas) o caminos (líneas). Cargue los datos en R, explore los sistemas de referencia, recorte los datos para su zona de estudio y calcule áreas (para polígonos) o longitud (para líneas).
Has llegado al final de este material!