logo标签列表

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

rdictionaryggplot2
5
空间数据：空间数据 收成数据：收成数据 代码：
    ## 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）绘制地图并将“ID_1”变量条目添加为地图上的标签，从而识别每个区域；

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

希望我清楚地表达了我的问题。

谢谢。

我下载了你的多边形数据，但无法在R或QGIS 2.0.1中读取。你确定文件没有损坏吗？ - jlhoward
@jlhoward：我非常确定它正在工作。我进行了双重检查，并且能够在R中加载文件并执行所有其他操作，将其转换为数据框以便与收益数据合并。我还能够复制上面的相同地图。 - iouraich
@smailov63 - 我相信你的原始文件没问题，但是你重新下载了上传的文件吗？有时候上传过程中会出现损坏。 - jlhoward
@jlhoward - 我重新上传了多边形数据文件，只是为了确保。说实话，我没有使用最初在这里上传的文件运行代码，因为我已经在项目目录中了。所以也许你是对的，当我上传它时出了些问题。如果重新上传解决了问题，请告诉我。谢谢。 - iouraich
@smailov63 - 上传的问题在于： "shapefile" 不是一个文件，而是一组文件，包括 .shp, .prj, .dbf, .csv, .shx。 readOGR（...） 查找所有这些文件。 我在这里找到了摩洛哥行政区域 shapefile，但是当我使用它时，您的代码无法运行。 您能否上传所有 Morocco_adm1.* 文件？ - jlhoward
@jlhoward - 我刚刚更新了空间数据。我上传了一个包含 .shp、.dbf 和 .shx 文件的压缩文件，以绘制多边形。希望这解决了问题。.csv 文件仅包含产量数据，我会在将其转换为数据框后与空间数据合并，就像代码中所做的那样。 - iouraich
1个回答
9

首先，让我为回复耽搁了这么久向您道歉——我错过了您的评论。这是您所想要的吗？

这个图是用以下代码生成的。在解释之前，您应该意识到创建图例是最小的问题。请注意两张地图中颜色的不同。您上面的代码没有将CO2变化分配给正确的区域。例如，根据MoroccoYields.csv，最大的变化（改善？）是在Region 4，达到了-0.205，但在您的地图上，最大的（最深的红色）在摩洛哥东北角，实际上是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行（地图上每条线端点的一行），并将其与具有14行（每个区域1行）的MoroccoYield组合。因此，cbind(...)将这些14行复制了7700次以填满所需的107,800行。表达式by =“ID_1”仅添加了另一列名为“by”，并将"ID_1"重复了107,800次。运行上面的语句，输入head(MoroccoReg.df)并查找最后一列。

那么如何进行合并？在R中有许多函数可以轻松完成此操作，但我无法让它们中的任何一个起作用。这就是实际工作的方法：

形状文件中的每个多边形都有一个ID。形状文件数据部分也有一个“ID_1”字段，但这些字段是不同且不相关的。[顺便说一下：形状文件数据部分中的ID_1字段和您的csv文件中的ID_1字段也不同：后者在区域编号前加上了"TR"，因此也必须处理。]使用以下方式重新排序形状文件：

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

这将改变多边形的顺序，但不会更改它们的ID。事实证明，多边形ID与shapefile数据部分中的行名称匹配，因此我使用cbind(...)将其添加到您的MoroccoYeild数据框中。

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

现在MoroccoYield有一个id字段，它映射到多边形ID，还有一个ID_1字段，用于识别地区。现在我们可以使用fortify(...)和merge(...)方法了。merge(...)方法需要一个by=参数。

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

追加所有 MoroccoYield 列到 MoroccoReg.df 相应的行中。

创建图例:

显而易见的问题是如何定位标签？理想情况下，我们会将来自 MoroccoYield$ID_1 的区域编号放置在每个区域的中心点，并创建一个基于 MoroccoYield$Label 标识区域的图例。

那么中心点在哪里找呢？这些信息存储在形状文件的多边形部分的一个不为人知的位置。长话短说，我创建了一个实用程序函数 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列。为了处理这个问题，在geom_text(...)中必须显式地将group设置为其他值。总之，这似乎有效。

没有必要道歉，实际上我应该向您道歉，因为从一开始我就提供了不完整的数据。这就是我所想的，如果不是更多的话：） 不过我有一个备注。当您检查使用修改后的代码绘制的地图版本时，我发现它上面存在一些差异，比如一些直线穿过多个区域（例如，一条直线穿过第二个区域和其他几个区域）。那么这是一个可以解决的问题吗？您认为是什么原因导致了这个问题？ - iouraich
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接