Visualizing Data in Map with R

R is a powerful software environment for graphics and data visualization. Here I have used R to produce data visualization in the Map of Nepal.  I have plotted the poverty rate by districts in the map as shown below: 

 

Here is the video for creating the data visualization explained below.

Here is the link to YouTube : https://youtu.be/f26U2kwAWkQ

The datafiles and R script can be downloaded from here. 

The first step here is to install the necessary R packages. This can be done by running the following codes in R/R studio:

install.packages("cartography")

install.packages("sf")

install.packages("tidyverse")

Now let us bring the packages into the memory of R by running :

library(cartography)
library(sf)
library(tidyverse)

The second step is to set the working directory. Mine is 'map' folder in my desktop.
setwd("C:/Users/siddhabhatta/Desktop/map")

The third step is to download the  shape files for Nepal and extract them in the working directory.  We can simply find the shape files by googling. 

 # I found them for Nepal at https://codefornepal.carto.com/tables/shape_files_of_districts_in_nepal/public

Alternatively, the files can be downloaded from here.

# Simply download the files and put them in the working directory.  There may be more than one files with the format .cpg, .dbf, .shp, .prj and .shx. We need the files with .dbf and .shp extension. 

The next step is to import the data and save it as some object 'data' in my case.

 data<-st_read("shape_files_of_districts_in_nepal.shp")

The above data contains all the information for drawing a map of Nepal.

We next need the data that we want to visualize in the map. The next two line imports the data containing district wise poverty rates in R.  The poverty data file can be downloaded by clicking here.

library(haven)
poverty <- read_dta("poverty.dta")
View(poverty)
poverty$pov<-poverty$poverty*100 # converts the poverty rate into percentages
View(poverty)

Now, we combine the poverty data with the main data set that we have imported from the shape file of Nepal. This can be done by running the following command. The combined data is saved as 'mapdata'
mapdata<-merge(data, poverty, by="dist_name")
View(mapdata)
names(mapdata)

Now, the map can be produced with the following command. 

 plot(st_geometry(mapdata)) # plots the map of Nepal
choroLayer(x=mapdata, var="pov", method="quantile", nclass=5) # fills the districts with color as per the poverty rate.
layoutLayer(title="Poverty by Districts: Nepal ", tabtitle = TRUE, frame=TRUE, scale=6)

Here are minor changes in the plot. 

 plot(st_geometry(mapdata))
choroLayer(x=mapdata, var="pov", method="quantile", nclass=8, legend.title.txt = "Poverty")
layoutLayer(title="Poverty in Nepal by Districts", tabtitle = TRUE, frame=TRUE, scale=6)

#using ggplot package, the graph can be enhanced.

ggplot(data=mapdata)+geom_sf(aes(fill=pov), color="white")+
  scale_fill_viridis_c(option = "viridis", trans = "sqrt")+
  xlab("Longitute")+ ylab("latitude")+
  ggtitle("Poverty Rate of Nepal by Districts")

#To plot the labels, we extract the centroids and save them as X and Y by the following command.

 points<-cbind(mapdata, st_coordinates(st_centroid(mapdata$geometry)))

library(ggthemes) # used for changes background theme, title theme and others
ggplot(data = points) +
  geom_sf(aes(fill=pov), color="black", size=0.2) +
  scale_fill_viridis_c(option = "viridis", trans = "sqrt")+
  geom_text(data= points,aes(x=X, y=Y, label=paste(dist_name)),
            color = "darkblue", size=2.5, fontface = "italic", angle=0, vjust=-1, check_overlap = FALSE) +
  geom_text(data= points,aes(x=X, y=Y, label=paste(pov)),
            color = "white", size=2.0, fontface = "bold", angle=0, vjust=+1, check_overlap = FALSE) +
  ggtitle("Poverty Rate of Nepal by Districts")+xlab("Longitude")+ylab("Latitude")+

  theme(
    panel.background = element_rect(fill='lightblue', size=0.5,
                                        linetype=0, colour='lightblue'),
  plot.title = element_text(color="red", size=16, hjust=0.5, face="bold.italic"),
  axis.title.x = element_text(color="blue", size=10, face="bold"),
  axis.title.y = element_text(color="red", size=10, face="bold")
)

We can save the plot by using ggsave command.  

ggsave("map.png", width=6, height=6, dpi='screen')

ggsave("map.pdf", width=6, height=6, dpi='screen')

...............................................................................................................................................................

 If you are interested in creating such map with stata, here is the link.