Overlaying ggplot, geom_polygon choropleth on ggmap giving errors

Question

I am struggling to overlay a choropleth on a map using ggmap, ggplot, and geom_polygon. I am very new to R, so I am not sure how to proceed. I'm using R Studio. I've uploaded my data (which is just fake test data), shapefiles, and R code here

I have a data file of subject counts within census dissemination areas for Alberta (actually just random numbers for testing). I also have a shapefile for these areas. I am only interested in the City of Edmonton area, and it's difficult to see what part of the city we're looking at without an underlying roadmap.

My code is below:

#################################
#enter long and lat for appx area of City of Edmonton (this is so we can change it easier later to narrow down if we want)
lat1 <- 53.65 
lat2 <- 53.44 
long1 <- -113.68 
long2 <- -113.35
#################################

xcoords <- c(long2, long1) 
ycoords <- c(lat2, lat1) 

#load data 
data <-read.csv("C:.../fakedata.csv") #data file location 
names(data) 
data

#Load in the shape file 
shp=readShapeSpatial('C:.../gda_048b06a_e.shp') #shape file location

plot(shp) 

Which gives this (as expected), all the DAs for the province. However, this is too big, and there are too many DAs for me to manually go through and see which ones are in Edmonton and which aren't so I decided to colour them in according to # of subjects, and restrict to particular long and lat:

#Making the data types match 
data$amount <-as.numeric(as.character(data$amount))

is.factor(shp$DAUID) 
data$DAUID<-as.factor(data$DAUID)
is.factor(data$DAUID) #should be True 
is.factor(shp$DAUID) #should be True

shp@data <- left_join(shp@data, data) 
head(shp@data) 
map2 <- fortify(shp, region = "DAUID")

#wait a bit for the fortify to work!
map2 <- rename(map2, DAUID = id) 
map2 <- left_join(map2, shp@data) 
ggplot2 <- ggplot(map2)

#DA map 
datamap <- ggplot2 + geom_polygon(aes(long, lat, group = group, fill = amount),color = "grey") +
scale_fill_gradient(low='white', high='red') + coord_map(xlim =
xcoords,ylim = ycoords) 
datamap

datamap is also what I expect it to look like.

#map of Edm 
edmonton <- get_map(location = 'edmonton',maptype="road") 
ggmap(edmonton)

What I want to do is overlay my DA map onto the road map, so we can get a better idea of where the dark vs light areas are. It's at this point I begin to receive errors.

#stuck on how to overlay datamap on top of ggmap(edmonton)

datamap2 <- ggmap(edmonton) + datamap 
datamap2

I receive this error:

Error in p + o : non-numeric argument to binary operator In addition: Warning message: Incompatible methods ("+.gg", "Ops.data.frame") for "+"

I've played around with it for a week trying to incorporate solutions I've seen to similar problems but I can't seem to get it to work. Since my shapefile is so much bigger than my city map, could that be the issue? I've never done anything like this in R before and GIS is outside of my skillset.

My question is similar to this previous question, but the solutions from that did not work for me. I apologize if it's too close in content to other questions, and appreciate any guidance!

(I also apologize if this is too long! I didn't want to leave out any useful information)

Answer 1

I looked into your case and ended up doing the work in my way. I hope you do not mind that. What you needed was to have a data set with amount in the fake data set in your polygon data. I guess that is why you were probably using left_join. Anyway, each polygon had different numbers of data points. I calculated the total of data points for each polygon and created num. Using this, you can create a vector with DAUID and replace id in mymap2. Once you arrange your data, you get a raster map with ggmap. Then, you draw polygons on top of the raster map. You ca control alpha value and see how colors come out. If you lower the value to 0.1, for instance, it is hard to see difference among colors. You may want to consider another way to show roads. One way would be to use Canadian road data set and draw some major roads in Edmonton. Hope this will help you.

library(readr)
library(dplyr)
library(ggplot2)
library(ggmap)
library(rgdal)

#load data
mydata <- read_csv("fakedata.csv") %>%
          rename(id = DAUID)

#Load in the shape file
mymap <- readOGR(dsn = ".", layer = "gda_048b06a_e")

mymap@data$DAUID <- as.numeric(as.character(mymap@data$DAUID))

mymap2 <- fortify(mymap)

# Get the number of data points for each polygon
group_by(mymap2, as.numeric(id)) %>%
summarize(total = n()) -> num

# Replace id with DAUID, and merge with mydata
mymap2 %>%
mutate(id = rep(unique(mymap$DAUID), num$total)) %>%
left_join(mydata, by = "id") -> mymap2


#map of Edm
edmonton <- get_map(location = "edmonton",maptype = "road")

ggmap(edmonton) + 
geom_map(data = mymap2, map = mymap2,
         aes(x = long, y = lat, group = group, map_id = id, fill = amount),
         color = "black", size = 0.2, alpha = 0.3) +
scale_fill_gradient(low = "white", high = "red") +
coord_map(xlim = c(-113.35, -113.68), ylim = c(53.44, 53.65))