In a previous article on ‘Forest and Tree Statistics in Kenya’ (https://medium.com/@wanjirumaggie45/forest-and-tree-statistics-in-kenya-using-r-6cfa9f67ea9b), we looked at Kenya’s current forest cover which stands at 7.6 %. The analysis was based on this data set: https://rainforests.mongabay.com/deforestation/archive/Kenya.htm. The results were okay but they got me thinking about delving deeper into the statistics such as looking at the data per county. Luckily, the same data source was neat for this case study as well.
Making beautiful data visualizations is something I enjoy doing. So I challenged myself to do an interactive map using the ‘plotly’ package in R as we will see below.
1. Load Packages
The rgdal package provides bindings for the ‘Geospatial’ Data Abstraction Library/ Reading shapefiles function into R (readOGR)(https://www.rdocumentation.org/packages/rgdal).
library(rgdal) #reading shapefiles
library(plotly) #interactive visualization
2. Loading the Data
Notice the function setNames used below. It is a convenience function that sets the names on an object and returns the object.
forest_county <- read.csv("C:/Users/MARGRET/Desktop/Blog/Data/forest_county.csv")forest_county<-setNames(forest_county, c("IEBCCty_No",
"COORD_X", "COORD_Y", "Percentage(%) Loss since (2000)" ))
3. Loading the Shapefiles
The fortify function converts a curves and points object to a data frame for ggplot2. *Warning* However, rather than using this function, the broom package is recommended, which implements a much wider range of methods.
fortify may be deprecated in the future. For now, let’s just use it since it’s still working well.
kenyaShp1<- readOGR("C:/Users/MARGRET/Desktop/shapefilesGEO")##Converting the Shp to a data framekenyaShpdf <- fortify(kenyaShp1)
4. Merging the Shapefiles and the data frame
By default, the data frames are merged on the columns with names they both have, but separate specifications of the columns can be given by
popDF <- merge(kenyaShpdf, forest_county, by.x = "id", by.y = "IEBCCty_No" )
5. And Finally the Interactive Map!
Map <- ggplot(popDF, aes(long, lat, group = group, fill = `Percentage(%) Loss since (2000)`,
label = County)) +
labs(fill = "% Loss since 2000") +
ggtitle("Tree Loss Since 2000 (2000 - 2018)\nKenyan Counties")+
theme_minimal(base_size = 12) +
plot.caption =element_text(size=10, hjust=0.5, face="italic", color="blue"))+
labs(caption = "Data Source: rainforests.mongabay.com")Map <- Map + scale_fill_gradient(low = "gray", high = "red")
Believe it or not the code below now produces the interactive map
Now, what if we wanted to compare tree cover(%) and tree percentage loss just to get an overview of the transition? Notice how some counties like Mandera, Wajir and Marsabit barely have any trees.
Nyeri has the highest percentage of tree cover at 53% followed by Nyamira at 50%.
Thanks for reading! If you would like to try this code and have no access to shapefiles, feel free to ask and I will gladly share.