Most of the time, we need to deal with complex systems that have several components, each with different properties and behavior. Usually, these properties vary with time and space, and understanding their spatiotemporal dynamics plays a key role in our understanding of the system performance and its vulnerabilities. Therefore, aggregation in time and space would hide variations that might be important in our analysis and understanding of the system behavior. In this blog post, I am going to show how we can add bar plots onto a map at different locations for better visualization of variables. Some examples for when this type of visualization is useful are irrigation districts that have different water rights or cropping patterns that affect their crop production and water requirements in dry or wet years; another example would be the sensitivity indices of specific model parameters that are clustered based on their location and varying based on the wetness and dryness of the system. To demonstrate this, I am going to create an animated graph that shows annual crop yield variations in different counties in the State of Washington. You can download the Washington counties’ data layer (WA_County_Boundaries-shp.zip) and NASS historical yield data for corn, winter wheat, and spring wheat (NASS.txt) from here.
library(rgdal)
shp_counties <- readOGR(dsn = file.path("(your directory)/WA_County_Boundaries-shp/WA_County_Boundaries.shp"))
“shp_countiesis” is a SpatialPolygonsDataFrame object and has a different format than a regular dataframe. It usually has a few slots that contain different type of information. For example, “data” has non-geographic properties. We can explore the information for each polygon with:
head(shp_counties@data)
#The "JURISDIC_2" and "JURISDIC_3" columns both contain the names of counties.
To visualize it with ggplot, it has to first be converted into a dataframe. Then, we can use it with geom_polygon function:
library(ggplot2)
counties <- fortify(shp_counties, region="JURISDIC_2")
head(counties)
#long lat order hole piece id group
#1 -13131351 5984710 1 FALSE 1 Adams Adams.1
#2 -13131340 5983102 2 FALSE 1 Adams Adams.1
ggplot() + geom_polygon(data = counties, aes(x = long, y = lat, group = group), colour = "dark red", fill = NA)
Now, I am going to extract the center of each polygon so that I can later add the bar plots to these coordinates:
library(rgeos)
counties_centroids<- as.data.frame(gCentroid(shp_counties, byid=T))
We are going to extract just a few years of data from the yield dataset to show on the map:
library(data.table)
yield<- fread("(your directory)/NASS.txt")
# Year Ag_District Ag_District_Code Data_Item Value County JURISDIC_5 Group
#1: 1947 EAST CENTRAL 50 CORN, GRAIN 41.0 Adams 53001 CORN
#2: 1948 EAST CENTRAL 50 CORN, GRAIN 40.0 Adams 53001 CORN
years<- as.data.frame(unique(yield$Year))
years<- as.data.frame(years[34:42,])
For the final map, I am going to need a common legend for all of the bar plots in all of the counties and years in which I am interested. So, I need to know all of the categories that are available:
unique(yield$Group)
unique(yield$Data_Item)
Based on the “Group” column, we know that there are three main groups of crops: corn, winter wheat and spring wheat. Based on the “Data_Item” column, we know that there are three different types of winter and spring wheat (non-irrigated, non-irrigated with continuous cropping (CC), and non-irrigated with summer fallow (SF)). Note that we do not have all of these crop types for all of the counties and all the years, and the common legend should be created from a location that all the crop types are available, so that it can be applicable for all of the counties and years that I am eventually going to plot. For this reason, I am going to subset a dataset and create a bar plot for one county and year when all of the crop types are available:
sample_yield<- subset(yield, Year=="1981" & County=="Adams")
I want to use custom colors so that each crop Group has a different color and each Data_item corresponding to a Group share the same Group color theme, which gradually changes from darker to brighter. First, I create three color functions, each corresponding to 3 Groups in my dataset.
colfunc1 <- colorRampPalette(c("darkRED"))
colfunc2 <- colorRampPalette(c("darkBLUE", "lightBLUE"))
colfunc3 <- colorRampPalette(c("darkGREEN", "lightgreen"))
Now I am going to use these functions to create a color code for each Group and save it in “colors”:
colors<-c(colfunc1(nrow(subset(sample_yield,sample_yield$Group=="CORN"))),colfunc2(nrow(subset(sample_yield,sample_yield$Group=="SW"))),colfunc3(nrow(subset(sample_yield,sample_yield$Group=="WW"))))
Next, in the template bar plot, I can use the customized “colors” in scale_fill_manual() function:
ggplot(sample_yield,aes(x=Group,y=Value,fill=factor(Data_Item)))+
scale_fill_manual(values=colors)+
geom_bar(stat='identity', color = "black")
From this plot, we just need to extract the legend. Also, we need to change its font size since we are going to add it to another plot. You should try different sizes to find out which one is more appropriate for your dataset/figure. I am saving this plot in “tmp” while I remove the legend title and change the font size. Then, I extract the legend section by using the get_legend() function and convert it to a ggplot by using the as_ggplot() function.
tmp <- ggplot(sample_yield,aes(x=Group,y=Value,fill=factor(Data_Item)))+
scale_fill_manual(values=colors)+
geom_bar(stat='identity', color = "black")+theme(legend.title = element_blank(),legend.text =element_text(size=19,face="bold",colour="black") )
library(cowplot)
legend<- get_legend(tmp)
library(ggpubr)
tmp_legend<- as_ggplot(legend)
Now, I am going to save the counties map with the appropriate font size and title and add the extracted legend from the yield data to it with the geom_subview() funtion. You need to adjust the coordinates of that you want to show the legend on the map based on your data.
tmp_map <- ggplot(counties)+
geom_polygon(aes(long, lat, group=group), fill="white")+
geom_path(color="black", mapping=aes(long, lat, group=group), size=1.5)+
theme(axis.text.y=element_text(size=17,face="bold",colour="black"),
axis.text.x=element_text(size=17,face="bold",colour="black"),
axis.title.y =element_text(size=17,face="bold",colour="black"),
axis.title.x =element_text(size=17,face="bold",colour="black"),
plot.title =element_text(size=25,face="bold",colour="black"))+
labs(title = "NASS Yield Data (1980-1995)")
library(ggimage)
tmp_map_final<- tmp_map+ geom_subview(x=-13830000, y=5730000, subview=tmp_legend)
In the below section, I am going to show how we can loop through all of the years and then all the counties in the county map, use the center coordinate of each county (polygon), plot the corresponded bar plot, and print it in the center of each polygon. We can use the ggplotGrob() funtion to extract different pieces of the bar plot created with gpplot. With this function, we can treat each part of the bar plot (background, axis, labels, main plot panel, etc.) as a graphical object and move it to the coordinates that are in our interest. For example, if we just want to use the main panel and not any other components, we can extract the panel, and adjust the other components as we wish to present in the final graph.
In this example, I am adding all of the bar plots for all counties in one year as a graphical object in a list “barplot_list”. Then, by using the annotation_custom() function I add each item in the “barplot_list” to the coordinates at the center of the polygons (counties) that I already extracted. Note that the orders of center coordinates and plots in the “barplot_list” are the same.
At the end, I just add the base map with the customized legend (“tmp_map_final”) and with the list of all bar plots with their customized locations (“barplot_annotation_list”). Then, I add them all in a list (“all_list“) and repeat this process for every year. The last step is to save this list with saveGIF() to create an animated gif. Note that we can use the same procedure but replace the bar plot with other types of plots such as pie charts.
counties_list<- as.data.frame(unique(counties$id)) #list of counties
all_list<- list()
for (y in 1:nrow(years)){ #loop through all of the years
nass_oneyear<- subset(yield,yield$Year==years[y,]) #extract one year
barplot_list <-
##create bar plot for each county
lapply(1:length(shp_counties$JURISDIC_2), function(i) {
#extract one county
nass_oneyear_onecounty<- subset(nass_oneyear,nass_oneyear$County==counties_list[i,])
# As I mentioned, for each county and year the number of crop types might be different. So, I need to customize the color for each sub-dataset using the manual color ramp that I previously defined for each item.
nass_oneyear_onecounty$itemcolor<- "NA"
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="CORN, GRAIN"]<- colors[1]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="SW, NON-IRRIGATED"]<- colors[2]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="SW, NON-IRRIGATED, CC"]<- colors[3]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="SW, NON-IRRIGATED, SF"]<- colors[4]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="WW, NON-IRRIGATED"]<- colors[5]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="WW, NON-IRRIGATED, CC"]<- colors[6]
nass_oneyear_onecounty$itemcolor[nass_oneyear_onecounty$Data_Item=="WW, NON-IRRIGATED, SF"]<- colors[7]
plots_comp <- ggplotGrob(
ggplot(nass_oneyear_onecounty,aes(x=Group,y=Value,group=(itemcolor),fill=factor(Data_Item)))+
scale_fill_manual(values=nass_oneyear_onecounty$itemcolor)+
geom_bar(stat='identity', color = "black") +
labs(x = NULL, y = NULL) +
theme(legend.position = "none", rect = element_blank(), axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x= element_blank(),
axis.ticks = element_blank(),
axis.text.y=element_text(size=14,face="bold",colour="black")) + coord_cartesian(expand=FALSE)
)})
barplots_list <- lapply(1:length(shp_counties), function(i)
annotation_custom(barplot_list[[i]], xmin = counties_centroids[i,1]- 28000, ymin = counties_centroids[i,2]- 28000, xmax = counties_centroids[i,1]+ 28000, ymax = counties_centroids[i,2]+ 28000))
# xmin, ymin, xmax and ymax can be used to adjust are horizontal and vertical location of the bar plots
all_list[[y]] <- list(tmp_map_final + barplots_list)
}
library(animation)
saveGIF(
{lapply(all_list, print)}
, "(your directory)/final.gif", interval = 2, ani.width = 1600, ani.height = 1200)
If we want to have labels for our bar plots (instead of having yield values at the y-axis), we may want to show the yield value corresponding to each item in a group. In this case I can use the below command lines within a loop before we create a graphical objects of ggplot (before plots_comp <- ggplotGrob(….) ), to add a label column that shows the cumulative yield value in each group.
nass_oneyear_onecounty <- nass_oneyear_onecounty %>%
group_by(Group) %>%
mutate(label_y = cumsum(Value)-10) #I subtracted 10 from these cumulative values just to print them inside of the bar plot sections.
Or we can just have one label that shows the total yield in each group:
nass_oneyear_onecounty <- nass_oneyear_onecounty %>%
group_by(Group) %>%
mutate(total = sum(Value))
Then we can add these labels to the bar plots by adding the geom_text() function in the ggplot section within ggplotGrob(), and specifying the column of interest:
geom_text(aes(y = label_y, label = round(Value)),colour = "white",size=3,fontface='bold')
Instead of manually adjusting the position the of label (such as in the first example), “vjust” can be added to the geom_text() for modifying text alignment:
geom_text(aes(y = total, label = round(total)),colour = "black",size=3,fontface='bold',vjust = -0.35)
Water Programming: A Collaborative Research Blog 