www.tuhocr.com
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
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_raw"
data_folder dir.create(data_folder) ## tạo folder lưu file trên máy tính
<- FAOsearch()
fao_metadata 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
<- get_faostat_bulk(code = "QCL", data_folder = data_folder)
crop_production
## 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")
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)
<- readRDS("data_raw/crop_production_e_all_data.rds")
crop_production 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
<- subset(crop_production, area == "Viet Nam")
vietnam_data 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ố
<- subset(vietnam_data, flag == "A")
vietnam_data_a dim(vietnam_data_a)
## [1] 4422 13
## chọn những cột cần thiết
<- vietnam_data_a[, c("item", "element", "year", "unit", "value")]
vietnam_crop
## sắp xếp theo năm và thông tin nông sản
<- vietnam_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
vietnam 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.
<- subset(vietnam, item == "Rice")
rice <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)
rice_ok
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_ok %>% tidyr::spread(key = "element", value = "value")
rice_final
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
library(ggplot2)
options(scipen = 99)
<- ggplot(rice_final, aes(
p 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)
<- readPNG("rice.png")
image <- readPNG("logor.png")
logor grid.raster(image, x = 0.56, y = 0.24, width = 0.23)
grid.raster(logor, x = 0.9, y = 0.9, width = 0.1)
<- subset(crop_production, area == "Thailand")
thailand_data 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ố
<- subset(thailand_data, flag == "A")
thailand_data_a dim(thailand_data_a)
## [1] 6441 13
## chọn những cột cần thiết
<- thailand_data_a[, c("item", "element", "year", "unit", "value")]
thailand_crop
## sắp xếp theo năm và thông tin nông sản
<- thailand_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
thailand 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.
<- subset(thailand, item == "Rice")
rice <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)
rice_ok
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_ok %>% tidyr::spread(key = "element", value = "value")
rice_final
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
<- subset(crop_production, area == "China, mainland")
china_data dim(china_data)
## [1] 30038 13
## subset dữ liệu Flag A từ nguồn chính thức do China công bố
<- subset(china_data, flag == "A")
china_data_a dim(china_data_a)
## [1] 5835 13
## chọn những cột cần thiết
<- china_data_a[, c("item", "element", "year", "unit", "value")]
china_crop
## sắp xếp theo năm và thông tin nông sản
<- china_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
china 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.
<- subset(china, item == "Rice")
rice <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)
rice_ok
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_ok %>% tidyr::spread(key = "element", value = "value")
rice_final
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
<- subset(crop_production, area == "India")
india_data dim(china_data)
## [1] 30038 13
## subset dữ liệu Flag A từ nguồn chính thức do China công bố
<- subset(india_data, flag == "A")
india_data_a dim(india_data_a)
## [1] 9846 13
## chọn những cột cần thiết
<- india_data_a[, c("item", "element", "year", "unit", "value")]
india_crop
## sắp xếp theo năm và thông tin nông sản
<- india_crop %>% dplyr::arrange(desc(year), desc(item), desc(element))
india 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.
<- subset(india, item == "Rice")
rice <- subset(rice, element != "yield")
rice_ok <- subset(rice_ok, select = -unit)
rice_ok
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_ok %>% tidyr::spread(key = "element", value = "value")
rice_final
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
<- function(input_title, input_subtitle, rice_final, input_country){
bubble_chart
<- ggplot(rice_final, aes(
p 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
<- bubble_chart(input_title = "Tình hình sản xuất lúa gạo ở Thái Lan từ 1960–2021",
q 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")
<- readPNG("thailand.png")
image <- readPNG("logor.png")
logor grid.raster(image, x = 0.56, y = 0.24, width = 0.23)
grid.raster(logor, x = 0.9, y = 0.9, width = 0.1)
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/