4

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: enter image description 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?

4

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

enter image description here

  • Thank you! This looks good too. I have previously done maps with different packages and ended up having to transform it all to use ggplot, and it was a major headache for me (still pretty novice at R). But I do need the learning experience and this looks great, so I think I'll give it a go! – Lauren Dec 17 '15 at 15:32
  • I have tried using your code and do not get the same result that you've posted. The color scale and point size scale are plotted as 2 different legends. Additionally, the Ethiopia polygon over which the points are plotted is actually a custom filled in polygon in my latest code (check my latest question), and it won't let me use this code over top of that. Any idea why it won't work? – Lauren Dec 31 '15 at 20:54
  • See @hrbrmstr's comment in your followup question for an explanation and @MLavoie's comment for a potential solution. – cengel Jan 4 '16 at 2:12
  • Yes, I tried them! It worked in a way, but still plotted the legend incorrectly with the graduated circles and color scale being 2 separate legend items, as shown here. The same thing happened when I ran the exact code that you suggested, shown here. Should I post another question for this solution? – Lauren Jan 4 '16 at 17:46
  • @cengel sorry to necro your post - but is it possible to do proportional symbols for VARIABLE1 and graduated scale for VARIABLE2 on the same ggplot map? – Reputable Misnomer Mar 25 at 23:42
3

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

mapview_viewer

  • This looks great! I will certainly take a closer look at this. :) – Lauren Dec 16 '15 at 17:55
2

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)
  • Care to explain the key.entries portion of that code? – Lauren Dec 15 '15 at 21:24
  • @Lauren, I avoid lattice at all cost (which includes spplot and bubble). As such, I expanded my answer to accomplish this in base plot. – Jeffrey Evans Dec 15 '15 at 22:57
  • I'm having a hard time working my way through the first bit of code, but I'm certainly more familiar with the second option you suggested (reclassing the vector sizes/colors). Is there a way to use this same type of code structure but with real numbers? So instead of branches$share == .03, x, I would need to use branches$share <= 0.3, x. Unfortunately, this type of structure gives me an error saying that <= is not useful for factors. My data is listed as factors, but they should be considered real numbers. It appears that as.real is now defunct. – Lauren Dec 16 '15 at 17:54
  • If R is automatically coercing a column to a factor it is because of something contained in the data (eg., special character, quoted text). If the data is quoted and coercion is not because there is a special character in the column you can just defactorize the data using: df$x <- as.numeric(as.character(df$x)) – Jeffrey Evans Dec 16 '15 at 18:33
  • @Lauren, the mapview looks like it may be a more tractable route cartographically but, please work through the code I provided as it is what we commonly refer to as a "teachable moment" and will help bolster your understanding of R in general. – Jeffrey Evans Dec 16 '15 at 20:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.