De datos Geo-Espaciales, Mapas y demás

Vamos a ver un proyecto muy interesante que hicimos el máster de Ciencia de Datos de la Universidad de Valencia, en la asignatura de Visualización Avanzado de Datos.

La idea del proyecto es aprender a manipular datos geo-espaciales, archivos formato shape, paquetes de mapas y obtener información geo-espacial.

Se utilizaron diferentes librerías de mapas (tmap, ggmap y leaflet). La mayoría de los mapas se hicieron con leaflet por el atractivo de los mapas que genera así como la posibilidad de hacerlos interactivos.

Las librerías que se utilizaron son:

library(ggmap)
library(sp)
library(rgdal)
library(rgeos)
library(maptools)
library(dplyr)
library(tidyr)
library(tmap)
library(leaflet)
library(spbabel)
library(geosphere)
library(htmlwidgets)
library(stringi)
library(stringdist)
library(raster)
options(scipen=999, digits = 2)

Como primera tarea tenemos el visualizar un mapa de Valencia delimitado por las zonas correspondientes a cada uno de los códigos postales. Colorear los polígonos por código postal y escribir dentro de cada zona el código postal correspondiente con un tamaño proporcional al área de dicha zona.

El ficheron cpval.geojson lo descargamos de la pagina http://www.cartociudad.es/portal/

Este contiene la información espacial de los códigos postales en el archivo con extensión .shp

direccion <- "Plaza del ayuntamiento, Valencia, Spain"
coords <- geocode(direccion)
coords

##     lon lat
## 1 -0.38  39

# cp <- readOGR("cpval.geojson")
cp <- readOGR("C:\\Users\\PCPCPCC\\Desktop\\Master\\VAD\\Geolocalizacion\\PracticaVisualizacionEspacial\\CARTO\\CODIGO_POSTAL.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\PCPCPCC\Desktop\Master\VAD\Geolocalizacion\PracticaVisualizacionEspacial\CARTO\CODIGO_POSTAL.shp", layer: "CODIGO_POSTAL"
## with 488 features
## It has 3 fields

plot(cp)

Podemos ver con la función plot() de R, el contenido de este archivo.

Prueba con leaftet en Valencia Capital.

Para representar solo Valencia Capital seleccionamos aquellos codigos postales que empiezan con 460.

Para mejorar la visualización eliminamos el CP 46012 que corresponde a la Albufera.

A continuación extraemos la información de los polígonos con sptable y obtenemos los centroides (en longitud y latitud) de los polígonos a partir de esta información, los cuales utilizaremos mas adelante.

cpval <- cp[grep("^460", cp$COD_POSTAL),]
cpval <- cpval[-15,]#Quitamos 46012 que se va muy al sur y queda fea la grafica!
geo <- sptable(cpval)[,c(1,5,6)] #Extraemos las coordenadas de los poligonos
means <- aggregate(geo[,2:3], geo[1], mean) #Obtenemos el centroide de los poligonos
colnames(means) <- c("ID","lat","lon")
means$ID <- cpval$COD_POSTAL

cpval@data$lat <- means$lat
cpval@data$lon <- means$lon

leaflet() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addPolygons(data=cpval,fillColor=c("red","green","blue","yellow"), fillOpacity = 0.2,weight = 2, stroke = TRUE,smoothFactor = 0.5, color = "black", opacity = 1) %>% addMarkers(cpval$lat,cpval$lon, label = as.character(cpval$COD_POSTAL),labelOptions = labelOptions(noHide = F)) %>% setView(coords$lon,coords$lat ,zoom=12)

Prueba con TMap en Valencia Capital.

areas <- numeric()
for (i in 1:length(cpval$ID_CP)){
    areas <- append(areas, cpval@polygons[[i]]@area)
} #Utilizamos el area de los poligonos para generar el area de las regiones.
cpval@data$areas <- areas
map <- tm_shape(cpval)+tm_fill("areas", legend.show = F,palette="RdYlGn")+tm_text("COD_POSTAL",size = "areas", legend.size.show = F)+tm_borders(alpha=0.8, col="black")
map

La siguiente tarea es: Dada una variable dirección, dibújala sobre el mapa, en la posición correspondiente, con un punto rojo.

Utilizamos la función geocode() del paquete ggmap para obtener las coordenadas de la dirección.

leaflet() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addPolygons(data=cpval,fillColor=c("red","green","blue","yellow","white","orchid"), fillOpacity = 0.2,weight = 2, stroke = TRUE,smoothFactor = 0.5, color = "black", opacity = 1) %>% addCircles(coords$lon,coords$lat, color="red", opacity=1, weight=10)

Extrae del mapa la zona (polígono) correspondiente a dicho código postal, y dibuja el mapa de dicha delimitación en otra figura, con el punto que indica la ubicación y el código postal superpuestos.

poligonos <- numeric()
poligonos <- array(dim=c(length(cpval),2))
for (i in 1:length(cpval@polygons)){
  for (j in 1:2){
    poligonos[i,j] <- cpval@polygons[[i]]@labpt[j]
    }
}

matdist <- distm(c(coords$lon,coords$lat), poligonos, fun=distVincentyEllipsoid)# Matriz de distancias entre el punto y los centroides de los poligonos.
posicion <- numeric(0)
posicion <- which.min(matdist)
poly <- rep(0,length(cpval$areas))
poly[posicion] <- 10
cpval@data$elegido <- poly

Usamos Leaftet.

customap <- leaflet() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addPolygons(data=cpval, fillColor=~ifelse(cpval$elegido == "10","green", "grey"), weight = 2, stroke = TRUE,smoothFactor = 0.5, color = "black", opacity = ~ifelse(cpval$elegido == "10",1, 0.2)) %>% addCircles(coords$lon,coords$lat, color="red", opacity=1, weight=10) %>% addLabelOnlyMarkers(coords$lon,coords$lat, label=cpval$COD_POSTAL[posicion], labelOptions = labelOptions(noHide = T))
customap

El siguiente apartado requería combinar información de la población de cada CP con su localización geográfica.

El fichero poblacion.csv se obtuvo de la pagina del INE, contiene la población en 2016 de cada CP de Valencia.

Las primeras lineas son para limpiar el .csv y darle un formato optimo.

Algo a tener en cuenta es que el CP que utiliza el INE no se corresponde con el CP “normal”. Por lo que debemos obtener estos CP a partir de queries a través de geocode() y revgeocode() que nos permite obtener información de un punto de coordenadas.

pobl <- readr::read_lines("C:\\Users\\PCPCPCC\\Desktop\\Master\\VAD\\Geolocalizacion\\PracticaVisualizacionEspacial\\poblacion.csv", skip= 7)
pobl <- pobl[1:266]#Las últimas lineas no contienen información relacionada.
pobl <- iconv(pobl, to="UTF-8", from="LATIN1") #Ajustar formato.
# pobl <- stri_trans_general(pobl,"Latin-ASCII") #Eliminar la mayoria de acentos.
pobl1 <- list()
for (i in 1:length(pobl)){
  pobl1 <- append(pobl1,paste(strsplit(pobl[i], " ")[[1]][-1],collapse=" "))
} #Quitamos el "Codigo Postal Falso"

localidad <- list()
poblacion <- list()
municipio <- list()
for (i in 1:length(pobl1)){
  localidad <- append(localidad,paste(strsplit(pobl1[[i]],",")[[1]][1]," Valencia"))#Agregamos la palabra Valencia para que el geocode sea mas fiable.
  poblacion <- append(poblacion,strsplit(pobl1[[i]],",")[[1]][length(strsplit(pobl1[[i]],",")[[1]])])
}
for (i in 1:length(pobl1)){
  municipio <- append(municipio,strsplit(pobl1[[i]],",")[[1]][1])
}

Obtenemos las coordenadas de las localidades.

# lonlat_sample <- data.frame()
# for (i in 1:length(localidad)){
#   lonlat_sample <- rbind(lonlat_sample,geocode(localidad[[i]]))
#   Sys.sleep(0.3) #Delay para evitar el limite de consultas por segundo
# }
# write.csv(file="coordenadas.csv",lonlat_sample)

Obtenemos el codigo postal a partir de las coordenadas.

El guardar y cargar los ficheros después de una consulta masiva es recomendable para no realizar las consultas cada vez que ejecutemos el código.

rm(lonlat_sample)

## Warning in rm(lonlat_sample): objeto 'lonlat_sample' no encontrado

lonlat_sample <- read.csv("C:\\Users\\PCPCPCC\\Desktop\\Master\\VAD\\Geolocalizacion\\PracticaVisualizacionEspacial\\coordenadas.csv", skip=1, header = F)
colnames(lonlat_sample) <- c("1","lon","lat")
lonlat_sample <- lonlat_sample[,2:3]
lonlat_sample$CP[1] <- 0 

#Agregamos la columna de "localidad" al dataframe
for (i in 1:nrow(lonlat_sample)){
  lonlat_sample$localidad[i] <- strsplit(localidad[[i]],"  ")[[1]][1]
}

# for (i in 1:nrow(lonlat_sample)){
#   vector <- c(lonlat_sample$lon[i],lonlat_sample$lat[i])
#   res <- revgeocode(vector, output = "more")
#   if (length(res$postal_code)>0){
#   lonlat_sample$CP[i] <- as.character(res$postal_code)
#   }
#   Sys.sleep(0.3)
# }
# write.csv(file="coords y cp.csv",lonlat_sample)#Guardamos para evitar hacer las consultas otra vez.

rm(lonlat_sample)
lonlat_sample <- read.csv("C:\\Users\\PCPCPCC\\Desktop\\Master\\VAD\\Geolocalizacion\\PracticaVisualizacionEspacial\\coords y cp.csv", skip=1, header = F, stringsAsFactors = F)
lonlat_sample <- lonlat_sample[2:ncol(lonlat_sample)]
colnames(lonlat_sample) <- c("lon","lat","CP","Localidad")

Hemos utilizado el fichero de 2016 para obtener los códigos postales y la población en un mismo data.frame. A continuación vamos a aprovechar la estructura de ese data.frame y verter la información de los últimos 20 años.

Obtenemos también del INE la información de la población para los últimos 20 años.

pobla <- readr::read_lines("C:\\Users\\PCPCPCC\\Desktop\\Master\\VAD\\Geolocalizacion\\PracticaVisualizacionEspacial\\poblacion años.csv", skip= 7)
pobla <- pobla[1:266]
# pobla <- iconv(pobla, to="UTF-8", from="LATIN1") #Ajustar formato.
# pobla <- stri_trans_general(pobla,"Latin-ASCII") #Eliminar la mayoria de acentos.
pobla1 <- list()
for (i in 1:length(pobla)){
  pobla1 <- append(pobla1,paste(strsplit(pobla[i], " ")[[1]][-1],collapse=" "))
} #Quitamos el "Codigo Postal Falso"

years <- list()
tablayears <- array(dim=c(length(pobla1),20))
colnames(tablayears) <- as.character(c(1997:2016))
for (i in 1:length(pobla1)){
  spl <- regexpr("[(0-9)]",pobla1[[i]])
  years <- append(years, substring(pobla1[[i]],spl,nchar(pobla1[[i]])-1)) #lista de años.
   for (j in 1:20){
     tablayears[i,j] <- strsplit(years[[i]],",")[[1]][j] #Matriz de años.
   }
}

df <- lonlat_sample[,c(1,2,3)]
for (i in 1:nrow(df)){
  df[i,4] <- municipio[[i]]
}
colnames(df)[4] <- "localidad"
df <- cbind(df, tablayears)
for (i in 5:ncol(df)){
df[,i] <- as.numeric(as.character(df[,i]))
}

df  <- df[!is.na(df[1]),]# quitamos NA's

A continuación unimos nuestros datos de población con el SpatialPolygonDataFrame de la provincia de Valencia.

Para esto utilizamos la función over() del paquete sp que nos permite realizar uniones entre diversos objetos de tipo spatial data como (SpatialPolygonsDataFrame, SpatialPoints, etc)

SpainRegions4 <- getData('GADM', country='ESP', level=4)
val <- SpainRegions4[SpainRegions4$NAME_2=="Valencia",]
# geo <- sptable(val)[,c(1,5,6)] #Extraemos las coordenadas de los poligonos
# means <- aggregate(geo[,2:3], geo[1], mean) #Obtenemos el centroide de los poligonos
# colnames(means) <- c("ID","lon","lat")
# means$ID <- val$ID_4
areas <- numeric()
for (i in 1:length(val@data$OBJECTID)){
    areas <- append(areas, val@polygons[[i]]@area)
} #Utilizamos el area de los poligonos para generar el area de las regiones.
val@data$areas <- areas

for (i in 1:ncol(tablayears)){  
  val@data[ncol(val@data)+1] <- rep(0,269)
  names(val@data)[ncol(val@data)] <- paste0("Pob",1996+i)
}

for (i in 1:nrow(df)){
  a <- SpatialPoints(df[i,c(1,2)],proj4string = val@proj4string)
  d <- over(a, val)
  val@data[val$NAME_4 == d$NAME_4,] <- cbind(d[1:18],df[i,5:24]/d$areas)#Densidad de población
}

Vemos la densidad de población en 2016

pal <- ~colorQuantile("Greens", Pob2016)(Pob2016)

leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  addPolygons(data=val,fillColor=pal,fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2016,0)), weight = 2, stroke =
              T,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = 1,
                                                                                        weight = 2,
                                                                                        fillOpacity = .8, bringToFront =T,
                                                                                        sendToBack = T))

Representa varios mapas indicando como se ha modificado la población de la provincia en los años con registros disponibles.

Para este apartado la idea del profesor era que utilizaramos la función facet_grid() de ggmap, pero leaflet nos permite hacer un mapa interactivo en el cual podemos representar cada año como podemos ver a continuación.

mappobyears <- leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2016)(Pob2016),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2016)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2016") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2015)(Pob2015),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", as.character(val$Pob2015)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2015") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2014)(Pob2014),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2014)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2014") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2013)(Pob2013),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2013)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2013") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2012)(Pob2012),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2012)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2012") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2011)(Pob2011),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2011)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2011") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2010)(Pob2010),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2010)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2010") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2009)(Pob2009),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2009)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2009") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2008)(Pob2008),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2008)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2008") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2007)(Pob2007),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2007)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2007") %>%
  addPolygons(data=val,fillColor=~colorQuantile("Greens", Pob2006)(Pob2006),fillOpacity = 0.7,
              label=stringr::str_c(val$NAME_4,": ", round(val$Pob2006)), weight = 2, stroke =
              TRUE,smoothFactor = 0.5, color = "black",highlightOptions = highlightOptions(color='#00ff00', opacity = .8,
                                                                                           weight = 2,
                                                                                           fillOpacity = .8, bringToFront
                                                                                           = TRUE, sendToBack = TRUE),
              group="2006") %>%
  setView(coords$lon, coords$lat, zoom=8) %>%
  addLayersControl(
    baseGroups=c("2016", "2015","2014","2013","2012","2011","2010","2009","2008","2007"),
    position = "topright",
    options = layersControlOptions(collapsed = FALSE)
  )
mappobyears

#saveWidget(widget=mappobyears,file="mapa poblacion years1.html")