欢迎各位培训班小伙伴参加明晚 8 点的直播课：「使用 R 语言绘制填充地图、栅格地图 + 等降水量线」～
最近有个小伙伴问到了关于在栅格地图上绘制等高线的问题。今天我们就以降水量栅格数据为例讲解如何在 R 语言中绘制这种地图。
在学习本课程之前，需要预先学习之前的课程：
使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）：https://rstata.duanshu.com/#/brief/course/379d19770956478ebc4d919910fca74c
附件中提供了 2021 年各月的降水量栅格数据：pre_2021.nc，数据来源于 中国1km分辨率逐月降水量数据集（1901-2023）：https://data.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2 。
填充地图
首先我们来绘制各省份的年均降水量填充地图。
读取栅格数据和省级矢量数据：
library(tidyverse)
library(sf)
library(terra)
# 读取降水量栅格数据
rast("pre_2021.nc") -> rst
rst
#> class : SpatRaster
#> dimensions : 5146, 7849, 12 (nrow, ncol, nlyr)
#> resolution : 0.008333333, 0.008333333 (x, y)
#> extent : 71.28589, 136.6942, 15.75212, 58.63546 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source : pre_2021.nc
#> varname : pre (monthly precipitation)
#> names : pre_1, pre_2, pre_3, pre_4, pre_5, pre_6, ...
# 年均
app(rst, mean, na.rm = T) -> rst
rst
#> class : SpatRaster
#> dimensions : 5146, 7849, 1 (nrow, ncol, nlyr)
#> resolution : 0.008333333, 0.008333333 (x, y)
#> extent : 71.28589, 136.6942, 15.75212, 58.63546 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> name : mean
#> min value : 4.166667
#> max value : 6252.083333
# 提取中国范围的
read_sf("2021行政区划/省.shp") %>%
st_union() %>%
nngeo::st_remove_holes() %>%
st_sfc() -> cn
plot(cn)
rst %>%
terra::crop(vect(cn)) %>%
terra::mask(vect(cn)) -> rst
# 读取省级行政区划
read_sf("2021行政区划/省.shp") %>%
select(-contains("类型")) -> prov
prov
#> Simple feature collection with 34 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 73.50235 ymin: 3.83703 xmax: 135.0957 ymax: 53.56362
#> Geodetic CRS: WGS 84
#> # A tibble: 34 × 3
#> 省 省代码 geometry
#> <chr> <dbl> <MULTIPOLYGON [°]>
#> 1 安徽省 340000 (((116.426 34.64413, 116.4263 34.64, 116.4261 34.63816…
#> 2 澳门特别行政区 820000 (((113.5815 22.19861, 113.5812 22.19779, 113.5789 22.1…
#> 3 北京市 110000 (((116.6255 41.05897, 116.6258 41.05887, 116.6264 41.0…
#> 4 福建省 350000 (((117.689 23.53157, 117.6888 23.53144, 117.6875 23.53…
#> 5 甘肃省 620000 (((106.0718 35.44995, 106.0713 35.44955, 106.0707 35.4…
#> 6 广东省 440000 (((110.5657 20.33342, 110.5604 20.32783, 110.5565 20.3…
#> 7 广西壮族自治区 450000 (((109.2167 20.90438, 109.2144 20.90184, 109.2108 20.9…
#> 8 贵州省 520000 (((105.0948 24.92524, 105.0947 24.9252, 105.0946 24.92…
#> 9 海南省 460000 (((112.0522 3.838239, 112.0483 3.837329, 112.0181 3.83…
#> 10 河北省 130000 (((118.2715 38.98381, 118.2715 38.982, 118.2707 38.979…
#> # ℹ 24 more rows
分区域计算平均降水量：
terra::extract(rst, vect(prov), fun = mean, na.rm = T) %>%
as_tibble() %>%
bind_cols(st_drop_geometry(prov)) %>%
select(-ID) %>%
rename(年度降水量_mm = mean) -> df
df
#> # A tibble: 34 × 3
#> 年度降水量_mm 省 省代码
#> <dbl> <chr> <dbl>
#> 1 1119. 安徽省 340000
#> 2 1597. 澳门特别行政区 820000
#> 3 519. 北京市 110000
#> 4 1299. 福建省 350000
#> 5 286. 甘肃省 620000
#> 6 1271. 广东省 440000
#> 7 1227. 广西壮族自治区 450000
#> 8 1071. 贵州省 520000
#> 9 1279. 海南省 460000
#> 10 485. 河北省 130000
#> # ℹ 24 more rows
然后就可以绘制填充地图了，这部分代码类似之前的课程：
# 绘制填充地图
# 中国地图通常使用这样的坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
# 读取小地图版本的中国省级地图
read_sf("chinaprov2019mini/chinaprov2019mini.shp") %>%
filter(!is.na(省代码)) -> provmap
# 线条
read_sf("chinaprov2019mini/chinaprov2019mini_line.shp") %>%
filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
select(class) -> provlinemap
# 小地图的范围
small_bbox <- st_bbox(c(xmin = 120000,
xmax = 1766004.1,
ymax = 2557786.0,
ymin = 320000),
crs = st_crs(mycrs)) %>%
st_as_sfc()
# 绘图
library(ggspatial)
library(ggnewscale)
# 绘制填充地图
# 和地图数据合并
provmap %>%
left_join(df) -> provmap2
hist(provmap2$年度降水量_mm)
# 绘图
provmap2 %>%
mutate(年度降水量_mm = if_else(is.na(年度降水量_mm), -400,
年度降水量_mm),
group = cut(年度降水量_mm,
breaks = seq(-400, 2400, by = 400), include.lowest = T,
labels = c("No data", "0 ~ 400", "400 ~ 800", "800 ~ 1200",
"1200 ~ 1600", "1600 ~ 2000", "2000 ~ 2400"))) %>%
# st_drop_geometry() %>%
# count(group)
ggplot() +
geom_sf(aes(fill = group),
color = "gray", linewidth = 0.1) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA",
"小地图框格" = "black",
"省份" = "black")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3,
"小地图框格" = 0.3,
"省份" = 0.3)
) +
scale_fill_manual(values = c("white", c("#d53e4f", "#fc8d59",
"#fee08b", "#e6f598",
"#99d594", "#3288bd"))) +
new_scale_color() +
stat_sf_coordinates(data = provmap,
geom = "text", color = "gray",
aes(label = 省), family = cnfont,
fun.geometry = st_point_on_surface,
size = 3) +
annotation_scale(location = "bl",
width_hint = 0.3,
text_family = cnfont) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.12, 0.2),
plot.background = element_rect(fill = "white",
color = "white")) +
guides(fill = guide_legend(ncol = 2)) +
labs(fill = "Precipitation\n(thounds mm)") +
annotation_north_arrow(
location = "tr",
which_north = "false",
pad_y = unit(0.1, "cm"),
style = north_arrow_fancy_orienteering(
text_family = cnfont
)
) -> p1
ggsave("pic1.png", width = 10, height = 8.5, device = png)
栅格地图
栅格地图的绘制方法也和之前课程中讲解的一样。栅格数据绘制地图需要转换成点或者多边形，所以这里我们改用 raster 包处理栅格数据：
# 栅格数据
library(raster)
raster(rst) -> rst
# 计算等降水量线
rasterToContour(rst, levels = seq(0, 2400, by = 400)) %>%
st_as_sf() -> line
line
#> Simple feature collection with 6 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension: XY
#> Bounding box: xmin: 78.54866 ymin: 18.43559 xmax: 132.7324 ymax: 49.96448
#> Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
#> level geometry
#> C_1 400 MULTILINESTRING ((78.70524 ...
#> C_2 800 MULTILINESTRING ((79.12219 ...
#> C_3 1200 MULTILINESTRING ((89.0151 2...
#> C_4 1600 MULTILINESTRING ((92.68309 ...
#> C_5 2000 MULTILINESTRING ((117.7477 ...
#> C_6 2400 MULTILINESTRING ((120.5423 ...
plot(line)
# 栅格数据转点（高分辨率）或者多边形（低分辨率），这里使用多边形试试：
rst %>%
aggregate(fact = 10) %>% # 降低分辨率，提升绘图速度
rasterToPolygons() %>%
st_as_sf() -> rstpoints
rstpoints %>%
st_transform(mycrs) -> rstpoints
# 提取小地图内的点
rstpoints %>%
st_intersection(small_bbox) -> rstpoints_small
# 主图中的点
st_bbox(c(xmin = -2725586,
xmax = 2982768,
ymax = 6000000,
ymin = 1800655),
crs = st_crs(mycrs)) %>%
st_as_sfc() -> plotbbox
rstpoints %>%
st_intersection(plotbbox) -> rstpoints_main
# 移动小地图内的点到恰当位置：
rstpoints_small %>%
mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
sf::st_set_crs(mycrs) -> rstpoints_small
# 合并两个部分
bind_rows(rstpoints_main, rstpoints_small) -> rstpoints_all
line %>%
st_transform(mycrs) -> line
# 绘图
ggplot(provmap) +
geom_sf(data = rstpoints_all, aes(fill = mean),
linewidth = 0.001,
color = NA) +
scale_fill_gradientn(colours = c("#9e0142", "#d53e4f",
"#f46d43", "#fdae61",
"#fee08b", "#ffffbf",
"#e6f598", "#abdda4",
"#66c2a5", "#3288bd",
"#5e4fa2"),
name = "Precipitation\n(thounds mm)") +
new_scale_color() +
geom_sf(fill = NA, color = "gray",
linewidth = 0.03) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
geom_sf(data = line, color = "white", linewidth = 0.5) +
stat_sf_coordinates(data = provmap,
geom = "text", color = "gray",
aes(label = 省), family = cnfont,
fun.geometry = st_point_on_surface,
size = 3) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA",
"小地图框格" = "black")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3,
"小地图框格" = 0.3)
) +
annotation_scale(location = "bl",
width_hint = 0.3,
text_family = cnfont) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.09, 0.2),
legend.direction = "horizontal",
legend.title.position = "top") +
annotation_north_arrow(
location = "tr",
which_north = "false",
pad_y = unit(0.1, "cm"),
style = north_arrow_fancy_orienteering(
text_family = cnfont
)
) -> p2
ggsave("pic2.png", width = 10, height = 8.5, device = png)
这样我们就把等降水线添加到了地图上。
添加等降水线文本标签
如果想把等降水量线的文本标签标注到线上，可以借助 geomtextpath 包。
安装：
remotes::install_github("AllanCameron/geomtextpath")
然后就可以绘制了：
library(geomtextpath)
line %>%
st_cast("LINESTRING") %>%
mutate(level = paste0(level, " mm")) -> line2
ggplot(provmap) +
geom_sf(data = rstpoints_all, aes(fill = mean),
linewidth = 0.001,
color = NA) +
scale_fill_gradientn(colours = c("#9e0142", "#d53e4f",
"#f46d43", "#fdae61",
"#fee08b", "#ffffbf",
"#e6f598", "#abdda4",
"#66c2a5", "#3288bd",
"#5e4fa2"),
name = "Precipitation\n(thounds mm)") +
new_scale_color() +
geom_sf(fill = NA, color = "gray",
linewidth = 0.03) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
geom_textsf(data = line2, aes(label = level),
family = cnfont,
linecolour = "white", color = "gray30",
text_smoothing = 50, size = 3) +
stat_sf_coordinates(data = provmap,
geom = "text", color = "gray",
aes(label = 省), family = cnfont,
fun.geometry = st_point_on_surface,
size = 3) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA",
"小地图框格" = "black")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3,
"小地图框格" = 0.3)
) +
annotation_scale(location = "bl",
width_hint = 0.3,
text_family = cnfont) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = c(0.09, 0.2),
legend.direction = "horizontal",
legend.title.position = "top") +
annotation_north_arrow(
location = "tr",
which_north = "false",
pad_y = unit(0.1, "cm"),
style = north_arrow_fancy_orienteering(
text_family = cnfont
)
) -> p3
ggsave("pic3.png", width = 10, height = 8.5, device = png)
直播信息
为了让大家更好的理解本文内容，欢迎各位培训班会员参加明晚 8 点的直播课：「使用 R 语言绘制填充地图、栅格地图 + 等降水量线」。
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
