如何使用ggplot2为描述相关标签的图例添加区域地图图例？

Question

SpatialPoly数据:SpatialData

产量数据:产量数据

码:

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

结果:

题:

在Yield数据中,我有一列描述了与"ID_1"列中每个条目对应的标签.我想要实现的是两件事:

1)绘制地图并在地图上添加"ID_1"变量条目作为标签,从而识别每个区域;

2)除了捕获数据中的值之外,还生成第二个图例,以及"ID_1"中的条目及其在数据帧的"标签"列中的相应描述.

我希望我能清楚地提出问题.

谢谢.

Answer 1

首先,让我为花了这么长时间才能回来道歉 - 我在所有其他人中错过了你的评论.这是你的想法吗？

这是使用以下代码生成的.在进行解释之前,您应该意识到创建一个图例是您遇到的最少问题.注意两个地图中的颜色是如何不同的.上面的代码不会将CO2更改分配给正确的区域.例如,根据MoroccoYields.csv,最大的变化(改进？)是-0.205在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(...)有效.cbind(...)只是按列方式结合它的论点.它不是合并功能.所以你有一个数据框,MoroccoReg.df有107,800行(地图上每个行端点都有一行),你要把它组合起来MoroccoYield,它有14行(每个Region有1行).因此,cbind(...)将这14行重复7700次以填充所需的107,800行.该表达式by="ID_1"仅添加了另一个名为" by"with "ID_1"replicated 107,800次"的列.运行上面的语句并键入head(MoroccoReg.df)并查找最后一列.

那么合并怎么做？有许多功能R可以使这很容易,但我无法让它们中的任何一个工作.这是做了什么工作:

shapefile中的每个多边形都有一个ID.shapefile数据部分中还有一个"ID_1"字段,但这些字段不同且不相关.[BTW: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)

所以现在MoroccoYield有一个id映射到多边形ID的ID_1字段和一个标识Region 的字段.现在我们可以fortify(...)和merge(...).merge(...)确实需要by=争论.

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

这会将所有MoroccoYield列附加到适当的行MoroccoReg.df.

创建图例:

显而易见的问题是如何定位标签？理想情况下,我们会将Region编号从MoroccoYield$ID_1每个区域的质心放置,然后创建一个标识区域的图例,基于MoroccoYield$Label.

那么在哪里可以找到质心？它们存储在shapefile的polygon部分中的一个模糊位置.总而言之,我创建了一个实用函数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传奇"做到这一点,所以我不得不使用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)

你得到这个:

我认为这解决了这个问题.地图中额外行的原因是ggplot必须按列分组MoroccoReg.df$group(因此,aes(..., group=group) 不是 aes(...,group=id)).但是,执行此操作时,会ggplot尝试"group"在所有图层中进行分组.在geom_text(...)我们使用新的本地数据集(labs数据框)的地方,没有group列.为了解决这个问题,我们必须明确地设置group其他内容geom_text(...).底线:这似乎有效.