How do I plot points as graduated/proportional circles in R?

Question

I'm back again with more R mapping questions! This one is concerning symbolizing points by a certain attribute with graduated/proportional circles.

DATA: My full code can be found here if you're interested. My CSV of interest is here. Ethiopia shapefile can be downloaded here

I've plotted points from the CSV over a shapefile. For the sake of brevity, here's the heavily abridged version of my code:

library(raster)
library(rgdal)

#set your working directory

#read in eth shapefile
eth <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "ETH_adm0")

#read in total branch location csv
branches <- read.csv("Branches_Africa.csv", header = TRUE)

#transform branches data frame to spdf for mapping
coordinates(branches) <- ~Lon + Lat

#plot shapefile with point overlay
plot(eth)
points(branches$Lon, branches$Lat, col = "blue", pch = 16, cex = .5)

As you can see in the branches CSV, there are several different attributes. In my case, I would like to symbolize the "share" attribute with graduated/proportional circles, as seen here:

Ideally I would also like to change the colors of each point as well to match the QGIS symbology as closely as possible.

If you look at my full code, you'll notice that I've not used ggplot to map, so I would prefer the solution not come from that package (or would ggplot make things easier once I've transformed all the data to fit the package?). I tried using spplot as per this post with:

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

but the code gave me the following 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

I also tried to work with mapCircles from the rCarto package, but I had no luck there either.

SO, how can I change the symbology of branches$share so that graduated/proportional (and preferably differently colored) circles are displayed instead of just one size/color?

Answer 1

I know you prefer the solution not to use the ggplot package, but I thought it might be helpful to add a quick example, just for comparison.

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

#data
eth <- readOGR(".", "ETH_adm0")
branches <- read.csv("Branches_Africa.csv", header = T)
branches <- branches[branches$CO == "ET", ] # select Ethiopia subset
branches$share <- as.numeric(as.character(branches$share)) # convert to numbers

# fortify Ethiopia boundary for ggplot
eth_df <- fortify(eth)
# get a basemap
eth_basemap <- get_map(coordinates(eth), zoom = 5)

# plot
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

If you are looking for a more interactive interface, you should possibly have a look at the mapview package which has been released on CRAN only recently. The package also features a spplot method for 'mapview' objects, but the latter does not support multi-layer objects yet. If you are looking for a short introduction, just have a look at the package vignette.

## required packages
# install.packages("mapview")
library(mapview)
library(sp)
library(RColorBrewer)


### polygon

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


### points

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

## select points from ethiopia
branches <- branches[eth, ]


### visualize

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

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

Answer 3

You can map graduated circles using the bubble function in sp.

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

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

Now, let's break down the function to provide more flexibility using base plotting.

First, we create a vector that will be passed to the "cex" argument in plot. This will define the size of our circles.

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))

Now, using the classInt package to find breaks and assign colors, we can create a vector that can be passed to the "col" argument in 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)

Now we can put it all together by plotting the x and y coordinates of our point object and passing our vectors defining point size and color to the "cex" and "col" arguments. The legend is constructed by pulling information from the classInt object(s) attributes.

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="") 

If you have a factor or nominal variable that you are plotting, an alternative is to just reclass the vector to the the sizes and colors you would like. Here is an example using the the "soil" attribute, which has three levels.

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))) 

Plot the nominal sizes and colors for soil.

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

The bubble function actually uses the xyplot function from lattice. More information on using xyplot can be found in the functions help invoked using ?xyplot. The key.space argument defines the labeling/legend attributes. Here is a very simple example using the data we created above. Although, I opted to not create a legend.

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