使用ggplot/mapping在R中绘制美国县的问题——可视化形状方面的问题

Question

使用ggplot/mapping在R中绘制美国县的问题——可视化形状方面的问题

9

我有一个在R中的数据框，名为obesity_map，它基本上给出了每个县的州、县和肥胖率。它看起来更或多或少像这样：

obesity_map = data.frame(state, county, obesity_rate)

我正在尝试通过以下方式在地图上显示美国各县的肥胖率，以便进行可视化：

us.state.map <- map_data('state')
head(us.state.map)
states <- levels(as.factor(us.state.map$region))
df <- data.frame(region = states, value = runif(length(states), min=0, max=100),stringsAsFactors = FALSE)

map.data <- merge(us.state.map, df, by='region', all=T)
map.data <- map.data[order(map.data$order),]
head(map.data)

map.county <- map_data('county')
county.obesity <- data.frame(region = obesity_map$state, subregion = obesity_map$county, value = obesity_map$obesity_rate)
map.county <- merge(county.obesity, map.county, all=TRUE)
ggplot(map.county, aes(x = long, y = lat, group=group, fill=as.factor(value))) + geom_polygon(colour = "white", size = 0.1)

它基本上创建了一个看起来像这样的图像：

正如您所见，美国被分成奇怪的形状，颜色不是一致的渐变色，您无法从中获得太多信息。但我真正想要的是以下内容，但每个县都填充了：

我对此还相当新，因此我将感激任何和所有的帮助！

编辑：

这是dput的输出：

dput(obesity_map)

structure(list(X = 1:3141, FIPS = c(1L, 3L, 5L, 7L, 9L, 11L, 
13L, 15L, 17L, 19L, 21L, 23L, 25L, 27L, 29L, 31L, 33L, 35L, 37L, 
39L, 41L, 43L, 45L, 47L, 49L, 51L, 53L, 55L, 57L, 59L, 61L, 63L, 
65L, 67L, 69L, 71L, 73L, 75L, 77L, 79L, 81L, 83L, 85L, 87L, 89L, 
91L, 93L, 95L, 97L, 99L, 101L, 103L, 105L, 107L, 109L, 111L, 
113L, 115L, 117L, 119L, 121L, 123L, 125L, 127L, 129L, 131L, 133L, 
13L, 16L, 20L, 50L, 60L, 68L, 70L, 90L, 100L, 110L, 122L, 130L, 
150L, 164L, 170L, 180L, 185L, 188L, 201L, 220L, 232L, 240L, 261L, 
270L, 280L, 282L, 290L, 1L, 3L, 5L, 7L, 9L, 11L, 12L, 13L, 15L, 
17L, 19L, 21L, 23L, 25L, 27L, 1L, 3L, 5L, 7L, 9L, 11L, 13L, 15L, 
17L, 19L, 21L, 23L, 25L, 27L, 29L, 31L, 33L, 35L, 37L, 39L, 41L,

由于涉及每个美国县，所以这是一个庞大的数字数量，我缩写了结果并放在了前几行。

基本上，数据框看起来像这样:

print(head(obesity_map))


  X FIPS state_names county_names obesity
1 1    1     Alabama      Autauga    24.5
2 2    3     Alabama      Baldwin    23.6
3 3    5     Alabama      Barbour    25.6
4 4    7     Alabama         Bibb     0.0
5 5    9     Alabama       Blount    24.2
6 6   11     Alabama      Bullock     0.0

我尝试按照示例使用 ggcounty，但一直出现错误。我不确定我哪里做错了：

library(ggcounty)

# breaks
obesity_map$obese <- cut(obesity_map$obesity, 
                  breaks=c(0, 5, 10, 15, 20, 25, 30), 
                  labels=c("1", "2", "3", "4", 
                           "5", "6"),
                  include.lowest=TRUE)

# get the US counties map (lower 48)
us <- ggcounty.us()

# start the plot with our base map
gg <- us$g

# add a new geom with our population (choropleth)
gg <- gg + geom_map(data=obesity_map, map=us$map,
                aes(map_id=FIPS, fill=obesity_map$obese), 
                color="white", size=0.125)

但是我总是收到一个错误，说：“错误：参数必须可强制转换为非负整数”

有什么想法吗？再次感谢您的所有帮助！我非常感激。

- user3648073

我不确定我是否理解你的意思？数据框中的数据是每个美国县的平均值，因此我正在尝试在地图上显示每个比率。例如，如果洛杉矶县的肥胖率为40％，橙县的肥胖率为50％，则会显示两个不同的县以不同的渐变方式着色其百分比。如果我的原始帖子缺乏清晰度，我深表歉意！ - user3648073

我明白了。你使用了“每个县的肥胖率不同”。你能添加一些数据样本吗？ - Paulo E. Cardoso

@PauloCardoso，没问题！对于可怕的格式化我很抱歉。还在摸索中！`%%Robesity = read.csv('./data/cleaned/obesitymap.csv', header = T)`print(head(obesity))

  X state_names county_names obesity 1 1     阿拉巴马州      奥托加县    24.5 2 2     阿拉巴马州      鲍尔德温县    23.6 3 3     阿拉巴马州      巴伯县    25.6 4 4     阿拉巴马州         比布县 -1111.1 5 5     阿拉巴马州       布朗特县    24.2 6 6     阿拉巴马州      布洛克县 -1111.1

- user3648073

1

也许在原帖中粘贴一个 dput(obesity) 的输出结果 :-) - hrbrmstr

@hrbrmstr，你在吗？我有一些关于ggcounty的问题要问你。顺便说一下，非常感谢你的帮助！ - user3648073

显示剩余2条评论

6个回答

16

这个例子与你的 obesity_map 数据集格式类似，但是尽可能适应它的格式。它还使用数据表联接，比 merge(...) 更快，特别是在像你的数据集这样大的情况下。

library(ggplot2)
# this creates an example formatted as your obesity.map - you have this already...
set.seed(1)    # for reproducible example
map.county <- map_data('county')
counties   <- unique(map.county[,5:6])
obesity_map <- data.frame(state_names=counties$region, 
                          county_names=counties$subregion, 
                          obesity= runif(nrow(counties), min=0, max=100))

# you start here...
library(data.table)   # use data table merge - it's *much* faster
map.county <- data.table(map_data('county'))
setkey(map.county,region,subregion)
obesity_map <- data.table(obesity_map)
setkey(obesity_map,state_names,county_names)
map.df      <- map.county[obesity_map]

ggplot(map.df, aes(x=long, y=lat, group=group, fill=obesity)) + 
  geom_polygon()+coord_map()

此外，如果你的数据集中有FIPS代码（看起来确实是这样的），我强烈建议您使用美国人口普查局的TIGER / Line县级地图文件（也包含这些代码）并将其合并在一起。这更加可靠。例如，在您提取的obesity_map数据框中，州和县名的首字母大写，而在R内置的县级数据集中，它们不是，因此您需要处理这个问题。此外，TIGER文件是最新的，而内部数据集则不是。

所以这是一个很有趣的问题。事实证明，实际的肥胖数据可以在USDA网站上下载，以MSExcel格式的文件形式出现，您可以从这里下载。人口普查局网站上还有美国县的shapfile，在这里。两个文件都具有FIPS信息。在R中，可以相对简单地将它们组合在一起：

library(XLConnect)    # for loadWorkbook(...) and readWorksheet(...)
library(rgdal)        # for readOGR(...)
library(RcolorBrewer) # for brewer.pal(...)
library(data.table)

setwd(" < directory with all your files > ")
wb <- loadWorkbook("DataDownload.xls")   # from the USDA website
df <- readWorksheet(wb,"HEALTH")         # this sheet has the obesity data

US.counties <- readOGR(dsn=".",layer="gz_2010_us_050_00_5m")
#leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US.counties <- US.counties[!(US.counties$STATE %in% c("02","15","72")),]  
county.data <- US.counties@data
county.data <- cbind(id=rownames(county.data),county.data)
county.data <- data.table(county.data)
county.data[,FIPS:=paste0(STATE,COUNTY)] # this is the state + county FIPS code
setkey(county.data,FIPS)      
obesity.data <- data.table(df)
setkey(obesity.data,FIPS)
county.data[obesity.data,obesity:=PCT_OBESE_ADULTS10]

map.df <- data.table(fortify(US.counties))
setkey(map.df,id)
setkey(county.data,id)
map.df[county.data,obesity:=obesity]

ggplot(map.df, aes(x=long, y=lat, group=group, fill=obesity)) +
  scale_fill_gradientn("",colours=brewer.pal(9,"YlOrRd"))+
  geom_polygon()+coord_map()+
  labs(title="2010 Adult Obesity by Country, percent",x="",y="")+
  theme_bw()

要生成这个：

- jlhoward

哇，谢谢你啊！我很感激！这真的解决了我的问题。最后一个问题。你知道我怎么把值为0的单独设置成一个颜色吗？我把所有的NA值都设置成了0，我想让它显示为灰色。再次感谢你的帮助！ - user3648073

实际上，如果你将它们保留为NA，它们应该默认显示为灰色。由于您是SO的新手，请阅读此内容。 - jlhoward

非常感谢您提供这个很好的使用data.table的示例。我仍在努力理解您在第一个map中如何使用map.county[obesity_map]，特别是它是用于合并吗？感谢您的分享。 - Paulo E. Cardoso

@jlhoward 这段代码的第二部分使用 FIPS 在脚本中能够运行，但在函数内（即使是一个包装器）却失败了，并出现了这个错误：

Error in :=(FIPS, paste0(STATE, COUNTY)) :  Check that is.data.table(DT) == TRUE. Otherwise, := and :=(...) are defined for use in j, once only and in particular ways. See help(":=").

您有任何想法吗？ - miguelmorin

8

我可以通过管理映射变量来使其正常工作。将其重命名为“区域”。

library(ggplot2)
library(maps)
m.usa <- map_data("county")
m.usa$id <- m.usa$subregion
m.usa <- m.usa[ ,-5]
names(m.usa)[5] <- 'region'


df <- data.frame(region = unique(m.usa$region),
                 obesity = rnorm(length(unique(m.usa$region)), 50, 10),
                 stringsAsFactors = F)

head(df)
region  obesity
1 autauga 44.54833
2 baldwin 68.61470
3 barbour 52.19718
4    bibb 50.88948
5  blount 42.73134
6 bullock 59.93515

ggplot(df, aes(map_id = region)) +
  geom_map(aes(fill = obesity), map = m.usa) + 
  expand_limits(x = m.usa$long, y = m.usa$lat) +
  coord_map()

geom_map

- Paulo E. Cardoso

非常感谢！我真的很感激你的帮助！ - user3648073

这个可行！问题在于 map_data(county) 使用县和州名称作为标识符，对于拼写不同的县（例如路易斯安那州的“Ste Genevieve”/“Ste Genevieve”）会造成问题。这篇 R-bloggers 帖子建议使用由人口普查局提供的 FIPS 代码形状文件。 - miguelmorin

1

在@jlhoward的回答基础上，使用data.table的代码以一种神秘的方式失败了：

 Error in `:=`(FIPS, paste0(STATE, COUNTY)) : 
  Check that is.data.table(DT) == TRUE. Otherwise, := and `:=`(...) are defined for use in j, once only and in particular ways. See help(":=").

这个错误发生过几次，但仅当代码位于函数内部时才会出现，即使只是一个最小的包装器也是如此。在脚本中运行正常。尽管现在我无法再现这个错误，但为了完整起见，我使用了merge()而不是data.table来调整代码：

library(rgdal)        # for readOGR(...)
library(ggplot2)      # for fortify() and plot()
library(RColorBrewer) # for brewer.pal(...)

US.counties <- readOGR(dsn=".",layer="gz_2010_us_050_00_5m")
#leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US.counties <- US.counties[!(US.counties$STATE %in% c("02","15","72")),]
county.data <- US.counties@data

county.data <- cbind(id=rownames(county.data),county.data)
county.data$FIPS <- paste0(county.data$STATE, county.data$COUNTY) # this is the state + county FIPS code

df <- data.frame(FIPS=county.data$FIPS,
                 PCT_OBESE_ADULTS10= runif(nrow(county.data), min=0, max=100))

# Merge county.data to obesity
county.data <- merge(county.data,
                     df,
                     by.x = "FIPS",
                     by.y = "FIPS")

map.df <- fortify(US.counties)

# Merge the map to county.data
map.df <- merge(map.df, county.data, by.x = "id", by.y = "id")

ggplot(map.df, aes(x=long, y=lat, group=group, fill=PCT_OBESE_ADULTS10)) +
  scale_fill_gradientn("",colours=brewer.pal(9,"YlOrRd"))+
  geom_polygon()+coord_map()+
  labs(title="2010 Adult Obesity by Country, percent",x="",y="")+
  theme_bw()

- miguelmorin

1

我认为你所需要做的就是像之前处理 map.data 变量一样重新排列 map.county 变量。

....
map.county <- merge(county.obesity, map.county, all=TRUE)

## reorder the map before plotting
map.county <- map.county[order(map.data$county),] 

## plot
ggplot(map.county, aes(x = long, y = lat, group=group, fill=as.factor(value))) + geom_polygon(colour = "white", size = 0.1)

- Scott Worland

0

我在使用TMAP和空间数据方面还有点新手，但是想跟进Martijn Tennekes的帖子。按照他的建议，我在第二张地图（带州边界）中遇到了一个错误。当运行以下代码时：

US_state <- unionSpatialPolygons(US,US$STATE)

我一直遇到这个错误：“Error in unionSpatialPolygons(US, US$STATE) : not a SpatialPolygons object”。

为了纠正这个问题，我不得不使用另一个变量，并将其作为一个空间多边形数据框运行：

US <- read_shape("gz_2010_us_050_00_20m.shp")
US2<-readShapeSpatial("gz_2010_us_050_00_20m.shp")

US <- US[!(US$STATE %in% c("02","15","72")),]  

US$FIPS <- paste0(US$STATE, US$COUNTY)
US <- append_data(US, med_inc_df, key.shp = "FIPS", key.data = "GEOID")

#the difference is here:
US_states <- unionSpatialPolygons(US2, US2$STATE)

tm_shape(US, projection="+init=epsg:2163") +
  tm_polygons("estimate", border.col = "grey30", title="") +
  tm_shape(US_states) +
  tm_borders(lwd=2, col = "black", alpha = .5) +
  tm_layout(title="2016 Median Income by County", 
            title.position = c("center", "top"),
            legend.text.size=1)

我的地图

- Andy White

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Martijn Tennekes · Accepted Answer

也许回答有点晚了，但我认为仍然值得分享。

数据的读取和预处理与jlhoward的回答类似，但有一些区别：

library(tmap)      # package for plotting
library(readxl)    # for reading Excel
library(maptools)  # for unionSpatialPolygons

# download data
download.file("http://www.ers.usda.gov/datafiles/Food_Environment_Atlas/Data_Access_and_Documentation_Downloads/Current_Version/DataDownload.xls", destfile = "DataDownload.xls", mode="wb")
df <- read_excel("DataDownload.xls", sheet = "HEALTH")

# download shape (a little less detail than in the other scripts)
f <- tempfile()
download.file("http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip", destfile = f)
unzip(f, exdir = ".")
US <- read_shape("gz_2010_us_050_00_20m.shp")

# leave out AK, HI, and PR (state FIPS: 02, 15, and 72)
US <- US[!(US$STATE %in% c("02","15","72")),]  

# append data to shape
US$FIPS <- paste0(US$STATE, US$COUNTY)
US <- append_data(US, df, key.shp = "FIPS", key.data = "FIPS")

当正确的数据附加到形状对象时，可以通过一行代码绘制色彩分段图：

qtm(US, fill = "PCT_OBESE_ADULTS10")

通过添加州边界、更好的地图投影和标题，可以对此进行增强：

# create shape object with state polygons
US_states <- unionSpatialPolygons(US, IDs=US$STATE)

tm_shape(US, projection="+init=epsg:2163") +
  tm_polygons("PCT_OBESE_ADULTS10", border.col = "grey30", title="") +
tm_shape(US_states) +
  tm_borders(lwd=2, col = "black", alpha = .5) +
tm_layout(title="2010 Adult Obesity by County, percent", 
          title.position = c("center", "top"),
          legend.text.size=1)