我正在尝试创建一个显示美国所有50个州的专题 map ,但是我无法以可靠的方式重新安置阿拉斯加和夏威夷。我有几个主意,但都没有一个很好的方法。我现在将演示它们。
首先,我们需要导入数据;仅使用maps
包中的数据是不够的,因为它不包括Hawaii和Alaska。
setwd(tempdir())
download.file("https://dl.dropbox.com/s/wl0z5rpygtowqbf/states_21basic.zip?dl=1",
"usmapdata.zip",
method = "curl")
# This is a mirror of http://www.arcgis.com/home/item.html?
# id=f7f805eb65eb4ab787a0a3e1116ca7e5
unzip("usmapdata.zip")
require(rgdal)
all_states <- readOGR("states_21basic/", "states")
require(ggplot2); require(maptools); require(rgeos); require(mapproj);
all_states <- fortify(all_states, region = "STATE_NAME")
现在,我们定义一些情节美学:p <- ggplot() + geom_polygon(
aes(x=long, y=lat, group = group, fill = as.numeric(as.factor(id))),
colour="white", size = 0.25
) + coord_map(projection="azequalarea") +
scale_fill_gradient(limits = c(1,50))
现在我们删除所有背景等,以便在我们重叠非连续状态时它们不会发生冲突:p <- p + theme(axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.background=element_blank(),
panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
plot.background=element_blank())
使用视口(viewport)我的第一个想法是使用视口(viewport):
AK <- p %+% subset(all_states, id == "Alaska") + theme(legend.position = "none")
HI <- p %+% subset(all_states, id == "Hawaii") + theme(legend.position = "none")
contiguous <- p %+% subset(all_states, id != "Alaska" & id != "Hawaii")
grid.newpage()
vp <- viewport(width = 1, height = 1)
print(contiguous, vp = vp)
subvp1 <- viewport(width = 0.25, height = 0.25, x = 0.18, y = 0.33)
print(AK, vp = subvp1)
subvp2 <- viewport(width = 0.12, height = 0.12, x = 0.32, y = 0.27)
print(HI, vp = subvp2)
这看起来不错,但并不令人满意,因为它对图形的细微变化非常敏感,例如调整大小或更改图例的大小和形状。
手动移动阿拉斯加和夏威夷
all_states_AKHImoved <- within(all_states, {
lat[id == "Alaska"] <- lat[id == "Alaska"] - 45
long[id == "Alaska"] <- long[id == "Alaska"] + 40
lat[id == "Hawaii"] <- lat[id == "Hawaii"] + 0
long[id == "Hawaii"] <- long[id == "Hawaii"] + 70
})
p %+% all_states_AKHImoved
这并不令人满意,因为阿拉斯加通常无法在大多数美国 map 上缩放,因此看起来很大。同样,将阿拉斯加和夏威夷搬迁也会改变 map 投影所引起的变形。
题
有谁有更好的方法?
最佳答案
这是通过投影和变换来实现的方法。你会需要:
require(maptools)
require(rgdal)
fixup <- function(usa,alaskaFix,hawaiiFix){
alaska=usa[usa$STATE_NAME=="Alaska",]
alaska = fix1(alaska,alaskaFix)
proj4string(alaska) <- proj4string(usa)
hawaii = usa[usa$STATE_NAME=="Hawaii",]
hawaii = fix1(hawaii,hawaiiFix)
proj4string(hawaii) <- proj4string(usa)
usa = usa[! usa$STATE_NAME %in% c("Alaska","Hawaii"),]
usa = rbind(usa,alaska,hawaii)
return(usa)
}
fix1 <- function(object,params){
r=params[1];scale=params[2];shift=params[3:4]
object = elide(object,rotate=r)
size = max(apply(bbox(object),1,diff))/scale
object = elide(object,scale=size)
object = elide(object,shift=shift)
object
}
然后读取您的shapefile。使用
rgdal
:us = readOGR(dsn = "states_21basic",layer="states")
现在转换为等面积,并运行fixup函数:
usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix)
这些参数分别是阿拉斯加和夏威夷的旋转,缩放,x和y位移,并通过反复试验获得。仔细调整它们。由于涉及大量数字,即使将Hawaii的scale参数更改为0.99999,也将其从地球上发送出去。
如果您想将其转回纬度:
usfixLL = spTransform(usfix,CRS("+init=epsg:4326"))
plot(usfixLL)
但是我不确定您是否需要使用
ggplot
中的转换,因为我们已经使用spTransform
完成了转换。现在,您可以跳过
ggplot2
加强业务。我不确定这是否对您重要,但是请注意usfix
版本中的州顺序有所不同-阿拉斯加和夏威夷现在是最后两个州。关于r - 使用ggplot2在美国专题 map 上重新定位阿拉斯加和夏威夷,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/13757771/