The Data
The code
#
# This is code for mapping of CGE_Morocco results
#
# rm(list = ls(all = TRUE)) # don't use this in code that others will copy/paste
## Loading packages
library(rgdal)
library(plyr)
library(maps)
library(maptools)
library(mapdata)
library(ggplot2)
library(RColorBrewer)
## Loading shape files administrative coordinates for Morocco maps
#Morocco <- readOGR(dsn=".", layer="Morocco_adm0")
MoroccoReg <- readOGR(dsn=".", layer="Morocco_adm1")
## Reorder the data in the shapefile based on the regional order
MoroccoReg <- MoroccoReg[order(MoroccoReg$ID_1), ]
## Add the yield impacts column to shapefile
MoroccoReg@data$SRESB2_CO2 <- c(0.003,0.100,0.116,-0.105,-0.010,0.048,0.006,-0.004,0.061,0.032,0.003,-0.016,-0.018,0.095)
## Check the structure and contents of shapefile
summary(MoroccoReg)
attributes(MoroccoReg)
## Plotting
MoroccoRegMap <- ggplot(data = MoroccoReg, aes(long, lat, group = group))
MoroccoRegMap <- MoroccoRegMap + geom_polygon()
MoroccoRegMap <- MoroccoRegMap + geom_path(colour = 'gray', linestyle = 2)
MoroccoRegMap <- MoroccoRegMap + scale_fill_brewer('ID_1')
MoroccoRegMap <- MoroccoRegMap + coord_equal() + theme_bw()
MoroccoRegMap
The plot
The Question is two fold:
First, I am looking to be able to color each region separately.
Second, I want to do that based on each region's projected yield impact as captured by the variable "SRESB2_CO2" in the data.
Thanks in advance for the help.
Here is the additional yield data in .csv format, code and results.
Yield data:
Data
Code:
## 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]))
}
## 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)
MoroccoYield$ID_1 <- substr(MoroccoYield$ID_1,3,10) # Eliminate the ID_1 column
## Reorder the data in the shapefile
MoroccoReg <- MoroccoReg[order(MoroccoReg$ID_1), ]
MoroccoYield <- cbind(id=rownames(MoroccoReg@data),MoroccoYield) # Add the column "id" for correct merging with the spatial data
## 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) # transform shapefile into dataframe
MoroccoReg.df <- merge(MoroccoReg.df,MoroccoYield, by="id") # merge the spatial data and the yield data
str(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.2, "lines"),
legend.position="right",
legend.text=element_text(size=rel(0.8)),
legend.title=element_text(size=rel(0.8), 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.2)),
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)
)
}
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 <- MoroccoRegMap1 + geom_text(data=labs, aes(x=long, y=lat, label=ID_1, group=ID_1), size=6)
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=5, hjust=0)
MoroccoRegMap1
Results: