用ggplot2在美国的专题地图上重新定位阿拉斯加和夏威夷 [英] Relocating Alaska and Hawaii on thematic map of the USA with ggplot2

查看:278
本文介绍了用ggplot2在美国的专题地图上重新定位阿拉斯加和夏威夷的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我试图制作一张显示美国所有50个州的专题地图,但我无法以可靠的方式重新安置阿拉斯加和夏威夷。我有几个想法,但没有一个能很好地工作。我会现在演示它们。



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

  setwd(tempdir())
download.file(https://dl.dropbox.com/s/wl0z5rpygtowqbf/states_21basic.zip?dl=1,
usmapdata .zip,
method =curl)
#这是http://www.arcgis.com/home/item.html的镜像?
#id = f7f805eb65eb4ab787a0a3e1116ca7e5
unzip(usmapdata.zip)

require(rgdal)
all_states< - readOGR(states_21basic /,states )

require(ggplot2);要求(maptools);要求(rgeos);要求(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))),
color =white,size = 0.25
)+ coord_map(projection =azequalarea)+
scale_fill_gradient(limits = c(1,50))

现在我们删除所有背景当我们重叠非连续状态时,它们不会发生冲突:

pre $ p < - p + theme(axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
轴。 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())



使用视口



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

 <$ c $ (all_states,id =none)
HI <-p%+%subset(all_states,id ==Alaska)+ theme(legend.position =none =Hawaii)+ theme(legend.position =none)
连续< -p%+%subset(all_states,id!=Alaska&
$ b 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<  - (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



这并不令人满意,因为阿拉斯加通常不会按比例大多数美国地图看起来非常大。

问题



有没有人会有更好的表现?方法?

解决方案

以下是如何通过投影和转换来完成的。您需要:

  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(夏威夷)< - proj4string(usa)

usa = usa [!美国$ STATE_NAME%%c(阿拉斯加,夏威夷),]
美国= rbind(美国,阿拉斯加,夏威夷)

美元(美元)


$ b $ fix1< - function(object,params){
r = params [1]; scale = params [2]; shift = params [3:4]
对象= 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,现在转换为等面积,并运行fixup函数:  


现在转换到等面积, (美国,CRS(+ init = epsg:2163))
usfix = fixup(usAEA, c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix)

这些参数分别是阿拉斯加和夏威夷的旋转,缩放,x和y偏移,并且通过试验和错误获得。仔细调整它们。即使将夏威夷的比例参数更改为0.99999,也会因为涉及大量数据而将其从地球上传送出去。



如果您想将其恢复为lat-long:



$ p $ us code usfixLL = spTransform(usfix,CRS(+ init = epsg:4326))
plot(usfixLL)

但我不确定是否需要在 ggplot ,因为我们已经使用 spTransform 完成了这项工作。



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


I am trying to create a thematic map showing all 50 US states, but I am having trouble relocating Alaska and Hawaii in a reliable way. I have a couple of ideas but none of them work well. I will demonstrate them now.

First we need to import the data; using the data in the maps package is not enough because it does not include Hawaii and 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")

Now we define some plot aesthetics:

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

Now we remove all background etc so they don't clash when we overlap the non-contiguous states:

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

Using viewports

The first idea I had was to use viewports:

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)

This looks nice but it is not satisfactory because it is very sensitive to slight changes in the figure, for example resizing or changes in the size and shape of the legend.

Manually moving Alaska and Hawaii

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

This is not satisfactory because Alaska is usually not to scale on most US maps so it looks very big. Also relocating Alaska and Hawaii changes the distortion introduced by the map projection.

Question

Does anyone have any better approaches?

解决方案

Here's how to do it by projecting and transforming. You will need:

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
}

Then read in your shapefile. Use rgdal:

us = readOGR(dsn = "states_21basic",layer="states")

Now transform to equal-area, and run the fixup function:

usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix)

The parameters are rotations, scaling, x and y shift for Alaska and Hawaii respectively, and were obtained by trial and error. Tweak them carefully. Even changing Hawaii's scale parameter to 0.99999 sent it off the planet because of the large numbers involved.

If you want to turn this back to lat-long:

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

But I'm not sure if you need to use the transformations in ggplot since we've done that with spTransform.

You can now jump through the ggplot2 fortify business. I'm not sure if it matters for you but note that the order of the states is different in the usfix version - Alaska and Hawaii are now the last two states.

这篇关于用ggplot2在美国的专题地图上重新定位阿拉斯加和夏威夷的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆