#Call package
pacman::p_load(sf,
                rio,
               here,
               janitor,
               tidyverse,
               dplyr,
               magrittr,
               ggplot2,
               purrr,
               lubridate,
               mice,
               plotly)

Import data:

I download VietNam shapefile in level 1 (adm1) from the page HDX. If you want to analyst deeply, I recommend this post ThangLeQuoc. You can combine R with SQL to implement your study.

Simple plot of Vietnam population:

#Import Viet Nam shapefile in level 1 (include 63 provinces):
getwd()
## [1] "C:/Users/locca/Documents/Xuân Lộc/R/Project R/Amazon"
#Take new workdirectory in clipboard

vn_adm1<-sf::st_read(r"(C:\Users\locca\Documents\Xuân Lộc\R\Data R\Viet Nam shapefile\vnm_admbnda_adm1_gov_20201027.shp)")
## Reading layer `vnm_admbnda_adm1_gov_20201027' from data source 
##   `C:\Users\locca\Documents\Xuân Lộc\R\Data R\Viet Nam shapefile\vnm_admbnda_adm1_gov_20201027.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 63 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.144 ymin: 7.180931 xmax: 117.8355 ymax: 23.39221
## Geodetic CRS:  WGS 84
vn_pop<-import("C:/Users/locca/Downloads/vnm_admpop_adm1_2022_v2.csv")


#Join data:
pop<-right_join(vn_adm1,
                vn_pop %>% select(ADM1_NAME,F_TL,M_TL,T_TL),
                by = c("ADM1_VI" = "ADM1_NAME")) 
              
           
#Plot map:
ggplot(data = pop)+
  geom_sf(aes(fill=T_TL))

Interactive Viet Nam with {leaflet}:

#Use leaflet to plot map:
library(leaflet)
pal <- colorBin(
  palette = "viridis", domain = pop$T_TL,
  bins = 10
)

labels<- paste0("<strong> Province: </strong> ",
               pop$ADM1_VI, "<br/> ",
               "<strong> Total population: </strong> ",
               pop$T_TL, "<br/> ") %>% 
         lapply(htmltools::HTML)

leaflet(pop) %>% 
  addPolygons(
     stroke = FALSE,
     fillColor = ~pal(T_TL),
     fillOpacity = 0.7,
     smoothFactor = 0.5,
     label = ~labels,
     highlight = highlightOptions(
      color = "black",
      bringToFront = TRUE
    )
  ) %>% 
  addTiles() %>% 
  addLegend("bottomright",  
            pal = pal,
            values =~T_TL,
            title = "Population per province in year 2022")

Plot the road map for three provinces in Viet Nam

You can download file from Giao thông Việt Nam

#Add road map:
road_vn<-st_read("C:\\Users\\locca\\Downloads\\giaothong\\Giao_Thong_WGS84.shp")
## Reading layer `Giao_Thong_WGS84' from data source 
##   `C:\Users\locca\Downloads\giaothong\Giao_Thong_WGS84.shp' using driver `ESRI Shapefile'
## Simple feature collection with 458 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 102.9073 ymin: 8.641914 xmax: 109.3985 ymax: 23.2791
## Geodetic CRS:  WGS 84
#Select 3 provinces to analysis:
province <- vn_adm1 %>% filter( ADM1_VI %in% c("Phú Yên", "Bình Định", "Khánh Hòa"))

#Filter a lot roads intersects the geometry of 3 provinces:
road_intersect<-st_intersection(province,road_vn)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
ggplot()+
  geom_sf(data = province)+
  geom_sf(data = road_intersect, color = "blue")