使用美国县级数据创建 Choropleth 地图

Posted

技术标签:

【中文标题】使用美国县级数据创建 Choropleth 地图【英文标题】:Creating a Choropleth map with US county level data 【发布时间】:2020-07-05 12:49:20 【问题描述】:

我正在尝试使用 R 生成有关 COVID-19 感染的县级数据的等值线图。我是 R 的相对新手,所以......

我已经用 ggmap 做了一些相当基本的东西来绘制空间数据,但从来没有像这样的东西。通常我只有需要在地图上叠加的兴趣点,所以我可以使用 geom_point 及其纬度/经度。在这种情况下,我需要构建底层地图,然后填充区域,而我在 ggplot 世界中正在努力做到这一点。

我已经遵循了一些我发现的在线示例:

library(ggplot2)
library(broom)
library(geojsonio)

#get a county level map geoJSON file
counties <- geojson_read("https://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_050_00_500k.json", what = "sp")

#filter our alaska and Hawaii
lower48 <- counties[(counties@data$STATE != "02" & counties@data$STATE != "15") ,]

#turn it into a dataframe for ggmap
new_counties <- tidy(lower48)

# Plot it
print(ggplot() +
  geom_polygon(data = new_counties, aes( x = long, y = lat, group = group), fill="#69b3a2", color="white") +
  theme_void() +
  coord_map())

这会产生这个情节:

到目前为止一切顺利。但我的 new_counties 数据框现在看起来像这样:

head(new_counties)
# A tibble: 6 x 7
   long   lat order hole  piece group id  
  <dbl> <dbl> <int> <lgl> <chr> <chr> <chr>
1 -85.4  33.9     1 FALSE 1     0.1   0    
2 -85.4  33.9     2 FALSE 1     0.1   0    
3 -85.4  33.9     3 FALSE 1     0.1   0    
4 -85.4  33.9     4 FALSE 1     0.1   0    
5 -85.4  33.9     5 FALSE 1     0.1   0    
6 -85.4  33.8     6 FALSE 1     0.1   0 

因此,我丢失了任何可以与我县级感染数据相关联的信息。

我的数据中每个县都有一个 5 位数的 FIPS 代码。前两位数字是州,后三位是县。我的 geoJSON 文件有更详细的 FIPS 代码。我尝试只抓取前 5 个并创建我自己的数据元素,我可以映射回

library(ggplot2)
library(broom)
library(geojsonio)

#get a county level map geoJSON file
counties <- geojson_read("https://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_050_00_500k.json", what = "sp")

#filter our alaska and Hawaii
lower48 <- counties[(counties@data$STATE != "02" & counties@data$STATE != "15") ,]

#add my own FIPS code
lower48@data$myFIPS <- substr(as.character(lower48@data$GEO_ID),1,5)  

#turn it into a dataframe for ggmap
new_counties <- tidy(lower48, region = "myFIPS")


# Plot it
print(ggplot() +
  geom_polygon(data = new_counties, aes( x = long, y = lat, group = group), fill="#69b3a2", color="white") +
  theme_void() +
  coord_map())

但这会产生这个情节

我不得不说我对 broom::tidy 不够熟悉,无法确切知道原因。 我还注意到,当我键入此内容时,我需要过滤掉波多黎各!

如果有人能指出我有用的方向....我不拘泥于当前的方法,但我想坚持使用 ggplot2 或 ggmap。我的老板最终希望我覆盖某些功能。最终的目标是效仿here 的示例,并生成一个显示随时间变化的数据的动画地图,但我显然还有很长的路要走。

【问题讨论】:

使用sf::st_read() 读取geojson,因此您将其作为sf 对象获取。然后使用ggplot2::geom_sf() 绘制它。 - sfsp 的继承者,所以我建议任何人离开 sp 其中sf 对象已经是一个data.frame。因此,所有标准过滤和子集操作都“正常工作”。 【参考方案1】:

有很多方法可以做到这一点,但概念基本相同:查找包含*** FIPS 代码的地图并使用它们与数据源链接,也包含相同的 FIPS 代码以及用于绘图的变量 (这里是每天的 covid-19 病例数)。

#devtools::install_github("UrbanInstitute/urbnmapr")
library(urbnmapr) # For map
library(ggplot2)  # For map
library(dplyr)    # For summarizing
library(tidyr)    # For reshaping
library(stringr)  # For padding leading zeros

# Get COVID cases, available from:
url <- "https://static.usafacts.org/public/data/covid-19/covid_confirmed_usafacts.csv
             ?_ga=2.162130428.136323622.1585096338-408005114.1585096338"

COV <- read.csv(url, stringsAsFactors = FALSE)
names(COV)[1] <- "countyFIPS"  # Fix the name of first column. Why!?

数据以宽格式存储,每个县的每日病例分布在列中。这需要在与地图合并之前收集。日期需要转换为正确的日期。 FIPS 代码存储为整数,因此需要将这些代码转换为带有前导 0 的字符,以便与地图数据合并。我为地图数据使用 urbnmap 包。

Covid <- pivot_longer(COV, cols=starts_with("X"), 
                     values_to="cases",
                     names_to=c("X","date_infected"),
                     names_sep="X") %>%                
  mutate(date_infected = as.Date(date_infected, format="%m.%d.%Y"),
         countyFIPS = str_pad(as.character(countyFIPS), 5, pad="0"))

# Obtain map data for counties (to link with covid data) and states (for showing borders)
states_sf <- get_urbn_map(map = "states", sf = TRUE)
counties_sf <- get_urbn_map(map = "counties", sf = TRUE)

# Merge county map with total cases of cov
counties_cov <- inner_join(counties_sf, group_by(Covid, countyFIPS) %>%
       summarise(cases=sum(cases)), by=c("county_fips"="countyFIPS"))

counties_cov %>%
  ggplot() +
  geom_sf(mapping = aes(fill = cases), color = NA) +
  geom_sf(data = states_sf, fill = NA, color = "black", size = 0.25) +
  coord_sf(datum = NA) +   
  scale_fill_gradient(name = "Cases", trans = "log", low='pink', high='navyblue', 
                      na.value="white", breaks=c(1, max(counties_cov$cases))) +
  theme_bw() + theme(legend.position="bottom", panel.border = element_blank())


对于动画,您可以使用 gganimate 包和过渡。除了不应汇总 covid 数据外,这些命令与上述类似。

library(gganimate)

counties_cov <- inner_join(counties_sf, Covid, by=c("county_fips"="countyFIPS"))

p <- ggplot(counties_cov) + ... # as above

p <- p + transition_time(date_infected) +
  labs(title = 'Date: frame_time')

animate(p, end_pause=30)

【讨论】:

太棒了...我使用choroplethr 包取得了一些进展,但这让我感到羞耻。我会研究这个... 由于某种原因,它不喜欢它从 URL 中读取的countyFIPS 列名,并将 i 元音变音(或任何你称之为的)粘贴在那里,所以我手动删除了它。你知道在地图上覆盖一些纬度/经度 geom_points 需要什么转换吗?照原样,他们都在南达科他州结束时发出警告Transformation introduced infinite values in discrete y-axis。我尝试使用另一个包时遇到了这个问题,它与地图使用的投影有关。该软件包提供了一个 transform() 函数来使纬度/经度的可绘制... 忽略...使用 usmap_transform 包中的 usmap_transform 函数解决了这个问题 我修复了关于 date_infected 的小错误。对于灰色县,在整个期间至少有一天有病例,而白色县在任何时候都没有病例。我认为这是一个不错的“功能”,但应该有一种方法可以禁用它。 ... 如果不喜欢灰色,可以在 scale_fill_gradient 函数中添加 na.values="white" 或 na.values="transparent"。

以上是关于使用美国县级数据创建 Choropleth 地图的主要内容,如果未能解决你的问题,请参考以下文章

对于美国县级地图,绘图可视化太慢

choropleth地图数据未加载到由bing提供的excel地图中的状态IN和MN

Plotly:如何在 Choropleth 地图上显示州线

D3 choropleth 状态图数据更新按钮单击

如何在 Plotly choropleth 地图上创建符号/按钮

使用 Altair 显示英国的 Choropleth 地图时出现问题