首先，让我为花了这么长时间才回来表示歉意 - 我错过了您在所有其他评论中的评论。这就是你的想法吗？
这是使用以下代码生成的。在进行解释之前，您应该意识到创建图例是最不重要的问题。请注意两张地图中的颜色有何不同。您上面的代码没有将二氧化碳变化分配到正确的区域。例如，根据MoroccoYields.csv，最大的变化（改进？）是-0.205 in Region 4，但在你的地图上最大的（最深红色的）位于摩洛哥的东北端，实际上是l'Oriental (Region 6)。代码后面有解释。
## Loading packages
library(rgdal)
library(plyr)
library(maps)
library(maptools)
library(mapdata)
library(ggplot2)
library(RColorBrewer)
library(foreign)
library(sp)
# get.centroids: function to extract polygon ID and centroid from shapefile
get.centroids = function(x){
poly = MoroccoReg@polygons[[x]]
ID = poly@ID
centroid = as.numeric(poly@labpt)
return(c(id=ID, long=centroid[1], lat=centroid[2]))
}
#setwd("Directory where shapefile and Yields are stored")
## Loading shapefiles and .csv files
MoroccoReg <- readOGR(dsn=".", layer="Morocco_adm1")
MoroccoYield <- read.csv(file = "Morocco_Yield.csv", header=TRUE, sep=",", na.string="NA", dec=".", strip.white=TRUE)
MoroccoYield$ID_1 <- substr(MoroccoYield$ID_1,3,10)
## Reorder the data in the shapefile based on the category variable "ID_1" and change to dataframe
MoroccoReg <- MoroccoReg[order(MoroccoReg$ID_1), ]
MoroccoYield <- cbind(id=rownames(MoroccoReg@data),MoroccoYield)
# build table of labels for annotation (legend).
labs <- do.call(rbind,lapply(1:14,get.centroids))
labs <- merge(labs,MoroccoYield[,c("id","ID_1","Label")],by="id")
labs[,2:3] <- sapply(labs[,2:3],function(x){as.numeric(as.character(x))})
labs$sort <- as.numeric(as.character(labs$ID_1))
labs <- labs[order(labs$sort),]
MoroccoReg.df <- fortify(MoroccoReg)
## This does NOT work...
## Add the yield impacts column to shapefile from the .csv file by "ID_1"
## Note that in the .csv file, I just added the column "ID_1" to match it with the shapefile
#MoroccoReg.df <- cbind(MoroccoReg.df,MoroccoYield,by = 'ID_1')
## Do it this way...
MoroccoReg.df <- merge(MoroccoReg.df,MoroccoYield, by="id")
## Check the structure and contents of shapefile
attributes(MoroccoReg.df)
## Plotting
MoroccoRegMap1 <- ggplot(data = MoroccoReg.df, aes(long, lat, group=id))
MoroccoRegMap1 <- MoroccoRegMap1 + geom_polygon(aes(fill = A2Med_noCO2))
MoroccoRegMap1 <- MoroccoRegMap1 + geom_path(colour = 'gray', linestyle = 2)
MoroccoRegMap1 <- MoroccoRegMap1 + scale_fill_gradient2(name = "%Change in yield",low = "#CC0000",mid = "#FFFFFF",high = "#006600")
MoroccoRegMap1 <- MoroccoRegMap1 + labs(title="SRES_A2, noCO2 Effect")
MoroccoRegMap1 <- MoroccoRegMap1 + coord_equal() #+ theme_map()
MoroccoRegMap1 <- MoroccoRegMap1 + geom_text(data=labs, aes(x=long, y=lat, label=ID_1), size=4)
MoroccoRegMap1 <- MoroccoRegMap1 + annotate("text", x=max(labs$long)-5, y=min(labs$lat)+3-0.5*(1:14),
label=paste(labs$ID_1,": ",labs$Label,sep=""),
size=3, hjust=0)
MoroccoRegMap1
解释：
首先，将产量数据与地图区域合并：您使用
MoroccoReg.df <- cbind(MoroccoReg.df,MoroccoYield,by = 'ID_1')
这不是如何cbind(...) works. cbind(...)只是按列组合它的参数。它不是合并函数。所以你有一个数据框，MoroccoReg.df，有 107,800 行（地图上的每个线端点一行），并且您将其与MoroccoYield，有 14 行（每个区域 1 行）。所以cbind(...)将这 14 行复制 7700 次以填充所需的 107,800 行。表达方式by="ID_1"仅添加另一列名为“by" with "ID_1"重复107,800次。运行上面的语句并输入head(MoroccoReg.df)并查找最后一列。
那么如何进行合并呢？里面有很多函数R这应该会让这变得容易，但我无法让它们中的任何一个工作。这是有效的：
shapefile 中的每个多边形都有一个 ID。 shapefile 数据部分中还有一个“ID_1”字段，但这些是不同且不相关的。 [顺便说一句：ID_1shapefile 数据部分中的字段，以及ID_1你的领域csv文件也不同：后者有"TR"前置到区域编号；所以这也必须得到处理]。
使用以下命令重新排序 shapefile：
MoroccoReg <- MoroccoReg[order(MoroccoReg$ID_1), ]
将更改多边形的顺序，但不会更改其 ID。事实证明，多边形 ID 与 shapefile 数据部分中的行名称匹配，因此我在前面添加了它（使用cbind(...)！）到你的MoroccoYeild数据框。
MoroccoYield <- cbind(id=rownames(MoroccoReg@data),MoroccoYield)
So now MoroccoYield has an id映射到多边形 ID 的字段，以及ID_1字段，标识区域。现在我们可以fortify(...) and merge(...). merge(...)确实需要一个by=争论。
MoroccoReg.df <- fortify(MoroccoReg)
MoroccoReg.df <- merge(MoroccoReg.df,MoroccoYield, by="id")
这会附加您的所有MoroccoYield列到适当的行MoroccoReg.df.
创建传奇：
显而易见的问题是如何放置标签？理想情况下，我们会将区域编号放在MoroccoYield$ID_1在每个区域的质心，然后创建一个图例来标识区域，基于MoroccoYield$Label.
那么哪里可以找到质心呢？它们存储在 shapefile 的多边形部分的一个不起眼的位置。长话短说，我创建了一个实用函数get.centroid(...)从多边形中提取质心。然后我将该函数应用于所有多边形，以生成具有相应多边形 ID 的质心表。然后我将其与标签合并MoroccoYield。这创建了一个数据框labs其中有以下列：
id: polygon ID
long: centroid longitude
lat: centroid latitude
ID_1: region ID
label: region name
sort: a sortable (numeric) version of ID_1
然后，将以下代码添加到您的 ggplot...
...
MoroccoRegMap1 <- MoroccoRegMap1 + geom_text(data=labs, aes(x=long, y=lat, label=label.id), size=4)
MoroccoRegMap1 <- MoroccoRegMap1 + annotate("text", x=max(labs$long)-5, y=min(labs$lat)+3-0.5*(1:14),
label=paste(labs$label.id,": ",labs$Label,sep=""),
size=3, hjust=0)
...创建地图。我找不到干净的方法来使用正式的“ggplot legend”来做到这一点，所以我不得不使用annotate(...)。定位注释是一种黑客行为，但它似乎有效。
Edit:作为对 @smailov83 的评论的回应，如果您将创建 ggplot 的代码更改为此...
MoroccoRegMap1 <- ggplot(data = MoroccoReg.df, aes(long, lat, group=group))
MoroccoRegMap1 <- MoroccoRegMap1 + geom_polygon(aes(fill = A2Med_noCO2))
MoroccoRegMap1 <- MoroccoRegMap1 + geom_path(colour = 'gray', linestyle = 2)
MoroccoRegMap1 <- MoroccoRegMap1 + scale_fill_gradient2(name = "%Change in yield",low = "#CC0000",mid = "#FFFFFF",high = "#006600")
MoroccoRegMap1 <- MoroccoRegMap1 + labs(title="SRES_A2, noCO2 Effect")
MoroccoRegMap1 <- MoroccoRegMap1 + coord_equal() #+ theme_map()
MoroccoRegMap1 <- MoroccoRegMap1 + geom_text(data=labs, aes(x=long, y=lat, label=ID_1, group=ID_1), size=4)
MoroccoRegMap1 <- MoroccoRegMap1 + annotate("text", x=max(labs$long)-5, y=min(labs$lat)+3-0.5*(1:14),
label=paste(labs$ID_1,": ",labs$Label,sep=""),
size=3, hjust=0)
...你得到这个：
我相信这可以解决问题。地图中额外线条的原因是ggplot必须按列分组MoroccoReg.df$group (so, aes(..., group=group) not aes(...,group=id)）。然而，当你这样做时，ggplot尝试分组依据"group"在所有层中。在geom_text(...)，我们使用一个新的本地数据集 -labs数据框 - 没有group柱子。为了解决这个问题，我们必须明确设置group到其他东西geom_text(...)。底线：这似乎有效。
