¿Cómo trazo puntos como círculos graduados/proporcionales en R?

Question

¡Estoy de vuelta con más preguntas sobre mapeo en R! Esta se refiere a cómo simbolizar puntos por un cierto atributo con círculos graduados o proporcionales.

DATOS: Puede encontrar todo mi código aquí si estás interesado. Mi archivo CSV de interés está aquí. El shapefile de Etiopía se puede descargar aquí

He graficado puntos del CSV sobre un shapefile. Para ser breve, aquí está la versión muy resumida de mi código:

library(raster)
library(rgdal)

#establecer el directorio de trabajo

#leer el shapefile de Etiopía
eth <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "ETH_adm0")

#leer la ubicación total de las sucursales desde el csv
branches <- read.csv("Branches_Africa.csv", header = TRUE)

#transformar el data frame de las sucursales a un spdf para el mapeo
coordinates(branches) <- ~Lon + Lat

#graficar el shapefile con superposición de puntos
plot(eth)
points(branches$Lon, branches$Lat, col = "blue", pch = 16, cex = .5)

Como puedes ver en el CSV de sucursales, hay varios atributos diferentes. En mi caso, me gustaría simbolizar el atributo "share" con círculos graduados o proporcionales, como se ve aquí:

idealmente, también me gustaría cambiar los colores de cada punto para que coincidan lo más posible con la simbología de QGIS.

Si miras mi código completo, notarás que no he usado ggplot para el mapeo, por lo que preferiría que la solución no venga de ese paquete (¿o acaso ggplot facilitaría las cosas una vez que haya transformado todos los datos para que se ajusten al paquete?). Intenté usar spplot según este post con:

spplot(branches, "share", col.regions = brewer.pal(9, "Reds"), cuts = 8, scales = list(draw = T))

pero el código me dio el siguiente error: Error in fill.call.groups(args.xyplot, z = z, edge.col = edge.col, colorkey = colorkey, : number of colors smaller than number of factor levels

También intenté trabajar con mapCircles del paquete rCarto, pero tampoco tuve suerte allí.

ENTONCES, ¿cómo puedo cambiar la simbología de branches$share para que se muestren círculos graduados o proporcionales (y preferiblemente de colores diferentes) en lugar de solo un tamaño/color?

Answer 1

Sé que prefieres que la solución no use el paquete ggplot, pero pensé que podría ser útil agregar un ejemplo rápido, solo para comparar.

library(rgdal)
library(ggmap)
library(scales)

#datos
eth <- readOGR(".", "ETH_adm0")
branches <- read.csv("Branches_Africa.csv", header = T)
branches <- branches[branches$CO == "ET", ] # seleccionar subconjunto de Etiopía
branches$share <- as.numeric(as.character(branches$share)) # convertir a números

# fortificar el límite de Etiopía para ggplot
eth_df <- fortify(eth)
# obtener un mapa base
eth_basemap <- get_map(coordinates(eth), zoom = 5)

# trazar
ggmap(eth_basemap) + 
  geom_polygon(data=eth_df, aes(x=long, y=lat, group=group), fill="red", alpha=0.1) +
  geom_point(data=branches, aes(x=Lon, y=Lat, size=share, fill=share), shape=21, alpha=0.8) +
  scale_size_continuous(range = c(2, 9), breaks=pretty_breaks(7)) +
  scale_fill_distiller(breaks = pretty_breaks(7)) 

Answer 2

Si buscas una interfaz más interactiva, deberías echar un vistazo al paquete mapview que se ha lanzado en CRAN recientemente. El paquete también cuenta con un método spplot para objetos 'mapview', pero este último aún no admite objetos de múltiples capas. Si buscas una introducción corta, simplemente echa un vistazo al documento del paquete.

## paquetes requeridos
# install.packages("mapview")
library(mapview)
library(sp)
library(RColorBrewer)

### polígono

## obtener datos
eth <- raster::getData(country = "ETH", level = 0)

### puntos

## importar datos
branches <- read.csv("Branches_Africa.csv", header = TRUE)
coordinates(branches) <- ~ Lon + Lat
proj4string(branches) <- proj4string(eth)

## seleccionar puntos de Etiopía
branches <- branches[eth, ]

### visualizar

## colores
cols <- colorRampPalette(brewer.pal(9, "Blues"))

## abrir visor
mapview(branches, zcol = "share", cex = "share", color = "black",
        fillColor = cols(100), alpha.regions = .8) + eth

Answer 3

Puedes mapear círculos graduados usando la función bubble en sp.

library(sp)
data(meuse)
coordinates(meuse) <- c("x", "y")

bubble(meuse, "cadmium", main = "concentraciones de cadmio (ppm)", 
       key.entries = 2^(-1:4), col = c("blue","red"))

Ahora, desglosemos la función para proporcionar más flexibilidad usando la creación de gráficos básicos.

Primero, creamos un vector que se pasará al argumento "cex" en plot. Esto definirá el tamaño de nuestros círculos.

zcol = "cadmium"
q = sqrt( abs(quantile(meuse@data[,zcol])) )
az = sqrt( abs(meuse@data[,zcol]) )
maxsize = 4
pt.cex <- as.vector(maxsize * az / max(az, q))

Ahora, usando el paquete classInt para encontrar intervalos y asignar colores, podemos crear un vector que se puede pasar al argumento "col" en plot.

library(classInt)    
pal <- c("blue", "green", "red3")
    pt.class <- classIntervals(meuse@data[,zcol], n=10, style="hclust", method="complete")
    pt.col <- findColours(pt.class, pal)

Ahora podemos juntarlo todo representando las coordenadas x e y de nuestro objeto punto y pasando nuestros vectores que definen el tamaño y color del punto a los argumentos "cex" y "col". La leyenda se construye extrayendo información de los atributos del objeto classInt.

plot(coordinates(meuse)[, 1], coordinates(meuse)[, 2], asp = 1, 
    cex = pt.cex, col = pt.col, pch = 20) 
    legend("topleft", fill = attr(pt.col, "palette"),
           legend = names(attr(pt.col, "table")), bg="white")
           pt.tab <- attr(pt.col, "table")
           legtext <- paste(names(pt.tab), " (", pt.tab, ")", sep="")

Si tienes una variable nominal que estás representando, una alternativa es cambiar simplemente el vector a los tamaños y colores que desees. Aquí tienes un ejemplo usando el atributo "soil", que tiene tres niveles.

soil.cex <- ifelse(meuse$soil == 1, 0.75, 
              ifelse(meuse$soil == 2, 1.5,
                ifelse(meuse$soil == 3, 2.5, NA))) 
soil.col <- ifelse(meuse$soil == 1, "blue", 
              ifelse(meuse$soil == 2, "green",
                ifelse(meuse$soil == 3, "red", NA))) 

Representa los tamaños y colores nominales para el suelo.

plot(coordinates(meuse)[, 1], coordinates(meuse)[, 2], asp = 1, 
    cex = soil.cex, col = soil.col, pch = 20) 
    legend("topleft", legend = c("tipo de suelo 1", "tipo de suelo 2", "tipo de suelo 3"),
           pch=c(20,20,20), col = c("blue","green","red"), pt.cex = c(0.75,1.5,2.5))

La función bubble en realidad utiliza la función xyplot de lattice. Más información sobre el uso de xyplot se puede encontrar en la ayuda de las funciones invocadas utilizando ?xyplot. El argumento key.space define los atributos de etiquetado/leyenda. Aquí tienes un ejemplo muy simple usando los datos que creamos anteriormente. Aunque opté por no crear una leyenda.

library(lattice)    
xyplot(coordinates(meuse)[, 2] ~ coordinates(meuse)[, 1], col = pt.col, 
           cex = pt.cex, pch = 20)