How to overlay points on polygons in ggplot2?

Question

I am struggling to plot a city spatial points data frame over a shape file showing average household income in different census tracts. My income data is downloaded from CDPHE Open Data for Colorado and I'm using the city data available in the maps package. I specifically have to use ggplot2 to visualize the data. I read through similar questions and modified code from other answers, but still can't get it.

One of the ways I have tried to code this is below:

library(ggplot2)
library(maps)
library(sp)
library(raster)
library(spdplyr)

incp_prj <- shapefile("Income_Poverty_(Census_Tracts).shp")

data(us.cities)
coords <- cbind(us.cities$long, us.cities$lat)
us.cities <- SpatialPointsDataFrame(coords = coords, data = us.cities, proj = incp_prj@proj4string)
co.cities <- us.cities %>% filter(country.etc == "CO")

pt_data = as.data.frame(incp_prj)
grid_data = as.data.frame(co.cities)

ggplot(grid_data, aes(x = long, y = lat)) + geom_tile(aes(fill = incp_prj$Poverty_Me)) + 
                            geom_point(data = pt_data)

Which returns the error:

Aesthetics must be either length 1 or the same as the data (21): fill

If you are interested in downloading my specific income data you can find it here.

Answer 1

After reading more examples from the R graph gallery I was able to figure it out! I needed to combine the numerical information about poverty to the polygon data frame even though they were the same spatial object.

Here is the code to graph the points over polygons if anyone is interested:

library(ggplot2)
library(maps)
library(sp)
library(raster)
library(spdplyr)
library(broom)
library(mapproj)
library(viridis)

#Read in shape file and pull cities from {maps} package
incp_prj <- shapefile("Income_Poverty_(Census_Tracts).shp")

data(us.cities)
coords <- cbind(us.cities$long, us.cities$lat)
us.cities <- SpatialPointsDataFrame(coords = coords, data = us.cities, proj = incp_prj@proj4string)
co.cities <- us.cities %>% filter(country.etc == "CO")
grid_data = as.data.frame(co.cities)

#Fortify data to put into a dataframe that ggplot can use
spdf_fortified <- tidy(incp_prj, region = "OBJECTID") 

#coerce x and y IDs to match
incp_prj %>% mutate(`OBJECTID`= as.character(`OBJECTID`))

#Combine numerical data(income) with the polygons
graph <- spdf_fortified %>%
  left_join(. , incp_prj %>% mutate(`OBJECTID`= as.character(`OBJECTID`)), by = c("id"="OBJECTID"), copy = TRUE)

#Plot
ggplot() +
  geom_polygon(data = graph, aes(fill = Poverty_Me, x = long, y = lat, group = group)) +
  theme_void() +
   scale_fill_viridis(trans = "log", name="Income (USD)", guide = guide_legend( keyheight = unit(3, units = "mm"), keywidth=unit(12, units = "mm"), label.position = "bottom", title.position = 'top', nrow=1) ) +
  labs(title = "Average Household Income and Major Cities in Colorado") +
  theme(
    text = element_text(color = "#22211d"),
    plot.background = element_rect(fill = "#f5f5f2", color = NA),
    panel.background = element_rect(fill = "#f5f5f2", color = NA),
    legend.background = element_rect(fill = "#f5f5f2", color = NA),

    plot.title = element_text(size= 16, hjust=0.01, color = "#4e4d47", margin = margin(b = -0.1, t = 0.4, l = 2, unit = "cm")),
    legend.position = c(0.7, 0.09)) +
  coord_map()+
#Add the city points
  geom_point(data = grid_data, aes(x = grid_data$long, y = grid_data$lat), size = 3, 
        shape = 23, fill = "orange")

Which produced this: