Diagrama de barras e información geográfica con R
Diagrama de barras e información geográfica con R
¡Qué noche tan dichosa!
Sólo ella conoció el momento
en que Cristo resucitó del abismo. - Pregón Pascual 2021.
He aquí un ejemplo de información espacial y diagramas de barras que transmiten fácilmente la información más relevante de una base de datos geolocalizados.
Carlos Gil Bellosta acuñó el término huertograma para calificar los gráficos espaciales que confunden al lector. Aquí presentamos un gráfico que creemos muy informativo y sencillo de entender por parte de los lectores.
¿Qué información transmitir? Los datos no informan
Paula Casado actualiza en El Arte del Dato las enseñanzas de Edward Tufte. A ellos remitimos para quien desee aprender las técnicas básicas de cómo transmitir información estadística de forma visual.
Destaco que la labor de un especialista en datos consiste en condensar, en resumir, la información. Y los gráficos, cuanto más sencillos y comprensibles sean, mucho mejor. Desafortunadamente, percibo una marabunta de aplicaciones gráficas que delegan en el usuario todos los datos y éste es quién debe esforzarse por entenderlos.
Por ejemplo, las Shiny Apps, mediante un cuadro de mandos, permiten al usuario conocer los microdatos. En ese tipo de aplicaciones, es el usuario quien extrae las conclusiones. Ha de dedicar tiempo a manipular y extraer lo más relevante de la base de datos. Un ejemplo de este enfoque se muestra en el mapa interactivo de Antonio Díaz Muñiz sobre multas de tráfico en la ciudad de Gijón.
Desde el atrevimiento de mi creciente ignorancia, creo que la labor de un estadístico consiste en analizar los datos y convertirlos en conocimiento, incluyendo la incertidumbre si fuera necesario. Es decir, requiere un estudio previo que parece que se obvia en las visualizaciones actuales.
Inspirado por este gráfico sobre atascos de tráfico, creamos uno similar con R utilizando la base de datos del Ayuntamiento de Gijón y contamo de una forma distinta cómo se multa en esta ciudad.
Grid de R
El ecosistema de gráficos de R parece constreñirse al universo ggplot2. Muchos especialistas en gráficos desconocen que este popular paquete se basa en el sistema grid. Y cuando llega el momento de diseñar gráficos no estándar, el paquete gráfico lattice permite una interacción directa con el sistema grid que con ggplot2 no se consigue -o al menos este enano ignorante no ha descubierto cómo hacerlo de una forma sencilla-.
library(grid) # Grid library(lattice) # Top level functions for grid objects suppressPackageStartupMessages(library(data.table)) library(colorspace) # Palette of colors
Preprocesado de los datos
La parte tediosa consiste en preprocesar los datos. De una forma resumida, leemos los datos, observamos las diez calles con más multas y creamos una tabla con el resumen de esas localizaciones. Así dispondremos de dos bases de datos, que darán lugar a dos elementos gráficos, y que enlazaremos mediante segmentos superpuestos.
La primera base contiene los microdatos con las multas y sus ubicaciones geográficas. Su representación formará el contorno urbanístico de la ciudad de Gijón. Para que la escala del mapa sea manejable, eliminamos las ubicaciones muy lejanas.
datos <- fread(file.path("graphics","datawhour.csv"),dec=",") names(datos) <- tolower(names(datos)) ## Hour. Convert hour POSIXct and numeric variables. datos$hora <- as.POSIXct(datos$hour,format="%H:%M") datos$Hour <- as.numeric(datos$hora) ## Map ## To plot the map, we remove the outliers: In this way, the map shows with more precision the ## places where traffic violations happen. ## Since the asymmetry of the map, we remove the percentiles 5 in latitude, and 15 and 95 in longitude. datos$Latitude <- datos$latitude datos$Longitude <- datos$longitude limitslat <- quantile(datos$latitude,c(0.03,1),na.rm=TRUE) limitslong <- quantile(datos$longitude,c(0.01,0.999),na.rm=TRUE) f <- datos$latitude > limitslat[1] & datos$latitude < limitslat[2] & datos$longitude > limitslong[1] & datos$longitude < limitslong[2] datosgr <- na.omit(datos[f,]) ## Data set to plot the map datosgr[,.(place,latitude,longitude,fee)]
place latitude longitude fee
1: AVDA PRINCIPE DE ASTURIAS 43.53894 -5.692616 200
2: AVDA PRINCIPE DE ASTURIAS 43.53894 -5.692616 200
3: AVDA PRINCIPE DE ASTURIAS 43.53894 -5.692616 200
4: AVDA PRINCIPE DE ASTURIAS 43.53894 -5.692616 200
5: AVDA PRINCIPE DE ASTURIAS 43.53894 -5.692616 200
---
96001: AVDA OVIEDO 43.52270 -5.682149 100
96002: AVDA OVIEDO 43.52270 -5.682149 100
96003: AVDA OVIEDO 43.52270 -5.682149 100
96004: AVDA OVIEDO 43.52270 -5.682149 100
96005: AVDA OVIEDO 43.52270 -5.682149 100
Seleccionamos las diez calles donde se realizan más multas y las etiquetamos adecuadamente. Estos son los datos que se representarán como barras. Cada rectángulo se coloreará de forma distinta.
datoscalles <- datos[,.(fee=round(sum(fee)/(365*3))), by=place] datosdd <- datoscalles[order(-fee),][1:10] datosdd$place2 <- c("Avda. Oviedo", "Ctra. Siero", "Avda. Princ.", "Juan Carlos I", "Justo Castillo", "Ctra. Somio", "Gaspar Laviana", "Cimadevilla", "Avda. Llano", "Ctra. Carbonera") datosdd[,.(place2,fee)]
place2 fee
1: Avda. Oviedo 2699
2: Ctra. Siero 1350
3: Avda. Princ. 1223
4: Juan Carlos I 1166
5: Justo Castillo 915
6: Ctra. Somio 844
7: Gaspar Laviana 464
8: Cimadevilla 435
9: Avda. Llano 409
10: Ctra. Carbonera 286
Armonización de escalas
Queremos superponer dos gráficos. El primero representa un gráfico de barras y el segundo grafica puntos mediante la latitud y longitud. Ajustamos las medidas de este último a las dimensiones del primer gráfico mediante una escala lineal. Ese ajuste se ha realizado a mano, reiterando el proceso hasta obtener un diseño conjunto aceptable.
maxfee <- max(datosdd$fee) minfee <- min(datosdd$fee) minlat <- min(datosgr$latitude,na.rm=TRUE) maxlat <- max(datosgr$latitude,na.rm=TRUE) minlong <- min(datosgr$longitude,na.rm=TRUE) maxlong <- max(datosgr$longitude,na.rm=TRUE) datosgrdd <- datosgr[, .(latitude,longitude,place,fee)] datosgrdd$lat <- 7*(datosgrdd$latitude - min(datosgrdd$latitude)) / diff(range(datosgrdd$latitude)) + 1.5 datosgrdd$long <- (maxfee-minfee)*(datosgrdd$longitude - min(datosgrdd$longitude)) / diff(range(datosgrdd$longitude)) +500 datosgrddresumen <- datosgrdd[, .(latcenter=median(lat,na.rm=TRUE),longcenter=median(long,na.rm=TRUE)),by=place] datosgrddresumen <- datosgrddresumen[place %in% datosdd$place,] datosdd <- merge(datosdd,datosgrddresumen,all.x=TRUE) datosdd <- datosdd[order(-fee)] ## Cimadevilla data is missing. So we add the coordinates datosdd[8,latcenter:=7.2] datosdd[8,longcenter:=1600] datosdd$pos <- 10:1 datosdd[,.(place2,fee,latcenter,longcenter,pos)]
place2 fee latcenter longcenter pos
1: Avda. Oviedo 2699 1.788799 1004.2576 10
2: Ctra. Siero 1350 2.232996 1848.5991 9
3: Avda. Princ. 1223 5.252121 732.6217 8
4: Juan Carlos I 1166 5.676419 1098.6869 7
5: Justo Castillo 915 4.360528 2525.4773 6
6: Ctra. Somio 844 6.216697 2609.8363 5
7: Gaspar Laviana 464 2.549424 1510.5530 4
8: Cimadevilla 435 7.200000 1600.0000 3
9: Avda. Llano 409 2.884444 1643.4182 2
10: Ctra. Carbonera 286 1.734828 1447.3211 1
Color de cada barra
Como queremos que cada barra tenga un color distinto, modificamos la función lattice::panel.barchart para que asigne el color de la barra en función del valor de y. Se ha añadido un bucle for que dibuja cada barra panel.rect según la paleta de colores que se escoja, y se etiqueta cada barra para usarla si fuera necesario mediante Grob.
## Modifiy panel.barchart to get a name for each bar, and color. Just horizontal case. panel.barchart2 <- function (x, y, box.ratio = 1, box.width = box.ratio/(1 + box.ratio), horizontal = TRUE, origin = NULL, reference = TRUE, stack = FALSE, groups = NULL, col = if (is.null(groups)) plot.polygon$col else superpose.polygon$col, border = if (is.null(groups)) plot.polygon$border else superpose.polygon$border, lty = if (is.null(groups)) plot.polygon$lty else superpose.polygon$lty, lwd = if (is.null(groups)) plot.polygon$lwd else superpose.polygon$lwd, ..., identifier = "barchart") { plot.polygon <- trellis.par.get("plot.polygon") superpose.polygon <- trellis.par.get("superpose.polygon") reference.line <- trellis.par.get("reference.line") keep <- (function(x, y, groups, subscripts, ...) { !is.na(x) & !is.na(y) & if (is.null(groups)) TRUE else !is.na(groups[subscripts]) })(x = x, y = y, groups = groups, ...) if (!any(keep)) return() x <- as.numeric(x[keep]) y <- as.numeric(y[keep]) if (!is.null(groups)) { groupSub <- function(groups, subscripts, ...) groups[subscripts[keep]] if (!is.factor(groups)) groups <- factor(groups) nvals <- nlevels(groups) groups <- as.numeric(groupSub(groups, ...)) } if (horizontal) { if (is.null(groups)) { if (is.null(origin)) { origin <- current.panel.limits()$xlim[1] reference <- FALSE } height <- box.width if (reference) panel.abline(v = origin, col = reference.line$col, lty = reference.line$lty, lwd = reference.line$lwd, identifier = paste(identifier, "abline", sep = ".")) for( i in y){ ## print(i) ## print(x) panel.rect(x = origin, y = length(y)-i + 1, height = height, width = x[i] - origin, border = border, col = col[i], lty = lty, lwd = lwd, just = c("left", "centre"), identifier = paste0(identifier,".",i)) } } } }
La paleta de color que usamos viene definida en R:
MyPalette <- sequential_hcl(10,"Heat") # used in barchat
Dibujando los gráficos
Ahora se dibujan bastantes elementos de una tacada: el gráfico de barras (panel.barchart2), el número en cada barra (panel.text), los puntos geolocalizados en color gris (panel.points) y los segmentos que unen las barras con los centros de las ubicaciones (panel.segments). Estos segmentos se han diseñado de tal forma que las líneas se corten lo menos posible.
graficobarras <- barchart(reorder(place2,fee)~fee,data=datosdd, xlab=list(label=""), main="Recaudación de 21 000 euros diarios en multas de tráfico", xlab.top =list(label=paste("Las 10 principales calles donde se multa"), just="center",font=2, x = grid::unit(55, "mm")), scales = list(x = list(draw = FALSE), y=list(font=2), cex=c(1.4,1.4)), par.settings = list(axis.line = list(col = "transparent"), par.main.text = list(font = 2, # make it bold just = "left", cex=2., x = grid::unit(5, "mm")), par.sub.text = list(font = 2, cex=1.3, just = "left", x = grid::unit(5, "mm"))), panel=function(...) { args <- list(...) panel.barchart2(...,border=NULL,col=MyPalette) panel.points(x=datosgrdd[["long"]],y=datosgrdd[["lat"]],pch=".",size=0.03,alpha=0.3,col="grey") panel.text(args$x / 1.5, args$y, paste0(args$x),col="white", offset=0,font=2,cex=1.4) ## Oviedo i <- 1 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- x1+100 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) # Carbonera i <- 10 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) # Llano i <- 9 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) ## panel.segments(x1,y1,x2,y2) # Cimadevilla i <- 8 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) # Gaspar i <- 7 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- x1 + 100 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) # Somio i <- 6 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) ## panel.segments(x1,y1,x2,y2) # Gaspar i <- 5 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- x1 + 100 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) ## panel.segments(x1,y1,x2,y2) # Juan Carlos i <- 4 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- x1 + 100 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) # Av prin i <- 3 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- x1 + 100 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) i <- 2 x1 <- datosdd[["fee"]][i] y1 <- datosdd[["pos"]][i] x2 <- 2750 y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- x1 y2 <- datosdd[["latcenter"]][i] panel.segments(x1,y1,x2,y2,col=MyPalette[i]) x1 <- x2 y1 <- y2 x2 <- datosdd[["longcenter"]][i] y2 <- y1 panel.segments(x1,y1,x2,y2,col=MyPalette[i]) panel.points(x2,y2,col=MyPalette[i], pch=16,cex=2) } )
Salvamos el gráfico, respetando la proporción de ancho y alto del mapa.
png(file.path("graphics","grmapabarrasmultas.png"),width=830,height=580) print(graficobarras) dev.off()