11

我们每天制作地图，显示我们地区 30 个不同区域的温度计算水平，每个区域根据水平填充不同的颜色。这张地图看起来像

在此处输入图像描述

现在我想将地图生成切换到 R。我已经下载了省和市边界（你可以找到整个西班牙的边界或在这里我所在地区的子集）并按照 Hadley 的示例设法用 ggplot2 绘制它们。

我还可以生成一个包含两列的 ascii 文件：标识符 (CODINE) 和每日级别。你可以在这里下载。

这是我第一个尝试使用 R 和 ggplot2 绘制 shapefile 的脚本，因此可能会出现错误，并且可以肯定它可以改进，欢迎提出建议。以下代码（基于前面提到的 Hadley）对我有用：

> require("rgdal")
> require("maptools")
> require("ggplot2")
> require("plyr")

# Reading municipal boundaries

esp = readOGR(dsn=".", layer="lineas_limite_municipales_etrs89")

muni=subset(esp, esp$PROV1 == "46" | esp$PROV1 == "12" | esp$PROV1 == "3")
muni@data$id = rownames(muni@data)
muni.points = fortify(muni, region="id")
muni.df = join(muni.points, muni@data, by="id")

# Reading province boundaries

prov = readOGR(dsn=".", layer="poligonos_provincia_etrs89")

pr=subset(prov, prov$CODINE == "46" | prov$CODINE == "12" | prov$CODINE == "03" )
pr@data$id = rownames(pr@data)
pr.points = fortify(pr, region="id")
pr.df = join(pr.points, pr@data, by="id")

ggplot(muni.df) + aes(long,lat,group=group) + geom_path(color="blue") +
+ coord_equal()+ geom_path(data=pr.df, + 
aes(x=long, y=lat, group=group),color="red", size=0.5)

这段代码绘制了一张带有所有边界的漂亮地图在此处输入图像描述

对于按级别填充的多边形，我尝试按照http://tormodboe.wordpress.com/2011/02/22/g%C3%B8y-med-kart-2/中的建议进行读取然后合并

level=read.csv("levels.dat",header=T,sep="")
munlevel=merge(muni.df,level,by="CODINE")

但它给出了一个错误

错误 en fix.by(by.x, x) : 'by' 必须指定唯一有效的列

我不熟悉 shapefile，也许我需要了解更多关于 shp 数据属性的信息，才能找到合并两个数据集的正确选择。如何合并数据以便绘制线条（市政边界）然后用级别填充？

4

1 回答 1

18

[注意：这个问题是在一个多月前提出的，所以 OP 可能已经找到了一种不同的方法来解决他们的问题。我在处理这个相关问题时偶然发现了它。包含此答案是希望对其他人有所帮助。]

这似乎是OP所要求的......

...并使用以下代码生成：

require("rgdal")
require("maptools")
require("ggplot2")
require("plyr")

# read temperature data
setwd("<location if your data file>")
temp.data        <- read.csv(file = "levels.dat", header=TRUE, sep=" ", na.string="NA", dec=".", strip.white=TRUE)
temp.data$CODINE <- str_pad(temp.data$CODINE, width = 5, side = 'left', pad = '0')

# read municipality polygons
setwd("<location of your shapefile")
esp     <- readOGR(dsn=".", layer="poligonos_municipio_etrs89")
muni    <- subset(esp, esp$PROVINCIA == "46" | esp$PROVINCIA == "12" | esp$PROVINCIA == "3")
# fortify and merge: muni.df is used in ggplot
muni@data$id <- rownames(muni@data)
muni.df <- fortify(muni)
muni.df <- join(muni.df, muni@data, by="id")
muni.df <- merge(muni.df, temp.data, by.x="CODIGOINE", by.y="CODINE", all.x=T, a..ly=F)
# create the map layers
ggp <- ggplot(data=muni.df, aes(x=long, y=lat, group=group)) 
ggp <- ggp + geom_polygon(aes(fill=LEVEL))         # draw polygons
ggp <- ggp + geom_path(color="grey", linestyle=2)  # draw boundaries
ggp <- ggp + coord_equal() 
ggp <- ggp + scale_fill_gradient(low = "#ffffcc", high = "#ff4444", 
                                 space = "Lab", na.value = "grey50",
                                 guide = "colourbar")
ggp <- ggp + labs(title="Temperature Levels: Comunitat Valenciana")
# render the map
print(ggp)

解释：

导入到 RreadOGR(...)中的 shapefile 具有类型SpacialDataFrame并具有两个主要部分：一个多边形部分，其中包含每个多边形上所有点的坐标，以及一个包含每个多边形信息的数据部分（因此，每个多边形一行）。这些可以被引用，例如，使用muni@polygons和muni@data。实用程序函数fortify(...)将多边形部分转换为组织用于绘图的数据框ggplot。所以基本的工作流程是：

[1] Import temperature data file (temp.data)
[2] Import polygon shapefile of municipalities (muni)
[3] Convert muni polygons to a data frame for plotting (muni.df <- fortify(...))
[4] Join columns from muni@data to muni.df
[5] Join columns from temp.data to muni.df
[6] Make the plot

连接必须在公共字段上完成，这是大多数问题出现的地方。原始 shapefile 中的每个多边形都有一个唯一的 ID 属性。在 shapefile 上运行fortify(...)会创建一个id基于此的列。但是数据部分没有 ID 列。相反，多边形 ID 存储为行名称。所以首先我们必须添加一个id列，muni@data如下所示：

muni@data$id <- rownames(muni@data)

现在我们有一个id字段 inmuni@data和一个对应的id字段 in muni.df，所以我们可以进行连接：

muni.df <- join(muni.df, muni@data, by="id")

要创建地图，我们需要根据温度水平设置填充颜色。为此，我们需要从 to加入LEVEL列。在其中有一个标识自治市的字段。现在， 中也有一个相应的字段。但是有一个问题：is ，带有前导零，而is integer 这意味着缺少前导零（也许是从 Excel 导入的？）。因此，仅加入这两个字段不会产生匹配项。我们必须先转换成前导零：temp.datamuni.dftemp.dataCODINECODIGOINEmuni.dfCODIGOINEchar(5)CODINECODINEchar(5)

temp.data$CODINE <- str_pad(temp.data$CODINE, width = 5, side = 'left', pad = '0')

现在我们可以temp.dat根据muni.df相应的字段加入。

muni.df <- merge(muni.df, temp.data, by.x="CODIGOINE", by.y="CODINE", all.x=T, a..ly=F)

我们使用merge(...)而不是join(...)因为连接字段具有不同的名称并且join(...)要求它们具有相同的名称。（请注意，但是这join(...)更快，如果可能的话应该使用）。因此，最后，我们有一个数据框，其中包含用于绘制多边形的所有信息和LEVEL可用于为每个多边形建立填充颜色的温度。

关于 OP 原始代码的一些注释：

  1. OP 的第一张地图（顶部的绿色地图）标识了“我们地区的 30 个不同区域......”。我找不到识别这些区域的 shapefile。市政府文件确定了 543 个市镇，我看不出有办法将它们分成 30 个区域。此外，温度级别文件有 542 行，每个市（或多或少）一个。

  2. OP 正在为市政当局导入线形文件以绘制边界。您不需要它，因为geom_polygon(...)它将绘制（并填充）多边形并geom_path(...)绘制边界。

于 2013-12-09T16:13:59.777 回答