使用 R 开发地理专题地图答案

【问题标题】：Developing Geographic Thematic Maps with R使用 R 开发地理专题地图
【发布时间】：2010-11-18 15:48:49
【问题描述】：

显然，R 中有许多用于各种空间分析的包。这可以在CRAN Task View: Analysis of Spatial Data 中看到。这些包数量众多且种类繁多，但我想做的只是一些简单的thematic maps。我有县和州 FIPS 代码的数据，我有县和州边界的 ESRI 形状文件以及允许加入数据的随附 FIPS 代码。如果需要，形状文件可以很容易地转换为其他格式。

那么用 R 创建专题地图最直接的方法是什么？

这张地图看起来像是用 ESRI Arc 产品创建的，但这是我想用 R 做的事情：

alt text http://www.infousagov.com/images/choro.jpg 映射copied from here。

【问题讨论】：

请注意，这种类型的地图称为等值线，存在一些相当大的问题，即较小的地理区域往往有更多的人（例如东海岸与蒙大拿州），因此视觉外观是偏向人口密度低的地区。
另外，如果您正在处理原始 ESRI 数据，您可能会发现它有太多顶点。可以在github.com/hadley/data-counties/tree/master 找到在 R 中进行泛化的粗略尝试
@hadley，我完全同意你对 choropleths 的“问题”的看法。这通常是数据空间表示的问题。
@JDLong @hadley 如果你研究的是玉米而不是人，这可能没什么大不了的。
@leoniedu 你是正确的玉米。我必须处理的是这样的事情：大地理区域的生产历史波动率低于小地理区域的生产波动率，这仅仅是因为样本量。因此，当地理区域大小不同时，很难确定真正的生产者级别的“风险”。但这让它很有趣。 :)

标签： r map geolocation geospatial

【解决方案1】：

以下代码对我很有帮助。稍微定制一下，你就完成了。
_{（来源：eduardoleoni.com）}

library(maptools)
substitute your shapefiles here
state.map <- readShapeSpatial("BRASIL.shp")
counties.map <- readShapeSpatial("55mu2500gsd.shp")
## this is the variable we will be plotting
counties.map@data$noise <- rnorm(nrow(counties.map@data))

热图函数

plot.heat <- function(counties.map,state.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
          floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
          ,
          ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
          ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- heat.colors(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map,border=gray(.8), lwd=bw,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat))
  if (!is.null(state.map)) {
    plot(state.map,add=TRUE,lwd=1)
  }
  ##with(counties.map.c,text(x,y,name,cex=0.75))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
  ##title("Cartogram")
}

画出来

plot.heat(counties.map,state.map,z="noise",breaks=c(-Inf,-2,-1,0,1,2,Inf))

【讨论】：

哦，看起来真的很棒！我希望有人有这样的代码示例。谢谢！
在哪里可以获得 .shp 文件？我需要一个去荷兰的，但找不到
您好，是否可以添加一个外部变量而不是生成一个随机变量？直截了当，我想将预计收益率的百分比变化绘制为我的变量。因此，我是否需要修改包含地图坐标的形状文件以使其包含产量数据？

【解决方案2】：

我想我会在此处添加一些新信息，因为自发布以来围绕此主题开展了一些活动。这是革命博客上“Choropleth Map R Challenge”的两个很好的链接：

Choropleth Map R Challenge

Choropleth Challenge Results

希望这些对查看此问题的人有用。

一切顺利，

杰

【讨论】：

谢谢京东。现在有大量通过革命博客链接的地图绘制信息。
在哪里可以获得 .shp 文件？我需要一个去荷兰的，但找不到..

【解决方案3】：

查看包裹

library(sp)
library(rgdal)

这对地理数据很有用，并且

library(RColorBrewer)

对于着色很有用。 This map 是用上面的包和这个代码制作的：

VegMap <- readOGR(".", "VegMapFile")
Veg9<-brewer.pal(9,'Set2')
spplot(VegMap, "Veg", col.regions=Veg9,
 +at=c(0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5),
 +main='Vegetation map')

"VegMapFile" 是一个 shapefile，"Veg" 是显示的变量。可能可以通过一点工作做得更好。我似乎不允许上传图片，这是图片的链接：

【讨论】：

非常好的信息！谢谢。当您获得更多声誉时，您应该能够嵌入图像。
断开的链接 - 如果您可以更新它，那么也许有足够代表的人可以嵌入。

【解决方案4】：

查看 PBSmapping 包（参见 borh the vignette/manual and the demo）和 thisO'Reilly Data Mashups in R 文章（不幸的是，它不是免费的，但下载价值 4.99 美元，据 Revolutions blog 报道）。

【讨论】：

它 5 美元而且没有 DRM，这让我非常乐意仅在原则上下载。写得好，代码好，强烈推荐！

【解决方案5】：

只有三行！

library(maps);
colors = floor(runif(63)*657);
map("state", col = colors, fill = T, resolution = 0)

完成！！只需将第二行更改为 63 个元素的任意向量（每个元素在 0 到 657 之间，它们是 colors() 的成员）

现在如果你想变得花哨，你可以写：

library(maps);
library(mapproj);
colors = floor(runif(63)*657);
map("state", col = colors, fill = T, projection = "polyconic", resolution = 0);

这63个元素代表了63个地区，你可以通过运行得到它们的名字：

map("state")$names;

【讨论】：

【解决方案6】：

R 图形库有一个非常好的similar map，这应该是一个很好的起点。代码在这里：www.ai.rug.nl/~hedderik/R/US2004。您需要使用 legend() 函数添加图例。

【讨论】：

不错。我一直忘记图形库是此类示例的一个非常简洁的资源。
请注意，第一个链接已损坏。