无法在r中构建区域分布地图

Question

无法在r中构建区域分布地图

3

我有一些人口统计数据，想用它来制作美国县份的等值线地图。我的工作流程没有出现任何错误，我能够创建最终地图，但是它所映射的数据是不正确的。我的工作流程使用了两个数据源——一个形状文件和一个数据框架。这个形状文件是一个县份形状文件，可以在这个链接https://www.dropbox.com/s/4ujxidyx42793j7/cb_2015_us_county_500k.zip?dl=1找到。数据框架文件可以在这个链接https://www.dropbox.com/s/qys6s6ikrs1g2xb/data.dem.csv?dl=1找到。

以下是我的代码：

#Load dependencies
library(sp)
library(spatialEco)
library(rgdal)
library(dplyr)
library(maptools)
library(taRifx.geo)
library(ggplot2)
library(USAboundaries)
library(splitstackshape)
library(maps)
library(cowplot)

#Read in shape and csv files
county.track<-readOGR("/path", "filename")
county.track@data$id = rownames(county.track@data)
data<-read.csv("/path/filename.csv")

#Convert data.frame (data) to points polygon file
data$y<-data$lat
data$x<-data$long
coordinates(data) <- ~ x + y
proj4string(data) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(county.track) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

#Overlay points onto polygons
county.track.data<-point.in.poly(data, county.track)

#Summarize point data by county
count<-select(as.data.frame(county.track.data), id, count)
count<-count %>%
  group_by(id) %>%
  summarize(count=sum(count))

#Merge with shape file data
county.track@data<-merge(county.track@data, count, by="id", all.x=T)

#Replace NA values with zeroes 
county.track@data$count[is.na(county.track@data$count)]<-0
county.track.points = fortify(county.track, region="id")
map.plot<-merge(county.track.points, county.track@data, by="id")

#Get rid of Hawaii and Alaska
map.plot<-map.plot %>%
  filter(lat<50 & lat>25) %>%
  filter(long>-130)

#Create choropleth map using ggplot2
 ggplot(map.plot) +
  geom_polygon(aes(long, lat, group=group, fill=log(count))) +
  coord_map()

输出结果应如下所示：

但这是错误的，有很多原因可以看出来。首先，显而易见的是大部分数据没有被映射。地图上的灰色区域表示NA。但我在上面的步骤中删除了NAs，而且在检查用于映射（map.plot）的数据时，填充变量（计数）中没有NAs。其次，映射值的分布不正确。洛杉矶县的计数值应该是793（对数值为6.675823），然而在地图上，许多颜色较浅的县表明其他空间单位的值更高，而一些排名前几位的县，如圣迭戈，根本没有被填充（地图的左下角）。

当我检查用于绘制地图（map.plot）的数据时，一切似乎都很好。洛杉矶县仍然是“count”变量最高价值的县，但地图表明情况并非如此（请参见此处的图像）。

我希望有人能在这里进行一些取证，并确定问题所在。我已经尽力检查了所有步骤，但似乎无法确定问题。提前致谢。

更新：我尝试使用来自同一来源的不同shapefile。上面链接中的shapefile与以下网址（https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html）标记为“cb_2015_us_county_500k.zip”的文件相同。当我选择不同的shapefile（例如cb_2015_us_county_5m.zip）时，我得到了一个不同的地图，但是出现了相同的问题：请参见以下地图作为例子：

我不确定发生了什么事情！在这张新地图中，洛杉矶县甚至没有被着色，但橙县却被着色了！非常感谢任何帮助。

- Cyrus Mohammadian

1个回答

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

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

- hrbrmstr · Accepted Answer

我不太确定你的合并操作是怎样的，但这个方法对我很有效：

library(albersusa) # devtools::install_github("hrbrmstr/albersusa)
library(readr)
library(dplyr)
library(rgeos)
library(maptools)
library(ggplot2)
library(ggalt)
library(ggthemes)
library(viridis)

df <- read_csv("data.dem.csv")

counties_composite() %>% 
  subset(state %in% unique(df$state)) -> usa

pts <- df[,2:1]
coordinates(pts) <- ~long+lat
proj4string(pts) <- CRS(proj4string(usa))

bind_cols(df, select(over(pts, usa), -state)) %>% 
  count(fips, wt=count) -> df

您拥有942个县：

glimpse(df)
## Observations: 942
## Variables: 2
## $ fips <chr> "01001", "01003", "01013", "01015", "01043", "01055", "01061", ...
## $ n    <int> 1, 2, 1, 3, 1, 3, 1, 1, 19, 6, 12, 7, 7, 1, 4, 4, 1, 5, 67, 19,...

美国有超过3千个县。

然而，并不是很多县拥有NA。

filter(df, is.na(fips))
## # A tibble: 1 x 2
##    fips     n
#3   <chr> <int>
## 1  <NA>    10

usa_map <- fortify(usa, region="fips")

gg <- ggplot()
gg <- gg + geom_map(data=usa_map, map=usa_map,
                    aes(long, lat, map_id=id),
                    color="#b2b2b2", size=0.05, fill="white")
gg <- gg + geom_map(data=df, map=usa_map,
                    aes(fill=n, map_id=fips),
                    color="#b2b2b2", size=0.05)
gg <- gg + scale_fill_viridis(name="Count", trans="log10")
gg <- gg + coord_proj(us_aeqd_proj)
gg <- gg + theme_map()
gg <- gg + theme(legend.position=c(0.85, 0.2))
gg