library(ggplot2)
library(sf)
library(dplyr)
library(haven)
library(stringi)
library(stringr)
library(rmarkdown)

ho3 <- read_dta("HO3.dta")
head(ho3) %>% 
  paged_table()

ho3 %>% 
  group_by(tinh) %>% 
  summarise(qt25 = quantile(thunhap, 0.25), 
            qt50 = quantile(thunhap, 0.50), 
            qt75 = quantile(thunhap, 0.75)) -> aggr_household_income

aggr_household_income <- aggr_household_income %>%
  mutate(ISO3166_2_CODE = paste("VN-", aggr_household_income$tinh, sep = ""))

paged_table(aggr_household_income)

vietnam <- read_sf("vietnam.json")

ggplot(vietnam) +
  geom_sf(fill = "white", color = "black", linewidth = 0.3) +
  theme_void()

vietnam_merged <- vietnam %>% 
  left_join(aggr_household_income, by="ISO3166_2_CODE")

head(vietnam_merged)

Simple feature collection with 6 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 104.7765 ymin: 9.021371 xmax: 107.0296 ymax: 22.74039
CRS:           NA
# A tibble: 6 × 9
  id    ISO3166_2_CODE Name_EN Name_VI                  geometry tinh       qt25
  <chr> <chr>          <chr>   <chr>              <MULTIPOLYGON> <dbl+lb>  <dbl>
1 <NA>  VN-89          An Gia… An Gia… (((105.1871 10.91317, 10… 89 [Tỉn… 20368.
2 <NA>  VN-24          Bac Gi… Bắc Gi… (((106.183 21.6053, 106.… 24 [Tỉn… 16226.
3 <NA>  VN-6           Bac Kan Bắc Kạn (((106.183 22.3883, 106.…  6 [Tỉn… 24015.
4 <NA>  VN-95          Bac Li… Bạc Li… (((105.3723 9.596911, 10… 95 [Tỉn… 40330.
5 <NA>  VN-27          Bac Ni… Bắc Ni… (((106.2974 21.12275, 10… 27 [Tỉn… 18187 
6 <NA>  VN-83          Ben Tre Bến Tre (((106.7725 10.1915, 106… 83 [Tỉn… 18870 
# ℹ 2 more variables: qt50 <dbl>, qt75 <dbl>

map <- ggplot(vietnam_merged) +
  geom_sf(aes(fill = qt50,
              text = paste("City:", Name_EN, "\n", "Avg. Household Income", qt50)),
          colour = NA) +
  labs(title = "Household Income in Vietnam, 2020",
       subtitle = "Data: Vietnam Household Living Standard Survey",
       fill = "millions VND") +
  theme_void() +
  theme(legend.position = "top",
        legend.margin = margin(t = 12, r = 0, b = 0, l = 0, unit = "pt"),
        legend.justification = c("left", "top")) +
  guides(fill = guide_legend(nrow = 1)) +
  theme(legend.key.height = unit(0.5, "cm")) + 
  theme(legend.key.width = unit(1.2, "cm"))

map

Acknowledgements

https://r-graph-gallery.com/327-chloropleth-map-from-geojson-with-ggplot2.html