Como importar Datos en Formato Shapefile en R?

Para importar archivos en formato shapefile en R debemos primero instalar las librerias raster y rgdal. Luego podemos usar la funciĆ³n shapefile y como argumento la ruta del archivo shapefile que vamos a cargar, en el ejemplo datos de los poligonos de las comunas de la ciudad de Cali y una visualizaciĆ³n de este mapa.

##Importar Datos Espaciales en R

require(raster)
require(rgdal)

comunas=shapefile("~/Desktop/shapes/Comunas.shp")
comunas
## class       : SpatialPolygonsDataFrame 
## features    : 22 
## extent      : 1053868, 1068492, 860190.2, 879441.5  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=3.441883333 +lon_0=-76.5205625 +k=1 +x_0=1061900.18 +y_0=872364.63 +a=6379137 +rf=298.257222101 +units=m +no_defs +type=crs 
## variables   : 4
## names       : OBJECTID, gid, comuna,   nombre 
## min values  :        1,  89,      1, Comuna 1 
## max values  :       22, 110,     22, Comuna 9
plot(comunas)

Dentro de los atributos del archivo contamos con la informaciĆ³n de 22 poligonos y 4 variables. Para realizar una consulta espacial y filtrar solo uno o varios poligonos de interes usamos la siguiente instruciĆ³n, similar a un data.frame:

comuna22=comunas[22,]
plot(comuna22)

Tambien podemos usar un campo de la tabla de atributos para realizar la consulta espacial, por ejemplo si queremos filtrar solo las comunas 20, 21 y 22 podemos usar lo siguiente:

comuna202122=comunas[comunas$comuna>=20,]
plot(comunas)
plot(comuna202122,col="red",add=TRUE)

Si queremos exportar este nuevo archivo shapefile podemos usar la misma funciĆ³n para importarlo de esta manera:

shapefile(comuna202122,"~/Desktop/shapes/Comunas_nuevo.shp",overwrite=TRUE)

Como crear un shapefile apartir de un archivo de coordenadas de puntos?

Lo primero es cargar la base de datos de ejemplo que corresponde a datos de ofertas de vivienda en la ciudad de Cali de la siguiente manera:

require(readxl)
Datos_Vivienda = read_excel("~/Desktop/shapes/Datos_Vivienda.xlsx")

Ahora vamos a crear el archivo shapefile con geometria de puntos con base en estos datos (tener en cuenta el archivo no puede tener valores faltantes)

Datos_Vivienda2=na.omit(Datos_Vivienda)
vivienda_map=SpatialPointsDataFrame(coords = Datos_Vivienda2[,11:12],
                       data = Datos_Vivienda2,
                       proj4string =crs(comunas))

vivienda_map
## class       : SpatialPointsDataFrame 
## features    : 8318 
## extent      : -76.58915, -76.463, 3.333, 3.4977  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=3.441883333 +lon_0=-76.5205625 +k=1 +x_0=1061900.18 +y_0=872364.63 +a=6379137 +rf=298.257222101 +units=m +no_defs +type=crs 
## variables   : 12
## names       :        Zona, piso, Estrato, precio_millon, Area_contruida, parqueaderos, Banos, Habitaciones,        Tipo,      Barrio, cordenada_longitud, Cordenada_latitud 
## min values  : Zona Centro,    1,       3,            58,             30,            1,     0,            0, Apartamento, 20 de julio,          -76.58915,             3.333 
## max values  :    Zona Sur,   NA,       6,          1999,           1745,           NA,    10,           10,        Casa,    zona sur,            -76.463,            3.4977
plot(vivienda_map)

De manera similar podemos hacer consultas sobre este nuevo archivo shapefile y filtrar solo casas de estrato 6:

pos_casas=which(vivienda_map@data$Tipo=="Casa"&
                  vivienda_map@data$Estrato==6)

vivienda_casas=vivienda_map[pos_casas,]
plot(vivienda_casas)

Podemos visalizarlo de manera interactiva utilizando la libreria leaflet y recordar que esta solo funciona con shapes que esten con sistema de referencia de coordenadas geografico:

require(leaflet)
## Loading required package: leaflet
leaflet() %>% addTiles() %>% 
  addCircleMarkers(lng = vivienda_casas$cordenada_longitud,
                   lat =vivienda_casas$Cordenada_latitud ,
                   radius = 0.1)

Aqui les recomiendo jugar un poco con otros argumentos que tiene la libreria leaflet para visualizar por ejemplo usando agregaciĆ³n con poligonos interactivos:

require(leaflet)
leaflet() %>% addTiles() %>% 
  addCircleMarkers(lng = vivienda_casas$cordenada_longitud,
                   lat =vivienda_casas$Cordenada_latitud ,
                   radius = 0.1,clusterOptions = markerClusterOptions() )

Como poder integrar datos a nivel espacial?

En este caso vamos a utilizar la funciĆ³n over para integrar los dos shapes y realizar calculos para contar por ejemplo el total de viviendas por comuna y visualizarlas en un mapa de coropletas:

Antes de iniciar es importante tener en cuenta que los sistemas de referencia deben ser iguales y para esto modelos usar transformaciĆ³n de coordenadas con la funciĆ³n spTransform:

comunas2=spTransform(comunas,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

vivienda_map2=vivienda_map
crs(vivienda_map2)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

total_ofertas=over(x =comunas2 ,y = vivienda_map2[,1],fn =length)
comunas2$ofertas=total_ofertas[,1]
head(comunas2@data)
OBJECTID gid comuna nombre ofertas
0 1 107 2 Comuna 2 1445
1 2 108 1 Comuna 1 132
2 3 109 3 Comuna 3 449
3 4 110 19 Comuna 19 1092
4 5 103 15 Comuna 15 106
5 6 104 17 Comuna 17 2164

Como pueden ver en la tabla de atributos ahora contamos con la variable total de ofertas por comuna la cual podemos visualizar con un mapa de coropletas:

spplot(comunas2[,5])

Se hace referencia a la columna 5 por que es la que se crea nueva con la variable total de ofertas.