结合ggplot和ggmap中的choropleth

Created a choropleth using ggplot2. Here's the ggplot code

使用ggplot2创建了一个等值区。这是ggplot代码

okc <- ggplot() +
  geom_polygon(data = mapdata, aes(x = long, y = lat, group = group,
                                   fill = B19013_001), color = "black", size = 0.5)+
  scale_fill_distiller(palette = "Reds", labels = comma,
                       breaks = pretty_breaks(n = 10), values = c(1,0)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme_nothing(legend = TRUE) +
  ggtitle('Map of 40109') 

Here's a sample of the data from mapdata:

以下是mapdata中的数据示例:

       long      lat order  hole piece         group          id
1 -97.54285 35.51951     1 FALSE     1 40109100100.1 40109100100
2 -97.54282 35.51954     2 FALSE     1 40109100100.1 40109100100
3 -97.54280 35.51963     3 FALSE     1 40109100100.1 40109100100
4 -97.54276 35.51976     4 FALSE     1 40109100100.1 40109100100
5 -97.54270 35.51993     5 FALSE     1 40109100100.1 40109100100
6 -97.54266 35.52016     6 FALSE     1 40109100100.1 40109100100
                                          NAME state county  tract B19013_001
1 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440
2 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440
3 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440
4 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440
5 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440
6 Census Tract 1001, Oklahoma County, Oklahoma    40    109 100100      33440

It produced this plot.

它产生了这个情节。

I also created a roadway map using ggmap. Here's the code:

我还使用ggmap创建了一个道路地图。这是代码:

map <- get_map(location = c(lon = mean(mapdata$lon), lat = mean(mapdata$lat))
               , zoom = 10
               , maptype = "roadmap"
               , color = "bw")
p <- ggmap(map) +
  scale_x_continuous(limits = c(min(mapdata$lon), max(mapdata$lon)), expand = c(0, 0)) +
  scale_y_continuous(limits = c(min(mapdata$lat), max(mapdata$lat)), expand = c(0, 0))
p

And here's the map it produces.

这是它产生的地图。

When I try to combine them though I get an error. Here's the code I use to combine them and the error:

当我尝试将它们组合起来时,虽然我收到了错误。这是我用来组合它们和错误的代码:

okc <- okc + p

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

I'm not sure why I'm getting this error. Is it because the maps are not scaled the same? I could not figure how else to scale ggmap other than using the very imprecise zoom function. If anyone has any ideas about how to layer the choropleth on top of the ggmap I will be very thankful.

我不知道为什么我会收到这个错误。是因为地图缩放不一样吗?除了使用非常不精确的缩放功能之外,我无法想象如何缩放ggmap。如果有人对如何在ggmap上叠加等值线有任何想法,我将非常感激。

Here's the rest of the code to recreate the ggplot choropleth.

这是重新创建ggplot等值区的其余代码。

    library(acs)
    library(ggplot2)
    library(ggmap)
    library(UScensus2010)
    library(RColorBrewer)
    library(dplyr)
    library(scales)

    #http://api.census.gov/data/key_signup.html
    api.key.install(key="c369cd6ed053a84332caa62301eb8afe98bed825")

    # Load in Shape File (You'll need to download this file from the census)
    #ftp://ftp2.census.gov/geo/tiger/TIGER2013/TRACT/tl_2013_40_tract.zip

    ## load, subset shapefile
    geodat<-readShapePoly("insert shapefile here", proj4string=CRS('+proj=longlat +datum=NAD83'))
    geodat<-geodat[geodat$COUNTYFP==109,]

    ## fortify for ggplot digestion
    geodat.f<-fortify(geodat,region="GEOID")

    # American Community Survey Data: Median HH Income for OK Census Tracts
    ok.counties=geo.make(state="OK", county="Oklahoma", tract="*")
    ok.income<-acs.fetch(geography=ok.counties, table.number="B19013", endyear=2013)


    # Merge Data Sets 
    geo_dat<-geography(ok.income)
    var_dat<-as.data.frame(estimate(ok.income))
    acs_data<-cbind(geo_dat,var_dat)
    acs_data$id<- paste("40109", acs_data$tract, sep = "")

    ## from dplyr
    mapdata<-left_join(geodat.f,acs_data)

    okc <- ggplot() +
      geom_polygon(data = mapdata, aes(x = long, y = lat, group = group,
                                       fill = B19013_001), color = "black", size = 0.5)+
      scale_fill_distiller(palette = "Reds", labels = comma,
                           breaks = pretty_breaks(n = 10), values = c(1,0)) +
      guides(fill = guide_legend(reverse = TRUE)) +
      theme_nothing(legend = TRUE) +
      ggtitle('Map of OKC')

1 个解决方案

#1

This is actually much better done in Leaflet. It looks better aesthetically and also code wise it is more intuitive.

这实际上在Leaflet中做得更好。它在美学上看起来更好,而且代码也更直观。

library(leaflet)
library(rgdal)
library(RColorBrewer)

pal <- colorNumeric("OrRd", domain = new$pct_minority_popn)

leaflet(mapdata) %>%
 addTiles %>%
 addPolygons(stroke=T, fillOpacity=.5, smoothFactor=.5, color=~pal(B19013_001)) %>%
 addLegend("bottomright", pal=pal, values=~B19013_001, title="Legend Title", opacity=.8)

You can change the bottom map by replacing the addTiles command with something like addProviderTiles("CartoDB.Positron"). You can see the rest of the options and more info on leaflet at: https://rstudio.github.io/leaflet/basemaps.html

您可以通过将addTiles命令替换为addProviderTiles(“CartoDB.Positron”)来更改底部映射。您可以在以下网址查看其他选项和有关传单的更多信息:https://rstudio.github.io/leaflet/basemaps.html