27

我正在尝试创建一张显示美国所有 50 个州的专题地图,但我无法以可靠的方式重新定位阿拉斯加和夏威夷。我有几个想法,但没有一个能很好地工作。我现在将展示它们。

首先我们需要导入数据;使用maps包中的数据是不够的,因为它不包括夏威夷和阿拉斯加。

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())

使用视口

我的第一个想法是使用视口:

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

在此处输入图像描述

这并不令人满意,因为阿拉斯加在大多数美国地图上通常不会按比例缩放,因此它看起来非常大。重新定位阿拉斯加和夏威夷也改变了地图投影引入的失真。

问题

有没有人有更好的方法?

4

3 回答 3

14

以下是如何通过投影和转换来做到这一点。你会需要:

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 位移,是通过反复试验获得的。仔细调整它们。甚至将夏威夷的比例参数更改为 0.99999 也将其从地球上发射出去,因为涉及的人数众多。

如果您想将其转回经纬度:

usfixLL = spTransform(usfix,CRS("+init=epsg:4326"))
plot(usfixLL)

但我不确定您是否需要使用转换,ggplot因为我们已经使用spTransform.

您现在可以跳过ggplot2强化业务。我不确定这对您是否重要,但请注意usfix版本中的州顺序不同 - 阿拉斯加和夏威夷现在是最后两个州。

于 2012-12-07T17:23:17.567 回答
14

fiftystater在 GitHub ( devtools::install_github("wmurphyrd/fiftystater")) 上发布了 R 包以提供一个简单的解决方案。它基于 Spacedman 的答案中的步骤(将链接但代表不足),并作为一个ggplot2::geom_map现成的形状数据框发布,命名fifty_states为消除安装依赖项、追踪源形状文件或调整省略值的需要。

library(ggplot2)
library(mapproj)
library(fiftystater)

crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests)

p <- ggplot(crimes, aes(map_id = state)) + 
  geom_map(aes(fill = Assault), map = fifty_states) + 
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  coord_map()
p

五十州地图 五十州地图

可以按通常的方式清理绘图噪声,并且fifty_states_inset_boxes包中还有添加插入框的功能:

p + scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  labs(x = "", y = "") +
  theme(panel.background = element_blank()) +
  fifty_states_inset_boxes()

五十个州的清洁工,带插盒 五十个带有插图框的州

于 2017-07-06T02:49:18.437 回答
3

一旦你开始像这样改变事物,你不妨将阿拉斯加和夏威夷表示为墨西哥湾某处的方盒。它还有一个额外的好处,那就是可以分辨夏威夷是什么颜色。

然后你可能会一直走下去,使用一个扭曲的系统,每个州都有相等的面积,然后你就可以看到罗德岛了。

谷歌图像上的美国制图示例显示了这种情况。不知道他们中有多少人有 shapefile 或数据。

你真的对状态的相对大小感兴趣,或者你想要一个让人们看到一个状态的价值的表示吗?

于 2012-12-07T08:32:27.303 回答