R ggplot2与shapefile和csv数据合并以填充多边形

Question

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

现在我想将地图生成切换为R.我已经下载了省级和市级边界(您可以找到整个西班牙的边界或这里是我所在地区的子集)并设法根据Hadley的示例使用ggplot2绘制它们.

我还可以生成一个包含两列的ascii文件:identifier(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数据属性的信息,以找到合并两个数据集的正确选择.我如何合并数据,以便我可以绘制线条(市政边界),然后用水平填充它？

Answer 1

[注意:这个问题在一个月前被问到,所以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)

说明:

导入到R中的Shapefile具有readOGR(...)类型SpacialDataFrame并且具有两个主要部分:包含每个多边形上所有点的坐标的ploygon部分,以及包含关于每个多边形的信息的数据部分(因此,每个多边形一行).这些可以参考,例如,使用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属性.fortify(...)在shapefile上运行会创建一个id基于此的列.但数据部分中没有ID列.相反,多边形ID存储为行名称.首先,我们必须添加一id列muni@data如下:

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

现在我们有一个id字段muni@data和相应的id字段muni.df,所以我们可以进行连接:

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

要创建地图,我们需要根据温度水平设置填充颜色.要做到这一点,我们需要加入LEVEL从列temp.data到muni.df.在temp.data这里有一个CODINE识别市政当局的领域.还有,现在,相应字段CODIGOINE在muni.df.但是有一个问题:CODIGOINE是char(5),前导零,而是CODINE整数,这意味着缺少前导零(从Excel导入,也许？).所以加入这两个字段就不会产生匹配.我们必须首先转换CODINE为char(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(...)绘制边界.