如何在 R 的空间叠加中包含夏威夷和阿拉斯加
How to include Hawaii and Alaska in spatial overlay in R
我在获取人口统计数据点以覆盖在美国县地图上时遇到了一些麻烦。我可以很好地绘制地图,但没有显示夏威夷和阿拉斯加的数据。我已经确定了问题的根源 - 它是在我的 over
命令之后。我的工作流程使用一个 csv 文件,可以在此处找到 (https://www.dropbox.com/s/0arazi2n0adivzc/data.dem2.csv?dl=0)。这是我的工作流程:
#Load dependencies
devtools::install_github("hrbrmstr/albersusa")
library(albersusa)
library(dplyr)
library(rgeos)
library(maptools)
library(ggplot2)
library(ggalt)
library(ggthemes)
library(viridis)
#Read Data
df<-read.csv("data.dem.csv")
#Retreive polygon shapefile
counties_composite() %>%
subset(df$state %in% unique(df$state)) -> usa #Note I've checked here and Alaska is present, see below
#Subset just points and create spatial points object
pts <- df[,4:1]
pts<-as.data.frame(pts)
coordinates(pts) <- ~long+lat
proj4string(pts) <- CRS(proj4string(usa)) #Note I've checked here as welland Alaska is present still, see here
#Spatial overlay
b<-over(pts, usa) #This is where the problem arises: see here
b<-select(b, -state)
b<-bind_cols(df, b)
bind_cols(df, select(over(pts, usa), -state)) %>%
count(fips, wt=count) -> df
usa_map <- fortify(usa, region="tips")
ggplot()+
geom_map(data=usa_map, map=usa_map,
aes(long, lat, map_id=id),
color="#b2b2b2", size=0.05, fill="grey") +
geom_map(data=df, map=usa_map,
aes(fill=n, map_id=fips),
color="#b2b2b2", size=0.05) +
scale_fill_viridis(name="Count", trans="log10") +
gg + coord_map() +
theme_map() +
theme(legend.position=c(0.85, 0.2))
您可能会怀疑,最终输出没有显示阿拉斯加或夏威夷的数据。我不确定发生了什么,但似乎 sp 包中的 over
命令是问题的根源。非常感谢任何建议。
请注意,这是一个与找到的问题不同的问题 Relocating Alaska and Hawaii on thematic map of the USA with ggplot2 and How do you create a 50 state map (instead of just lower-48)
题目相互之间没有任何关系。这不是重复的。第一个问题是关于夏威夷和阿拉斯加实际多边形的位置,正如您从我的地图上看到的那样,我没有这个问题。第二个 link 是关于获取包括夏威夷和阿拉斯加的地图。同样,我的地图包括这两个,但在我的数据处理工作流程的某个地方,这两个的数据被删除(特别是叠加功能)。请不要标记为重复。
你需要做的比之前的答案多一点,因为复合 shapefile - 根据它的定义 - 将阿拉斯加和夏威夷从它们的原始位置移动,这将使 over()
在尝试匹配点时错过它们多边形。这很容易解决:
library(albersusa) # devtools::install_github("hrbrmstr/albersusa)
library(readr)
library(dplyr)
library(rgeos)
library(rgdal)
library(maptools)
library(ggplot2)
library(ggalt)
library(ggthemes)
library(viridis)
df <- read_csv("data.dem2.csv")
# need this for the composite map & no need to subset
usa <- counties_composite()
# need this for the "over" since the composite map totally
# messes with the lon/lat positions of alaska & hawaii
URL <- "http://eric.clst.org/wupl/Stuff/gz_2010_us_050_00_500k.json"
fil <- basename(URL)
if (!file.exists(fil)) download.file(URL, fil)
orig_counties <- readOGR(fil, "OGRGeoJSON", stringsAsFactors=FALSE)
# your new csv has an extra column at the beginning
pts <- as.data.frame(df[,3:2])
coordinates(pts) <- ~long+lat
proj4string(pts) <- CRS(proj4string(orig_counties))
# don't need to select out the duplicate col name anymore
# but we do need to create the FIPS code
bind_cols(df, over(pts, orig_counties)) %>%
mutate(fips=sprintf("%s%s", STATE, COUNTY)) %>%
count(fips, wt=count) -> df
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
我在获取人口统计数据点以覆盖在美国县地图上时遇到了一些麻烦。我可以很好地绘制地图,但没有显示夏威夷和阿拉斯加的数据。我已经确定了问题的根源 - 它是在我的 over
命令之后。我的工作流程使用一个 csv 文件,可以在此处找到 (https://www.dropbox.com/s/0arazi2n0adivzc/data.dem2.csv?dl=0)。这是我的工作流程:
#Load dependencies
devtools::install_github("hrbrmstr/albersusa")
library(albersusa)
library(dplyr)
library(rgeos)
library(maptools)
library(ggplot2)
library(ggalt)
library(ggthemes)
library(viridis)
#Read Data
df<-read.csv("data.dem.csv")
#Retreive polygon shapefile
counties_composite() %>%
subset(df$state %in% unique(df$state)) -> usa #Note I've checked here and Alaska is present, see below
#Subset just points and create spatial points object
pts <- df[,4:1]
pts<-as.data.frame(pts)
coordinates(pts) <- ~long+lat
proj4string(pts) <- CRS(proj4string(usa)) #Note I've checked here as welland Alaska is present still, see here
#Spatial overlay
b<-over(pts, usa) #This is where the problem arises: see here
b<-select(b, -state)
b<-bind_cols(df, b)
bind_cols(df, select(over(pts, usa), -state)) %>%
count(fips, wt=count) -> df
usa_map <- fortify(usa, region="tips")
ggplot()+
geom_map(data=usa_map, map=usa_map,
aes(long, lat, map_id=id),
color="#b2b2b2", size=0.05, fill="grey") +
geom_map(data=df, map=usa_map,
aes(fill=n, map_id=fips),
color="#b2b2b2", size=0.05) +
scale_fill_viridis(name="Count", trans="log10") +
gg + coord_map() +
theme_map() +
theme(legend.position=c(0.85, 0.2))
您可能会怀疑,最终输出没有显示阿拉斯加或夏威夷的数据。我不确定发生了什么,但似乎 sp 包中的 over
命令是问题的根源。非常感谢任何建议。
请注意,这是一个与找到的问题不同的问题 Relocating Alaska and Hawaii on thematic map of the USA with ggplot2 and How do you create a 50 state map (instead of just lower-48)
题目相互之间没有任何关系。这不是重复的。第一个问题是关于夏威夷和阿拉斯加实际多边形的位置,正如您从我的地图上看到的那样,我没有这个问题。第二个 link 是关于获取包括夏威夷和阿拉斯加的地图。同样,我的地图包括这两个,但在我的数据处理工作流程的某个地方,这两个的数据被删除(特别是叠加功能)。请不要标记为重复。
你需要做的比之前的答案多一点,因为复合 shapefile - 根据它的定义 - 将阿拉斯加和夏威夷从它们的原始位置移动,这将使 over()
在尝试匹配点时错过它们多边形。这很容易解决:
library(albersusa) # devtools::install_github("hrbrmstr/albersusa)
library(readr)
library(dplyr)
library(rgeos)
library(rgdal)
library(maptools)
library(ggplot2)
library(ggalt)
library(ggthemes)
library(viridis)
df <- read_csv("data.dem2.csv")
# need this for the composite map & no need to subset
usa <- counties_composite()
# need this for the "over" since the composite map totally
# messes with the lon/lat positions of alaska & hawaii
URL <- "http://eric.clst.org/wupl/Stuff/gz_2010_us_050_00_500k.json"
fil <- basename(URL)
if (!file.exists(fil)) download.file(URL, fil)
orig_counties <- readOGR(fil, "OGRGeoJSON", stringsAsFactors=FALSE)
# your new csv has an extra column at the beginning
pts <- as.data.frame(df[,3:2])
coordinates(pts) <- ~long+lat
proj4string(pts) <- CRS(proj4string(orig_counties))
# don't need to select out the duplicate col name anymore
# but we do need to create the FIPS code
bind_cols(df, over(pts, orig_counties)) %>%
mutate(fips=sprintf("%s%s", STATE, COUNTY)) %>%
count(fips, wt=count) -> df
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