#setting working directory

setwd("C:\\Users\\USER\\Desktop\\Viec thay Ngoc\\ASEAN EM 2019, Bangkok\\Report\\Data of General Statistic Center")

#loading packages

library(ggplot2)
library(maptools)
## Loading required package: sp
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
library(rgdal)
## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1

#loading Vietnam provincial data

library(raster)
load(file="vietnam.RData")

#convert to dataframe

library(tidyverse)
## -- Attaching packages -------------------- tidyverse 1.2.1 --
## v tibble  2.1.3     v purrr   0.3.2
## v tidyr   0.8.3     v dplyr   0.8.3
## v readr   1.3.1     v stringr 1.4.0
## v tibble  2.1.3     v forcats 0.4.0
## -- Conflicts ----------------------- tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
if (!require(gpclib)) install.packages("gpclib", type="source")
## Loading required package: gpclib
## General Polygon Clipper Library for R (version 1.5-5)
##  Type 'class ? gpc.poly' for help
gpclibPermit()
## Warning in gpclibPermit(): support for gpclib will be withdrawn from
## maptools at the next major release
## [1] TRUE
vn_df <- vietnam %>% fortify(region = "GID_1")
head(vn_df)
##       long      lat order  hole piece      id     group
## 1 105.3745 10.24604     1 FALSE     1 VNM.1_1 VNM.1_1.1
## 2 105.3362 10.23442     2 FALSE     1 VNM.1_1 VNM.1_1.1
## 3 105.3154 10.26842     3 FALSE     1 VNM.1_1 VNM.1_1.1
## 4 105.3120 10.27393     4 FALSE     1 VNM.1_1 VNM.1_1.1
## 5 105.3065 10.26894     5 FALSE     1 VNM.1_1 VNM.1_1.1
## 6 105.3013 10.27698     6 FALSE     1 VNM.1_1 VNM.1_1.1

#loading physician dataset

pd <- read.csv("C:\\Users\\USER\\Desktop\\Viec thay Ngoc\\ASEAN EM 2019, Bangkok\\Report\\Data of General Statistic Center\\healthcare worker_2017.csv",
header = TRUE)

#merge datasets

vn_pd <- merge(vn_df, pd, by = "id", all.x = TRUE)

#plot choropleth map

p1 <- vn_pd %>% 
  ggplot(aes(x = long, y = lat, group = group))+
  geom_polygon(aes(fill = bs), color = "grey30") +  
  labs(x = NULL, y = NULL) 

p1

# upgrade the map
library(viridis)
## Loading required package: viridisLite
library(viridisLite)
p2 <- p1 + scale_fill_viridis(direction = -1, option = "D", "Physician density")
p2

p3 <- p2 + scale_fill_viridis(option = "C", direction = -1, "Physicians \nper 10,000 inhabitants") + 
  labs(title = "Physician Density in Vietnamese \nPublic Healthcare Facilities (lastest available data)", 
       caption = "Data Source: General Statistics Office of Vietnam") + 
  theme(axis.line = element_blank(), 
        axis.text = element_blank(),
        axis.ticks = element_blank(), 
        axis.title = element_blank()) +
  theme_bw() +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
p3

centroids <- setNames(do.call("rbind.data.frame", 
by(vn_pd, vn_pd$tinh, function(x)
{Polygon(x[c('long', 'lat')])@labpt})), c('long', 'lat'))
centroids$label <- vn_pd$tinh[match(rownames(centroids), vn_pd$group)]


p4 <- p3 +
with(centroids, annotate(geom="text", x = long, y = lat, label=label, size=2))

p4
## Warning: Removed 63 rows containing missing values (geom_text).