Last active
January 8, 2021 12:28
-
-
Save anjesh/64cd61e3f94352ce5249ef368082df0c to your computer and use it in GitHub Desktop.
Step-by-step plotting choropleth map of Nepal
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(rgdal) | |
library(ggplot2) | |
library(dplyr) | |
# clone NepalMaps from https://github.com/anjesh/NepalMaps | |
# read shapefile | |
nepal_shp <- readOGR(dsn="NepalMaps/baselayers/NPL_adm", layer="NPL_adm3", stringsAsFactors = FALSE) | |
# fortify shapefile data to data frame | |
shp_df <- fortify(nepal_shp, region = "NAME_3") | |
# update the HPI data in districts.csv as per http://data.opennepal.net/content/human-poverty-index-value-districts-2011 | |
hpi.csvdata <- read.csv("https://github.com/opennepal/odp-poverty/raw/master/Human%20Poverty%20Index%20Value%20by%20Districts%20(2011)/data.csv", stringsAsFactors = FALSE) | |
hpi.data <- hpi.csvdata %>% | |
filter(Sub.Group == "HPI") %>% | |
select(id = District, HPI = Value) | |
# identify the mismatched districts by uncommenting 2 lines below | |
# unique(hpi.data$id[!hpi.data$id %in% unique(nepal.adm3.shp.df$id)]) | |
# unique(unique(nepal.adm3.shp.df$id)[!unique(nepal.adm3.shp.df$id) %in% hpi.data$id]) | |
# fix the mismatched districts | |
hpi.data$id[hpi.data$id == "Darchaula"] <- "Darchula" | |
hpi.data$id[hpi.data$id == "Kavre"] <- "Kavrepalanchok" | |
# merge hpi data with shapefile dataframe | |
shp_df <- shp_df %>% | |
left_join(hpi.data, by="id") | |
# finding the centers of all the district polygons. | |
centroids <- rgeos::gCentroid(nepal_shp, byid = TRUE, id = unique(nepal_shp$NAME_3)) | |
centroids_df <- centroids %>% | |
as.data.frame() %>% | |
mutate(label = row.names(.)) | |
# showing district names for which HPI>40 | |
centroids.selected <- centroids_df[centroids_df$label %in% (hpi.data[hpi.data$HPI>40,]$id),] | |
# setting up bare theme | |
theme_bare <- theme( | |
axis.line = element_blank(), | |
axis.text.x = element_blank(), | |
axis.text.y = element_blank(), | |
axis.ticks = element_blank(), | |
axis.title.x = element_blank(), | |
axis.title.y = element_blank(), | |
legend.text=element_text(size=7), | |
legend.title=element_text(size=8), | |
panel.background = element_blank(), | |
panel.border = element_rect(colour = "gray", fill=NA, size=0.5) | |
) | |
# plot map | |
ggplot(data = shp_df,aes(x = long, y = lat, group = group)) + | |
geom_polygon(aes(fill = HPI), color = 'gray', size = 0.1) + | |
ggtitle("Human Poverty Index Map") + | |
guides(fill=guide_colorbar(title="HP Index")) + | |
scale_fill_gradient(high = "#e34a33", low = "#fee8c8", guide = "colorbar") + | |
coord_fixed(1.3) + | |
theme(legend.justification=c(0,-0.1), legend.position=c(0.05,0)) + | |
with(centroids.selected, annotate(geom="text", x = x, y = y, label=label, size=2)) + | |
theme_bare |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(rgdal) | |
library(ggplot2) | |
library(dplyr) | |
# clone NepalMaps from https://github.com/anjesh/NepalMaps | |
# read shapefile | |
nepal.adm3.shp <- readOGR(dsn="./NepalMaps/baselayers/NPL_adm", layer="NPL_adm3", stringsAsFactors = FALSE) | |
# fortify shapefile data to data frame | |
nepal.adm3.shp.df <- fortify(nepal.adm3.shp, region = "NAME_3") | |
# write districts to csv file | |
nepal.adm3.shp.df %>% | |
distinct(id) %>% | |
write.csv("districts.csv", row.names = FALSE) | |
# update the HPI data in districts.csv as per http://data.opennepal.net/content/human-poverty-index-value-districts-2011 | |
hpi.data <- read.csv("districts.csv") | |
colnames(hpi.data) <- c("id","HPI") | |
# merge hpi data with shapefile dataframe | |
nepal.adm3.shp.df <- merge(nepal.adm3.shp.df, hpi.data, by ="id") | |
# finding the centers of all the district polygons. Thanks http://stackoverflow.com/questions/28962453/how-can-i-add-labels-to-a-choropleth-map-created-using-ggplot2 | |
centroids <- setNames(do.call("rbind.data.frame", by(nepal.adm3.shp.df, nepal.adm3.shp.df$group, function(x) {Polygon(x[c('long', 'lat')])@labpt})), c('long', 'lat')) | |
centroids$label <- nepal.adm3.shp.df$id[match(rownames(centroids), nepal.adm3.shp.df$group)] | |
# showing district names for which HPI>40 | |
centroids.selected <- centroids[centroids$label %in% (hpi.data[hpi.data$HPI>40,]$id),] | |
# setting up bare theme | |
theme_bare <- theme( | |
axis.line = element_blank(), | |
axis.text.x = element_blank(), | |
axis.text.y = element_blank(), | |
axis.ticks = element_blank(), | |
axis.title.x = element_blank(), | |
axis.title.y = element_blank(), | |
legend.text=element_text(size=7), | |
legend.title=element_text(size=8), | |
panel.background = element_blank(), | |
panel.border = element_rect(colour = "gray", fill=NA, size=0.5) | |
) | |
# plot map | |
ggplot(data = nepal.adm3.shp.df,aes(x = long, y = lat, group = group)) + | |
geom_polygon(aes(fill = HPI), color = 'gray', size = 0.1) + | |
ggtitle("Human Poverty Index Map") + | |
guides(fill=guide_colorbar(title="HP Index")) + | |
scale_fill_gradient(high = "#e34a33", low = "#fee8c8", guide = "colorbar") + | |
coord_fixed(1.3) + | |
theme(legend.justification=c(0,-0.1), legend.position=c(0.05,0)) + | |
with(centroids.selected, annotate(geom="text", x = long, y = lat, label=label, size=2)) + | |
theme_bare |
Thankx a lot Anjesh ! I am a beginner and am self-learning and this helps me a lot.
Thanks so much Anjesh. I worked through the steps and I got it right to a point where I got stocked.. I am working on Africa and I am only selecting 20 countries out of Africa.. I used a percapita GDP of these countries for a particular year and I set the other countries which is not among my study to zero.. How ever when I tried to plot the cloroplethr with label name.. It couldn't imprint names on the countries .. Showing geom_text deleted 768 missing and it returned to the normal map without labels..
You might need to install extra libraries as well.
install.packages('maptools')
install.packages('rgeos')
install.packages('rgdal')
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I have simplified the script https://gist.github.com/anjesh/64cd61e3f94352ce5249ef368082df0c#file-map-choropleth-nepal-updated-r You should be able to make changes to the code to make the centroids work as per your need