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

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

ABCDEFGHIJ0123456789

tinh <dbl+lbl>	huyen <dbl>	xa <dbl>	diaban <dbl>	hoso <dbl>	thunhap <dbl>	thubq <dbl>	tongthu <dbl>	chisxkd <dbl>	chikhac <dbl>
1	1	1	37	1	56400	14100	487800	684400	3500
1	1	1	37	2	31450	7862	377400	54000	0
1	1	1	37	4	1288400	9329	726700	0	250700
1	1	1	37	5	260700	65175	0	291000	525700
1	1	1	37	6	34297	17148	411570	579000	144000
1	1	1	37	7	23430	11715	360000	75000	23430

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)

ABCDEFGHIJ0123456789

tinh <dbl+lbl>	qt25 <dbl>	qt50 <dbl>	qt75 <dbl>	ISO3166_2_CODE <chr>
1	16916.00	107800.0	355201.5	VN-1
2	28798.00	55378.0	112766.0	VN-2
4	24427.75	44290.0	108741.5	VN-4
6	24014.75	47402.5	109239.0	VN-6
8	18663.50	63722.5	159498.0	VN-8
10	27000.00	67802.5	166745.0	VN-10
11	24916.75	47629.5	107230.2	VN-11
12	25974.25	61515.0	127302.8	VN-12
14	22378.75	54388.0	124274.0	VN-14
15	19631.25	56066.5	132581.8	VN-15

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