UPdating:
2020-09-05，更新产量数据为2018年最新统计数据！
2019-10-05，庆祝建国70周年，地图修改为中国红！
最近，Boss交给了一个任务，把各国凡纳滨对虾的产量在世界地图上标注一下。
对于生成地图还是很有兴趣的。没想到竟然很有挑战性！
之前看过几个blog，印象中有不少包可以完成这个任务。
但最后发现，还是ggplot2包比较好用。为了画地图，用到的包主要有：
- maps 提供世界地图数据。
- ggrepel 重排地图上的国家名称，消除名称间的重叠现象。
- showtext和sysfonts 确保在生成的地图上可以使用中文字体如黑体。
- data.table 数据清洗，生成符合要求的地图数据集，供ggplot2调用。
1. 获取世界各国凡纳滨对虾产量数据
从FAO下载世界各国所有水产养殖种类的产量数据。选择下载csv格式数据集。
解压缩后，产量数据在TS_FI_AQUACULTURE.csv文件中。凡纳滨对虾的产量数据，可以通过其3位种类代码是“PNV”进行查找。当前可以获得的最新即2018年的产量数据。
world_aquaculture_production_dt <-
fread(
input = "datasets/Aquaculture_2019.1.0/TS_FI_AQUACULTURE.csv",
header = TRUE,
stringsAsFactors = FALSE,
encoding = "UTF-8"
)
species_3alpha_code <- "PNV"
world_PNV_production_dt <-
world_aquaculture_production_dt[SPECIES %chin% species_3alpha_code]
world_PNV_production_2018_dt <-
world_PNV_production_dt[YEAR == 2018]凡纳滨对虾在海水、淡水和半咸水中都可以养殖，因此在统计世界各个国家产量时，需要累加在不同养殖环境下的产量.
由于在TS_FI_AQUACULTURE.csv文件中，国家列（COUNTRY）是用代码表示。因此需要从CL_FI_COUNTRY_GROUPS.csv文件中匹配UN_Code列查询国家名称等信息。其中Name_En和Name_Cn分别表示英文和中文国家名称。
#读取国家信息
country_group_dt <-
fread(
input = "datasets/Aquaculture_2019.1.0/CL_FI_COUNTRY_GROUPS.csv",
header = TRUE,
stringsAsFactors = FALSE,
encoding = "UTF-8"
)
world_PNV_production_2018_country_info_dt <-
merge(
world_PNV_production_2018_dt,
country_group_dt[, c("UN_Code",
"Identifier",
"ISO2_Code",
"ISO3_Code",
"Name_En",
"Name_Cn")],
by.x = c("COUNTRY"),
by.y = c("UN_Code")
)
world_PNV_production_2018_country_info_dt[1:5]## COUNTRY PRODUCTION_AREA ENVIRONMENT SPECIES YEAR QUANTITY_UNIT QUANTITY
## 1: 76 41 3 PNV 2018 t 62000.000
## 2: 84 31 2 PNV 2018 t 500.000
## 3: 156 4 1 PNV 2018 t 642807.000
## 4: 156 61 2 PNV 2018 t 1117534.000
## 5: 158 4 1 PNV 2018 t 3721.258
## QUANTITY_SYMBOL VALUE VALUE_SYMBOL Identifier ISO2_Code ISO3_Code
## 1: F 407244.42 NA 21 BR BRA
## 2: F 2000.00 NA 23 BZ BLZ
## 3: 5818688.96 NA 41 CN CHN
## 4: 7810445.13 NA 41 CN CHN
## 5: 24906.66 NA 214 TW TWN
## Name_En Name_Cn
## 1: Brazil 巴 西
## 2: Belize 伯利兹
## 3: China 中华人民共和国
## 4: China 中华人民共和国
## 5: Taiwan Province of China# 累加每个国家的总产量，凡纳滨对虾可在不同环境下养殖，如海、淡水、半咸水等
PNV_production_2018_per_country_dt <-
world_PNV_production_2018_country_info_dt[, .(Production = sum(QUANTITY)), keyby =
c("Name_En")]
#按照产量降序排列
setorder(PNV_production_2018_per_country_dt,-Production)
setnames(PNV_production_2018_per_country_dt,
c("Name_En"),
c("Country"))
PNV_production_2018_per_country_dt## Country Production
## 1: China 1760341.000
## 2: Indonesia 708680.000
## 3: India 622000.000
## 4: Ecuador 510000.000
## 5: Viet Nam 475000.000
## 6: Thailand 347258.000
## 7: Mexico 157934.361
## 8: Brazil 62000.000
## 9: Saudi Arabia 56100.000
## 10: Iran (Islamic Rep. of) 47859.000
## 11: Malaysia 36007.250
## 12: Honduras 31500.000
## 13: Peru 29716.930
## 14: Nicaragua 29458.400
## 15: Venezuela, Boliv Rep of 24500.000
## 16: Guatemala 17273.000
## 17: Philippines 12527.000
## 18: Taiwan Province of China 7978.804
## 19: Cuba 6474.000
## 20: Panama 6409.000
## 21: Colombia 5397.000
## 22: Korea, Republic of 4480.000
## 23: Costa Rica 3000.000
## 24: United States of America 2035.000
## 25: El Salvador 1150.000
## 26: Belize 500.000
## 27: Dominican Republic 400.000
## 28: Egypt 155.000
## 29: Suriname 30.000
## 30: Northern Mariana Is. 23.000
## 31: Lebanon 15.000
## 32: Spain 15.000
## 33: Singapore 12.600
## 34: Guam 10.000
## 35: Tunisia 1.200
## Country Production2.绘制地图
2.1 获取地图数据
利用map_data()函数从maps包中获取地图数据。其中long和lat分别表示经、纬度。
地图数据中region列与国家对应。group这一列很重要，在绘图时需要指定。
#获取世界地图数据
map_world_dt <- as.data.table(map_data('world'))
map_world_dt[1:20]## long lat group order region subregion
## 1: -69.89912 12.45200 1 1 Aruba <NA>
## 2: -69.89571 12.42300 1 2 Aruba <NA>
## 3: -69.94219 12.43853 1 3 Aruba <NA>
## 4: -70.00415 12.50049 1 4 Aruba <NA>
## 5: -70.06612 12.54697 1 5 Aruba <NA>
## 6: -70.05088 12.59707 1 6 Aruba <NA>
## 7: -70.03511 12.61411 1 7 Aruba <NA>
## 8: -69.97314 12.56763 1 8 Aruba <NA>
## 9: -69.91181 12.48047 1 9 Aruba <NA>
## 10: -69.89912 12.45200 1 10 Aruba <NA>
## 11: 74.89131 37.23164 2 12 Afghanistan <NA>
## 12: 74.84023 37.22505 2 13 Afghanistan <NA>
## 13: 74.76738 37.24917 2 14 Afghanistan <NA>
## 14: 74.73896 37.28564 2 15 Afghanistan <NA>
## 15: 74.72666 37.29072 2 16 Afghanistan <NA>
## 16: 74.66895 37.26670 2 17 Afghanistan <NA>
## 17: 74.55899 37.23662 2 18 Afghanistan <NA>
## 18: 74.37217 37.15771 2 19 Afghanistan <NA>
## 19: 74.37617 37.13735 2 20 Afghanistan <NA>
## 20: 74.49796 37.05722 2 21 Afghanistan <NA>subregion貌似与州或者省份对应。我们来看一下美国的各个州。貌似subregion列并没有包括全部50个州。
unique(map_world_dt[region %chin% c("USA")]$subregion)## [1] "Hawaii" "Florida" "13"
## [4] "16" "Texas" "Louisiana"
## [7] "25" "Mississippi" "Alabama"
## [10] "30" "California" "37"
## [13] "Core Banks" "Ocracoke Island" "Pea Island"
## [16] "41" "Chincoteague Island" "Long Beach Island"
## [19] "44" "New York" "Massachusettes"
## [22] "48" "49" "52"
## [25] "Maine" "67" "69"
## [28] "70" "Washington" "73"
## [31] NA "Alaska"对地图数据进行清洗，主要包括：
- 将台湾明确为中国的一个省。
- 在地图中删除南极洲数据。因为凡纳滨对虾为高温物种，不可能在此养殖。
- 美国阿拉斯基地区靠近北极圈，温度太低，也不可能养殖凡纳滨对虾。因此在绘制地区时，不把该地区标注为可养殖地区。
#把台湾明确为台湾省
map_world_dt[region %chin% c("Taiwan"), region := c("Taiwan Province of China")]
#因为凡纳滨对虾不可能在南极洲生存，直接在地图数据中剔除南极洲，不绘制.
map_world_dt <- map_world_dt[!region %chin% c("Antarctica")]
#美国阿拉斯加靠近北极，不可能养凡纳滨对虾，这部分地区单独设置为一个国家，
#这样绘制美国地图时，不包括该地区
map_world_dt[subregion %chin% c("Alaska"), region:= c("USANoDraw")]2.2 统一FAO数据与maps包中的国家名称
从FAO下载的数据，部分国家名跟maps包中的国家名不一致。需要找出这些国家，并把FAO数据中国家名称修改为maps包国家名称。
把凡纳滨对虾产量数据合并到地图数据集中。
计算每个国家的经纬度均值，生成一个新的数据集，用于在地图中定位国家名字的显示位置。
#查找名称不一致的国家
FAO_nomatched_country_name_v <-
PNV_production_2018_per_country_dt$Country[!(PNV_production_2018_per_country_dt$Country %in% unique(map_world_dt$region))]
FAO_nomatched_country_name_v## [1] "Viet Nam" "Iran (Islamic Rep. of)"
## [3] "Venezuela, Boliv Rep of" "Korea, Republic of"
## [5] "United States of America" "Northern Mariana Is."#根据maps包中的国家名，重新对相应的国家重新命名
maps_matched_country_name_v <-
c(
"Vietnam",
"Iran",
"Venezuela",
"Korea",
"USA",
"Northern Mariana Islands"
)
PNV_production_2018_per_country_dt[Country %chin% FAO_nomatched_country_name_v,
Country := maps_matched_country_name_v]
#把凡纳滨对虾产量数据跟地图数据合并
map_world_PNV_production_2018_dt <-
merge(
map_world_dt,
PNV_production_2018_per_country_dt,
all.x = TRUE,
by.x = "region",
by.y = "Country"
)
#将产量为0的国家设置为NA，这样绘制地图时，这些国家不会上色
map_world_PNV_production_2018_dt[Production == 0, Production:=NA]
#计算国家经纬度均值，用于定位国家名称的显示位置
map_world_PNV_country_position_2018_dt <-
map_world_PNV_production_2018_dt[Production > 0,
.(MeanLong = mean(long),
MeanLat = mean(lat),
group = mean(group)),
by = c("region")]
map_world_PNV_country_position_2018_dt[1:5]## region MeanLong MeanLat group
## 1: Belize -88.46050 17.421947 218.5692
## 2: Brazil -55.28057 -9.802375 237.0159
## 3: China 106.84757 35.082435 434.6369
## 4: Colombia -72.76287 3.471105 437.9872
## 5: Costa Rica -84.03994 9.756017 450.00002.2 绘制凡纳滨对虾世界养殖分布图
地图中的文字使用黑体。因此需要利用showtext包加载黑体文件。
调用ggplot函数和geom_polygon函数来绘制地图。
为了消除可能存在的国家名重叠的问题，调用ggrepel包中geom_text_repel()函数优化国家名在地图中的布局。
在地图中，颜色的深浅，表示产量的高低。产量越高，颜色越深。
国家名，根据map_world_PNV_country_position_2018_dt数据集中的平均经纬度定位显示，并且加了一个蓝色圆点。
为了正确的画出地图，需要在aes函数中设置group=group，这个特别重要。
font_add(family = "heiti", regular = "simhei.ttf")
showtext_auto()
ggplot(map_world_PNV_production_2018_dt,
aes(x = long, y = lat, group = group)) +
geom_polygon(aes(fill = Production),color="grey50") +
#国家绘制为朱红色
scale_fill_gradient(
low = "#ffdfd5", ##f3f8f1 绿色
high = "#d83010", ##228b22 橙色
breaks = c(5000, 10000, 50000, 250000, 500000, 1600000),
space = "Lab",
na.value = "#a9a9a9",
guide = "colourbar",
aesthetics = "fill"
) +
guides(fill = guide_legend(reverse = T)) +
labs(
fill = "产量(吨)/年",
title = "凡纳滨对虾世界养殖分布",
subtitle = "主要养殖国家年产量",
x = NULL,
y = NULL
) +
geom_point(aes(x =MeanLong, y= MeanLat),
color = "green",
size = 1.5,
data=map_world_PNV_country_position_2018_dt) +
geom_text_repel(
aes(x = MeanLong, y = MeanLat, label = region),
data = map_world_PNV_country_position_2018_dt,
size = 8
) +
theme(
text = element_text(family = "heiti", color = "#000000"),
plot.title = element_text(family = "heiti", hjust=0.5, size = 36),
plot.subtitle = element_text(family = "heiti", hjust=0.5, size = 28),
axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid = element_blank(),
panel.background = element_rect(fill = "#FFFFFF"),
plot.background = element_rect(fill = "#FFFFFF"),
legend.position = c(0.18, 0.20),
legend.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(family = "heiti", size = 18),
legend.text=element_text(family="heiti", size=18)
) +
annotate(
geom = 'text',
label = '数据来源：FAO, 2020. http://www.fao.org/fishery/statistics/global-aquaculture-production/zh',
x = 10, y = -55,
size = 4,
family = 'heiti',
color = 'grey50',
hjust = 'left'
)