r - 如何使用ggplot2使用带有描述相关标签的图例添加区域 map 的图例？

Question

SpatialPoly数据:SpatialData

产量数据:Yield Data

代码:

    ## Loading packages
    library(rgdal)
    library(plyr)
    library(maps)
    library(maptools)
    library(mapdata)
    library(ggplot2)
    library(RColorBrewer)
    library(foreign)  
    library(sp)

    ## Loading shapefiles and .csv files
    #Morocco <- readOGR(dsn=".", layer="Morocco_adm0")
    MoroccoReg <- readOGR(dsn=".", layer="Morocco_adm1")
    MoroccoYield <- read.csv(file = "F:/Purdue University/RA_Position/PhD_ResearchandDissert/PhD_Draft/Country-CGE/RMaps_Morocco/Morocco_Yield.csv", header=TRUE, sep=",", na.string="NA", dec=".", strip.white=TRUE)

    ## Reorder the data in the shapefile based on the category variable "ID_1" and change to dataframe
    MoroccoReg <- MoroccoReg[order(MoroccoReg$ID_1), ]
    MoroccoReg.df <- fortify(MoroccoReg)

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

    ## Check the structure and contents of shapefile
    attributes(MoroccoReg.df)

    ## Define new theme for map
    ## I have found this function on the website
    theme_map <- function (base_size = 12, base_family = "") {
    theme_gray(base_size = base_size, base_family = base_family) %+replace% 
    theme(
    axis.line=element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_blank(),
    axis.ticks=element_blank(),
    axis.ticks.length=unit(0.3, "lines"),
    axis.ticks.margin=unit(0.5, "lines"),
    axis.title.x=element_blank(),
    axis.title.y=element_blank(),
    legend.background=element_rect(fill="white", colour=NA),
    legend.key=element_rect(colour="white"),
    legend.key.size=unit(1.5, "lines"),
    legend.position="right",
    legend.text=element_text(size=rel(1.2)),
    legend.title=element_text(size=rel(1.4), face="bold", hjust=0),
    panel.background=element_blank(),
    panel.border=element_blank(),
    panel.grid.major=element_blank(),
    panel.grid.minor=element_blank(),
    panel.margin=unit(0, "lines"),
    plot.background=element_blank(),
    plot.margin=unit(c(1, 1, 0.5, 0.5), "lines"),
    plot.title=element_text(size=rel(1.8), face="bold", hjust=0.5),
    strip.background=element_rect(fill="grey90", colour="grey50"),
    strip.text.x=element_text(size=rel(0.8)),
    strip.text.y=element_text(size=rel(0.8), angle=-90) 
    )   
    }

    ## Plotting 

    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)
    #MoroccoRegMap <- MoroccoRegMap + scale_fill_gradient(low = "#CC0000",high = "#006600")
    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

结果:

问题:

在良率数据中，我有一列描述与“ID_1”列中的每个条目相对应的标签。我想要实现的是两件事:

1)绘制 map ，并在 map 上添加“ID_1”变量条目作为标签，从而识别每个区域；

2)生成第二个图例，该图例除了捕获数据中的值外，还包括数据帧中“ID_1”中的条目及其在“标签”列中的对应描述。

我希望我清楚地提出问题。

谢谢。

Answer 1

首先，我很抱歉花了这么长时间才回来-我在所有其他人中都错过了您的评论。这就是您的想法吗？

这是用以下代码生成的。在进行解释之前，您应该知道，创建图例是最少的问题。请注意，两个图中的颜色是如何不同的。上面的代码未将CO2更改分配给正确的区域。例如，根据MoroccoYields.csv，最大的更改(改进？)是-0.205中的Region 4，但是在您的 map 上最大的更改(最暗的红色)在摩洛哥的东北端，实际上是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

说明:

首先，在将 yield 数据与 map 区域合并时:

MoroccoReg.df <- cbind(MoroccoReg.df,MoroccoYield,by = 'ID_1')

这不是 cbind(...)的工作方式。 cbind(...)仅将其参数按列组合。它不是合并功能。因此，您有一个数据框 MoroccoReg.df，其中包含107,800行( map 上每个线端点一行)，并将其与 MoroccoYield组合，后者有14行(每个Region 1)。因此， cbind(...)将这14行复制7700次，以填充所需的107,800行。表达式 by="ID_1"仅添加了另一个名为“ by"的列，其中 "ID_1"已重复107800次。运行上面的语句并键入 head(MoroccoReg.df)并查找最后一列。

那么如何进行合并呢？ R中有许多功能可以简化此操作，但我无法使用其中的任何功能。这是起作用的:

shapefile中的每个多边形都有一个ID。 shapefile数据部分中也有一个“ID_1”字段，但它们是不同的且不相关。 [BTW:shapefile数据部分中的 ID_1字段与 ID_1文件中的 csv字段也不同:后者的 "TR"放在区域号之前；因此也必须处理]。
使用以下命令对shapefile重新排序:

MoroccoReg    <- MoroccoReg[order(MoroccoReg$ID_1), ]

将更改多边形的顺序，但不会更改其ID。事实证明，多边形ID与shapefile数据部分中的行名匹配，因此我在 cbind(...)数据框中添加了该前缀(使用 MoroccoYeild!)。

MoroccoYield  <- cbind(id=rownames(MoroccoReg@data),MoroccoYield)

因此，现在 MoroccoYield具有一个映射到多边形ID的 id字段和一个标识Region的 ID_1字段。现在我们可以 fortify(...)和 merge(...)了。 merge(...)确实接受了 by=参数。

MoroccoReg.df <- fortify(MoroccoReg)
MoroccoReg.df <- merge(MoroccoReg.df,MoroccoYield, by="id")

这会将所有 MoroccoYield列追加到 MoroccoReg.df的相应行。

创建图例:

显而易见的问题是如何放置标签？理想情况下，我们将 MoroccoYield$ID_1中的Region号放在每个区域的质心处，然后根据 MoroccoYield$Label创建一个图例来标识Regions。

那么在哪里可以找到质心呢？它们存储在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)

...创建 map 。我找不到一种干净的方法来使用正式的“ggplot图例”来执行此操作，因此我不得不使用 annotate(...)。定位注释有点麻烦，但它似乎可以工作。

编辑:为了响应@ 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)

...你得到这个:

我相信可以解决此问题。 map 中多余线条的原因是 ggplot必须按列 MoroccoReg.df$group分组(因此， aes(..., group=group)而不是 aes(...,group=id))。但是，执行此操作时， ggplot尝试在所有图层中按 "group"分组。在 geom_text(...)中，我们在其中使用新的本地数据集- labs数据框-没有 group列。为了解决这个问题，我们必须在 group中显式地将 geom_text(...)设置为其他内容。底线:这似乎可行。