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"中的条目及其在数据帧的"标签"列中的相应描述.
我希望我能清楚地提出问题.
谢谢.
首先,让我为花了这么长时间才能回来道歉 - 我在所有其他人中错过了你的评论.这是你的想法吗?
这是使用以下代码生成的.在进行解释之前,您应该意识到创建一个图例是您遇到的最少问题.注意两个地图中的颜色是如何不同的.上面的代码不会将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_1
shapefile数据部分中的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(...)
.底线:这似乎有效.