Ex 1: Mark the corresponding stations on a map.

library(sf)
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(ggplot2)
# Read the Vietnam map
vietnam_map <- st_read("mapvietnam")
## Reading layer `mapvietnam' from data source 
##   `D:\USTH\Year_3\Analysis of Spatial & Temporal Data\Spatial Mapping\mapvietnam' 
##   using driver `GeoJSON'
## Simple feature collection with 63 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.1421 ymin: 6.953306 xmax: 116.9473 ymax: 23.3939
## Geodetic CRS:  WGS 84
# Load the station data
stations <- read.table("station.list.txt", header=TRUE, stringsAsFactors=FALSE)

# Convert to spatial data frame
stations_sf <- st_as_sf(stations, coords = c("LONG", "LAT"), crs = 4326)

# Plot the map
ggplot() +
  geom_sf(data = vietnam_map, fill = "lightgray", color = "black") +  # Vietnam boundary
  geom_sf(data = stations_sf, color = "red", size = 1.5) +  # Stations
  theme_minimal() +
  labs(title = "Meteorological Stations in Vietnam",
       x = "Longitude",
       y = "Latitude")

# Ex 2: Interactive map

library(leaflet)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# provinces as red dots 
# Load station data
stations <- read.table("station.list.txt", header=TRUE, stringsAsFactors=FALSE)

# Create leaflet map
leaflet(data = stations) %>%
  addProviderTiles(providers$OpenStreetMap) %>%  # Change map provider if needed
  addCircleMarkers(
    lng = ~LONG, lat = ~LAT, 
    popup = ~STATION_NAME,
    radius = 3, color = "red", stroke = FALSE, fillOpacity = 0.8
  ) %>%
  setView(lng = mean(stations$LONG), lat = mean(stations$LAT), zoom = 6) %>%
  addMiniMap()

Ex 3: PCI 2023 in Vietnam normal plot

library(ggplot2)
library(sf)
library(readxl)
library(dplyr)
library(stringr)
library(tidyr)
# Load Vietnam map
vietnam_map <- st_read("mapvietnam")
## Reading layer `mapvietnam' from data source 
##   `D:\USTH\Year_3\Analysis of Spatial & Temporal Data\Spatial Mapping\mapvietnam' 
##   using driver `GeoJSON'
## Simple feature collection with 63 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 102.1421 ymin: 6.953306 xmax: 116.9473 ymax: 23.3939
## Geodetic CRS:  WGS 84
# Load PCI data
pci_data <- read_excel("_data_pci_file_general_file_1715183666-Appendix_1_Table_Indicator_PCI2023_VN_web.xlsx", sheet = "Tổng hợp")

# Rename the first column explicitly to "Province"
colnames(pci_data)[1] <- "Province"

# Standardize column names (replace spaces with underscores)
colnames(pci_data) <- gsub(" ", "_", colnames(pci_data))

# Standardize province names in both datasets
pci_data <- pci_data %>%
  mutate(Province = str_trim(tolower(Province)))  # Clean province names

vietnam_map <- vietnam_map %>%
  mutate(Name = str_trim(tolower(Name)))  # Ensure Name column is formatted correctly

vietnam_map <- vietnam_map %>%
  left_join(pci_data, by = c("Name" = "Province"))

# Extract the 10 PCI factors (columns 2 to 11)
pci_factors <- colnames(pci_data)[2:11]
plot_pci_map <- function(factor_name, title_text) {
  ggplot() +
    geom_sf(data = vietnam_map, aes_string(fill = factor_name), color = "black", size = 0.2) + 
    scale_fill_viridis_c(option = "plasma", name = "PCI Score", guide = guide_colorbar(barwidth = 8, barheight = 0.5)) +
    theme_minimal() +
    theme(legend.position = "bottom", legend.title = element_text(size = 10),
          plot.title = element_text(size = 14, face = "bold"),
          axis.text = element_blank(), axis.ticks = element_blank(),
          panel.grid = element_blank()) +
    labs(title = paste("Vietnam PCI 2023:", title_text))
}

# Generate 10 separate plots
plot1 <- plot_pci_map("CSTP_1_Mark_entry", "Market Entry (CSTP 1)")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot2 <- plot_pci_map("CSTP_2_Land_access", "Land Access (CSTP 2)")
plot3 <- plot_pci_map("CSTP_3_Transparency", "Transparency (CSTP 3)")
plot4 <- plot_pci_map("CSTP_4_Time_costs", "Time Costs (CSTP 4)")
plot5 <- plot_pci_map("CSTP_5_Informal_Costs", "Informal Costs (CSTP 5)")
plot6 <- plot_pci_map("CSTP_6_Fair_competition", "Fair Competition (CSTP 6)")
plot7 <- plot_pci_map("CSTP_7_Government_Proactivity_and_Initiative", "Government Proactivity (CSTP 7)")
plot8 <- plot_pci_map("CSTP_8_Business_support_policies", "Business Support Policies (CSTP 8)")
plot9 <- plot_pci_map("CSTP_9_Labor_training", "Labor Training (CSTP 9)")
plot10 <- plot_pci_map("CSTP_10_Legal_Institutions_and_Pubic_order_and_security", "Legal Institutions & Public Order (CSTP 10)")

# Display each plot separately
print(plot1)

print(plot2)

print(plot3)

print(plot4)

print(plot5)

print(plot6)

print(plot7)

print(plot8)

print(plot9)

print(plot10)