Các bạn sử dụng package FAOSTAT để có thể download dataset về nông nghiệp về máy tính.

Hướng dẫn chi tiết: https://cran.r-project.org/web/packages/FAOSTAT/index.html

BƯỚC 1: DOWNLOAD DATASET

Thay vì download thủ công từ website của FAO https://www.fao.org/faostat/en/#data/QCL thì ta sẽ download trực tiếp bộ dataset nông nghiệp (mã QCL) về máy tính và dùng R phân tích dữ liệu.

library(FAOSTAT)
data_folder <- "data_raw"
dir.create(data_folder) ## tạo folder lưu file trên máy tính
fao_metadata <- FAOsearch()
FAOsearch(dataset = "crop", full = FALSE)

## dòng lệnh này sẽ thực thi truy vấn download full dataset nông nghiệp về máy tính
crop_production <- get_faostat_bulk(code = "QCL", data_folder = data_folder)

## bạn save object này ra thành 1 file riêng biệt để thuận tiện sử dụng sau này
saveRDS(crop_production, "data_raw/crop_production_e_all_data.rds")

BƯỚC 2: LÀM SẠCH DỮ LIỆU

Thực tế dataset của FAO đã được làm sạch rồi. Tuy nhiên để thuận tiện minh họa cho ví dụ này, mình sẽ lọc ra dữ liệu của lúa gạo từ số liệu chính thức của Việt Nam, flag A trong dataset.

library(dplyr)
library(tidyr)
library(DT)

crop_production <- readRDS("data_raw/crop_production_e_all_data.rds")
dim(crop_production)
## [1] 3761216      13
str(crop_production)
## 'data.frame':    3761216 obs. of  13 variables:
##  $ area_code      : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ area_code__m49_: chr  "'004" "'004" "'004" "'004" ...
##  $ area           : chr  "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
##  $ item_code      : int  221 221 221 221 221 221 221 221 221 221 ...
##  $ item_code__cpc_: chr  "'01371" "'01371" "'01371" "'01371" ...
##  $ item           : chr  "Almonds, in shell" "Almonds, in shell" "Almonds, in shell" "Almonds, in shell" ...
##  $ element_code   : int  5312 5312 5312 5312 5312 5312 5312 5312 5312 5312 ...
##  $ element        : chr  "area_harvested" "area_harvested" "area_harvested" "area_harvested" ...
##  $ year_code      : int  1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 ...
##  $ year           : int  1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 ...
##  $ unit           : chr  "ha" "ha" "ha" "ha" ...
##  $ value          : num  0 5900 6000 6000 6000 5800 5800 5800 5700 5700 ...
##  $ flag           : chr  "E" "E" "E" "E" ...
## tách riêng dataset của Việt Nam
vietnam_data <- subset(crop_production, area == "Viet Nam")
dim(vietnam_data)
## [1] 15070    13
## subset dữ liệu Flag A từ nguồn chính thức do Việt Nam công bố
vietnam_data_a <- subset(vietnam_data, flag == "A")
dim(vietnam_data_a)
## [1] 4422   13
## chọn những cột cần thiết
vietnam_crop <- vietnam_data_a[, c("item", "element", "year", "unit", "value")]

## sắp xếp theo năm và thông tin nông sản
vietnam <- vietnam_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
str(vietnam)
## 'data.frame':    4422 obs. of  5 variables:
##  $ item   : chr  "Tea leaves" "Tea leaves" "Tea leaves" "Swine / pigs" ...
##  $ element: chr  "yield" "production" "area_harvested" "stocks" ...
##  $ year   : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ unit   : chr  "hg/ha" "tonnes" "ha" "Head" ...
##  $ value  : num  97545 1073000 110000 23533400 125413 ...
## chúng ta có clean dataset để đưa vào phân tích giai đoạn sau
datatable(vietnam, options = list(pageLength = 30))

Đặc trưng bộ dữ liệu này bao gồm thông tin nông sản vật nuôi và cây trồng, gần như toàn diện về sản lượng nông nghiệp ở Việt Nam từ 1961 đến 2021. Các bạn có thể sử dụng công cụ vẽ đồ thị trên website của FAOSTAT hoặc dùng R để vẽ các loại đồ thị bạn quan tâm, trả lời cho câu hỏi nghiên cứu về tình hình nông nghiệp ở Việt Nam.

https://www.fao.org/faostat/en/#compare

Trong bài hướng dẫn này, mình sẽ sử dụng dữ liệu lúa gạo Rice trong dataset để vẽ đồ thị bubble chart nhằm đánh giá tổng quát tương quan giữa năng suất và diện tích canh tác theo thời gian ở Việt Nam.

BƯỚC 3: LỌC DỮ LIỆU VỀ LÚA GẠO ĐƯA VÀO ĐỒ THỊ

rice <- subset(vietnam, item == "Rice")
rice_ok <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)

head(rice_ok)
##     item        element year    value
## 22  Rice     production 2021 43852729
## 23  Rice area_harvested 2021  7219797
## 132 Rice     production 2020 42765000
## 133 Rice area_harvested 2020  7222403
## 293 Rice     production 2019 43495488
## 294 Rice area_harvested 2019  7451544
## spread dataset để ggplot nhận diện dữ liệu
rice_final <- rice_ok %>% tidyr::spread(key = "element", value = "value")

head(rice_final)
##   item year area_harvested production
## 1 Rice 1961        4744000    8997400
## 2 Rice 1962        4888860    9747040
## 3 Rice 1963        4496520    9622670
## 4 Rice 1964        4987800    9697030
## 5 Rice 1965        4826300    9369700
## 6 Rice 1966        4681300    8463500
na.omit(rice_final) -> rice_final

BƯỚC 4: VẼ ĐỒ THỊ

library(ggplot2)
options(scipen = 99)

p <- ggplot(rice_final, aes(
  x = year, y = production / 1000000,
  size = area_harvested / 1000,
  fill = production / area_harvested
)) +
  geom_point(alpha = 0.6, shape = 21, color = "black") +
  scale_size(range = c(1, 9), name = "Diện tích canh tác (nghìn ha)", 
             breaks = round(seq(min(rice_final$area_harvested)/1000, max(rice_final$area_harvested)/1000, length.out = 6),
                            digits = -2)) + 
  scale_fill_viridis_c(
    name = "Năng suất (tấn/ha)",
    option = "D",
    guide = "colourbar",
    breaks = round(seq(min(rice_final$production / rice_final$area_harvested), 
                 max(rice_final$production / rice_final$area_harvested), length.out = 6), digits = 0)
  ) +
  guides(
    size = guide_legend(order = 1),
    fill = guide_colourbar(barwidth = 1, barheight = 8, order = 2)
  ) +
  scale_x_continuous(expand = c(0, 0), limits = c(1955, 2025)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1.15 * max(rice_final$production) / 1000000)) +
  
  ## customize theme
  
  theme_classic() +
  ylab("Sản lượng thu hoạch (triệu tấn)") +
  xlab("Năm") +
  
  ## x and y lab
  theme(axis.title.x = element_text(face = "bold", size = 10)) +
  
  theme(axis.title.y = element_text(
    face = "bold",
    margin = margin(t = 1, r = 10, b = 1, l = 10)
  )) +
  
  ## tick text
  
  theme(
    axis.text.x = element_text(
      face = "bold", color = "#993333",
      size = 10
    ),
    axis.text.y = element_text(
      face = "bold", color = "#993333",
      size = 10
    )
  ) +
  
  ## title
  
  labs(title = "Tình hình sản xuất lúa gạo ở Việt Nam từ 1960–2021",
              subtitle = "Nguồn: FAOSTAT [crop dataset, Vietnam, flag A]",
              caption = "Đồ thị: Duc Nguyen | tuhocr.com") +
  
  theme(
  plot.title = element_text(size = 12, 
                            face = "bold", 
                            hjust = 0, margin = margin(t = 10, r = 0, b = 5, l = 0)),
  plot.subtitle = element_text(size = 11, hjust = 0, margin = margin(t = 0, r = 0, b = 5, l = 0)),
  plot.caption = element_text(size = 9, hjust = 1, margin = margin(t = 0, r = 0, b = 5, l = 0))
) +
  
  ## legend

  theme(legend.position = "right", 
        legend.title = element_text(size = 10, face = "bold"),
        legend.background = element_blank(),
        legend.justification = c("right", "bottom"),
        # legend.box.just = "right",
        legend.margin = margin(t = 0, r = 0, b = 0, l = 5),
        legend.spacing.x = unit(0.5, "lines"),
        legend.spacing.y = unit(0.4, "lines")) +
  
  ## margin
  
  theme(
    panel.background = element_rect(fill = "white"),
    plot.margin = margin(t = 3, r = 14, b = 2, l = 2),
    plot.background = element_rect(fill = "grey90", colour = "black", linewidth = 1)
  ) +

  ## bổ sung thêm
# p + 

#   # annotate(geom = "text", x = 2018, y = 0.85*max(rice_final$production)/1000,
#   #          label = "Năm 2021 \nđạt 43,8 nghìn tấn",
#   #          color = "blue",
#   #          size = 3) +
  geom_hline(yintercept = 48, colour = "#148c34", lwd = 1.2, lty = 2) +
  annotate(
    geom = "text", x = 2000, y = 1.1 * max(rice_final$production) / 1000000,
    label = "Ngưỡng cực hạn 48 triệu tấn/năm",
    color = "red"
  ) +
#   coord_cartesian(clip = "off")
  geom_label(aes(x = 1975, y = 0.9 * max(rice_final$production)/1000000, 
                 label = paste("Năm", max(rice_final$year), "Việt Nam", "\nthu hoạch được", round(rice_final[rice_final$year == max(rice_final$year), ]$production/1000000, digits = 0), "triệu tấn gạo", "\ntrên diện tích canh tác", round(rice_final[rice_final$year == max(rice_final$year), ]$area_harvested/1000, digits = 0), "nghìn ha", "\nnăng suất là", round(rice_final[rice_final$year == max(rice_final$year), ]$production/rice_final[rice_final$year == max(rice_final$year), ]$area_harvested, digits = 2), "tấn/ha")),
             fill = "#f5ec42", label.size = 0.15, size = 3.5)

## gọi đồ thị
p

## add hình
library(grid)
library(png)
image <- readPNG("rice.png")
logor <- readPNG("logor.png")
grid.raster(image, x = 0.56, y = 0.24, width = 0.23)
grid.raster(logor, x = 0.9, y = 0.9, width = 0.1)

BƯỚC 5: SETUP DATASET VỀ LÚA GẠO Ở CÁC QUỐC GIA KHÁC

THAILAND-RICE DATASET

thailand_data <- subset(crop_production, area == "Thailand")
dim(thailand_data)
## [1] 20992    13
## subset dữ liệu Flag A từ nguồn chính thức do Thái Lan công bố
thailand_data_a <- subset(thailand_data, flag == "A")
dim(thailand_data_a)
## [1] 6441   13
## chọn những cột cần thiết
thailand_crop <- thailand_data_a[, c("item", "element", "year", "unit", "value")]

## sắp xếp theo năm và thông tin nông sản
thailand <- thailand_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
str(thailand)
## 'data.frame':    6441 obs. of  5 variables:
##  $ item   : chr  "Sugar cane" "Sugar cane" "Sugar cane" "Sugar Crops Primary" ...
##  $ element: chr  "yield" "production" "area_harvested" "yield" ...
##  $ year   : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ unit   : chr  "hg/ha" "tonnes" "ha" "hg/ha" ...
##  $ value  : num  443211 66278506 1495417 443211 66278506 ...
# ## chúng ta có clean dataset để đưa vào phân tích giai đoạn sau
# datatable(thailand, options = list(pageLength = 30))

## từ đoạn này trở đi thì dùng code cũ, vì object sẽ tự động pass vào.
rice <- subset(thailand, item == "Rice")
rice_ok <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)

head(rice_ok)
##     item        element year    value
## 8   Rice     production 2021 33582000
## 9   Rice area_harvested 2021 11244000
## 41  Rice     production 2020 30231025
## 42  Rice area_harvested 2020 10401653
## 110 Rice     production 2019 28617948
## 111 Rice area_harvested 2019  9812614
## spread dataset để ggplot nhận diện dữ liệu
rice_final <- rice_ok %>% tidyr::spread(key = "element", value = "value")

head(rice_final)
##   item year area_harvested production
## 1 Rice 1963             NA   12171000
## 2 Rice 1964             NA   11600000
## 3 Rice 1965             NA   11164000
## 4 Rice 1966        7353000   13500000
## 5 Rice 1967        6400000   11198000
## 6 Rice 1968        6940000   12410000
na.omit(rice_final) -> rice_final_thailand

CHINA-RICE DATASET

china_data <- subset(crop_production, area == "China, mainland")
dim(china_data)
## [1] 30038    13
## subset dữ liệu Flag A từ nguồn chính thức do China công bố
china_data_a <- subset(china_data, flag == "A")
dim(china_data_a)
## [1] 5835   13
## chọn những cột cần thiết
china_crop <- china_data_a[, c("item", "element", "year", "unit", "value")]

## sắp xếp theo năm và thông tin nông sản
china <- china_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
str(china)
## 'data.frame':    5835 obs. of  5 variables:
##  $ item   : chr  "Wheat" "Wheat" "Wheat" "Unmanufactured tobacco" ...
##  $ element: chr  "yield" "production" "area_harvested" "production" ...
##  $ year   : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ unit   : chr  "hg/ha" "tonnes" "ha" "tonnes" ...
##  $ value  : num  58106 136946000 23568400 2127600 449224200 ...
## từ đoạn này trở đi thì dùng code cũ, vì object sẽ tự động pass vào.
rice <- subset(china, item == "Rice")
rice_ok <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)

head(rice_ok)
##     item        element year     value
## 24  Rice     production 2021 212843000
## 25  Rice area_harvested 2021  29921200
## 86  Rice     production 2020 211860000
## 87  Rice area_harvested 2020  30080000
## 164 Rice     production 2019 209614000
## 165 Rice area_harvested 2019  29690000
## spread dataset để ggplot nhận diện dữ liệu
rice_final <- rice_ok %>% tidyr::spread(key = "element", value = "value")

head(rice_final)
##   item year area_harvested production
## 1 Rice 1961       26250000   53640000
## 2 Rice 1962       26908000   62985008
## 3 Rice 1963       27687008   73765008
## 4 Rice 1964       29577008   83000000
## 5 Rice 1965       29795008   87720000
## 6 Rice 1966       30498000   95390000
na.omit(rice_final) -> rice_final_china

INDIA-RICE DATASET

india_data <- subset(crop_production, area == "India")
dim(china_data)
## [1] 30038    13
## subset dữ liệu Flag A từ nguồn chính thức do China công bố
india_data_a <- subset(india_data, flag == "A")
dim(india_data_a)
## [1] 9846   13
## chọn những cột cần thiết
india_crop <- india_data_a[, c("item", "element", "year", "unit", "value")]

## sắp xếp theo năm và thông tin nông sản
india <- india_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
str(india)
## 'data.frame':    9846 obs. of  5 variables:
##  $ item   : chr  "Wheat" "Wheat" "Wheat" "Watermelons" ...
##  $ element: chr  "yield" "production" "area_harvested" "yield" ...
##  $ year   : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ unit   : chr  "hg/ha" "tonnes" "ha" "hg/ha" ...
##  $ value  : num  34669 109590000 31610000 273445 3254000 ...
## từ đoạn này trở đi thì dùng code cũ, vì object sẽ tự động pass vào.
rice <- subset(india, item == "Rice")
rice_ok <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)

head(rice_ok)
##     item        element year     value
## 29  Rice     production 2021 195425000
## 30  Rice area_harvested 2021  46379000
## 182 Rice     production 2020 186500000
## 183 Rice area_harvested 2020  45070000
## 382 Rice     production 2019 178305480
## 383 Rice area_harvested 2019  43662300
## spread dataset để ggplot nhận diện dữ liệu
rice_final <- rice_ok %>% tidyr::spread(key = "element", value = "value")

head(rice_final)
##   item year area_harvested production
## 1 Rice 1961       34694000   53494496
## 2 Rice 1962       35695008   49825552
## 3 Rice 1963       35809008   55497008
## 4 Rice 1964       36462000   58962000
## 5 Rice 1965       35470000   45883504
## 6 Rice 1966       35251008   45657008
na.omit(rice_final) -> rice_final_india

BƯỚC 6: FUNCTION HÓA ĐỒ THỊ

bubble_chart <- function(input_title, input_subtitle, rice_final, input_country){

p <- ggplot(rice_final, aes(
  x = year, y = production / 1000000,
  size = area_harvested / 1000,
  fill = production / area_harvested
)) +
  geom_point(alpha = 0.6, shape = 21, color = "black") +
  scale_size(range = c(1, 9), name = "Diện tích canh tác (nghìn ha)", 
             breaks = round(seq(min(rice_final$area_harvested)/1000, max(rice_final$area_harvested)/1000, length.out = 6),
                            digits = -2)) + 
  scale_fill_viridis_c(
    name = "Năng suất (tấn/ha)",
    option = "D",
    guide = "colourbar",
    breaks = round(seq(min(rice_final$production / rice_final$area_harvested), 
                 max(rice_final$production / rice_final$area_harvested), length.out = 6), digits = 1)
  ) +
  guides(
    size = guide_legend(order = 1),
    fill = guide_colourbar(barwidth = 1, barheight = 8, order = 2)
  ) +
  scale_x_continuous(expand = c(0, 0), limits = c(1955, 2025)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1.15 * max(rice_final$production) / 1000000)) +
  
  ## customize theme
  
  theme_classic() +
  ylab("Sản lượng thu hoạch (triệu tấn)") +
  xlab("Năm") +
  
  ## x and y lab
  theme(axis.title.x = element_text(face = "bold", size = 10)) +
  
  theme(axis.title.y = element_text(
    face = "bold",
    margin = margin(t = 1, r = 10, b = 1, l = 10)
  )) +
  
  ## tick text
  
  theme(
    axis.text.x = element_text(
      face = "bold", color = "#993333",
      size = 10
    ),
    axis.text.y = element_text(
      face = "bold", color = "#993333",
      size = 10
    )
  ) +
  
  ## title
  
  labs(title = input_title,
              subtitle = input_subtitle,
              caption = "Đồ thị: Duc Nguyen | tuhocr.com") +
  
  theme(
  plot.title = element_text(size = 12, 
                            face = "bold", 
                            hjust = 0, margin = margin(t = 10, r = 0, b = 5, l = 0)),
  plot.subtitle = element_text(size = 11, hjust = 0, margin = margin(t = 0, r = 0, b = 5, l = 0)),
  plot.caption = element_text(size = 9, hjust = 1, margin = margin(t = 0, r = 0, b = 5, l = 0))
) +
  
  ## legend

  theme(legend.position = "right", 
        legend.title = element_text(size = 10, face = "bold"),
        legend.background = element_blank(),
        legend.justification = c("right", "bottom"),
        # legend.box.just = "right",
        legend.margin = margin(t = 0, r = 0, b = 0, l = 5),
        legend.spacing.x = unit(0.5, "lines"),
        legend.spacing.y = unit(0.4, "lines")) +
  
  ## margin
  
  theme(
    panel.background = element_rect(fill = "white"),
    plot.margin = margin(t = 3, r = 14, b = 2, l = 2),
    plot.background = element_rect(fill = "grey90", colour = "black", linewidth = 1)
  ) +

  ## bổ sung thêm
# p + 

#   # annotate(geom = "text", x = 2018, y = 0.85*max(rice_final$production)/1000,
#   #          label = "Năm 2021 \nđạt 43,8 nghìn tấn",
#   #          color = "blue",
#   #          size = 3) +
  # geom_hline(yintercept = 48, colour = "#148c34", lwd = 1.2, lty = 2) +
  # annotate(
  #   geom = "text", x = 2000, y = 1.1 * max(rice_final$production) / 1000000,
  #   label = "Ngưỡng cực hạn 48 triệu tấn/năm",
  #   color = "red"
  # ) +
#   coord_cartesian(clip = "off")
  geom_label(aes(x = 1975, y = 0.9 * max(rice_final$production)/1000000, 
                 label = paste("Năm", max(rice_final$year), input_country, "\nthu hoạch được", round(rice_final[rice_final$year == max(rice_final$year), ]$production/1000000, digits = 0), "triệu tấn gạo", "\ntrên diện tích canh tác", round(rice_final[rice_final$year == max(rice_final$year), ]$area_harvested/1000, digits = 0), "nghìn ha", "\nnăng suất là", round(rice_final[rice_final$year == max(rice_final$year), ]$production/rice_final[rice_final$year == max(rice_final$year), ]$area_harvested, digits = 2), "tấn/ha")),
             fill = "#f5ec42", label.size = 0.15, size = 3.5)
## gọi đồ thị

return(p)
## add hình
# library(grid)
# library(png)
# image <- readPNG("rice.png")
# logor <- readPNG("logor.png")

}

Vẽ đồ thị cho Thái Lan

q <- bubble_chart(input_title = "Tình hình sản xuất lúa gạo ở Thái Lan từ 1960–2021",
             input_subtitle = "Nguồn: FAOSTAT [crop dataset, Thailand, flag A]",
             rice_final = rice_final_thailand,
             input_country = "Thái Lan")
q +
    geom_hline(yintercept = 40, colour = "#148c34", lwd = 1.2, lty = 2) +
  annotate(
    geom = "text", x = 2000, y = 1.1 * max(rice_final_thailand$production) / 1000000,
    label = "Ngưỡng cực hạn 40 triệu tấn/năm",
    color = "red")

image <- readPNG("thailand.png")
logor <- readPNG("logor.png")
grid.raster(image, x = 0.56, y = 0.24, width = 0.23)
grid.raster(logor, x = 0.9, y = 0.9, width = 0.1)

Sơ kết

Trên đây là hướng dẫn vẽ đồ thị từ nguồn FAOSTAT. Để học R bài bản từ A đến Z, thân mời Bạn tham gia khóa học “HDSD R để xử lý dữ liệu” để có nền tảng vững chắc về R nhằm tự tay làm các câu chuyện dữ liệu của riêng mình!

ĐĂNG KÝ NGAY: https://www.tuhocr.com/register

Hướng dẫn cài đặt package tuhocr https://tuhocr.github.io/