# 低分辨率绘图
# install.packages("scico")
# library(gcookbook)
# library(ggplot2)
# library(scico)
library(tidyverse)
library(sf)
library(raster)
library(terra)
library(dplyr)
library(ggspatial)
library(ggnewscale)
library(ggplot2)
library(showtext)
#配置字体
showtext_auto(enable = TRUE)
font_add("times",
regular = "times.ttf",
bold ="timesbd.ttf",
italic="timesbi.ttf",
bolditalic="timesi.ttf"
)
cnfont <- "times"
# 中国地图通常使用这样的坐标系
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("F:/R_codedata/使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）/使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）/provmapdata/minishp/chinaprov2021mini/chinaprov2021mini.shp") %>%
filter(!is.na(省代码)) -> provmap
# 线条
read_sf("F:/R_codedata/使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）/使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）/provmapdata/minishp/chinaprov2021mini/chinaprov2021mini_line.shp") %>%
filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
dplyr::select(class) -> provlinemap #
##循环
lapply(fs::dir_ls("F:/A 毕业论文/original data/season_means_raster/"), function(f){
print(f)
f
# 从变量f中提取文件名的前六个字符
file_name <- basename(f) # 获取文件名部分
month <- substr(file_name, 1, 6) # 提取前六个字符
# 打印提取的结果
cat("提取的前六个字符为:", month, "\n")
paste0("F:/A 毕业论文/original data/season_means_raster/",month,".pdf") -> destfile # 生成文件名
# 为了防止重复运行，加一个判断，如果文件已经存在则不需要再次运行
if(!file.exists(destfile)){
raster(f) %>%
rasterToPolygons() %>%
st_as_sf() -> cn2022polygons # 数据转为sf格式
cn2022polygons %>%
st_transform(mycrs) -> cn2022polygons
#小地图的坐标范围
small_bbox <- st_bbox(c(xmin = 120000,
xmax = 1766004.1,
ymax = 2557786.0,
ymin = 320000),
crs = st_crs(mycrs)) %>%
st_as_sfc()
# 提取小地图内的点
cn2022polygons %>%
st_intersection(small_bbox) -> cn2022polygons_small
# 移动小地图内的点到恰当位置：
cn2022polygons_small %>%
mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
sf::st_set_crs(mycrs) -> cn2022polygons_small
# 合并两个部分
bind_rows(cn2022polygons, cn2022polygons_small) -> cn2022polygons_all
##在数据上乘科学计数法
cn2022polygons_all[, 1] <- cn2022polygons_all[, 1] * 1000000000
cn2022polygons_all <- cn2022polygons_all %>%
set_names(c("dust", "geometry")) ##给变量的列重命名，有几列写几列
#查看数据类型
class(cn2022polygons_all)
#一步break转换操作
# 使用 cut 函数切割 dust 列
cn2022polygons_all$dust <- cut(cn2022polygons_all$dust,
breaks = c(0, 25, 50, 100, 200, 300, 400, 500, 600, 700, Inf),
labels = c("0~25", "26~50", "51~100", "101~200", "201~300",
"301~400", "401~500", "501~600", "601~700", ">700"),
include.lowest = TRUE)
##
cn2022polygons_all <- st_as_sf(cn2022polygons_all)
# 提取中国范围的
cn2022polygons_all %>%
st_intersection(provmap) -> cn2022polygons_all
print(cn2022polygons_all)
ggplot(provmap) +
geom_sf(data = cn2022polygons_all,
aes(fill = dust), #比例尺
linewidth = 0.001,color="white") +
scale_fill_manual(values = c("0~25" = "white",
"26~50" = "#FFF7EC",
"51~100" = "#FEE8C8",
"101~200" = "#FDD49E",
"201~300" = "#FDBB84",
"301~400" = "#FC8D59",
"401~500" = "#EF6548",
"501~600" = "#D7301F",
"601~700" = "#B30000",
">700" = "#7F0000"),
name = "DUSTMASS\nµg/m3") +
new_scale_color() +
geom_sf(fill = NA, color = "gray", linewidth = 0.01) +
geom_sf(data = provlinemap,
aes(color = class, linewidth = class),
show.legend = F) +
scale_color_manual(
values = c("九段线" = "#A29AC4",
"海岸线" = "#0055AA",
"小地图框格" = "black")
) +
scale_linewidth_manual(
values = c("九段线" = 0.6,
"海岸线" = 0.3,
"小地图框格" = 0.3)
) +
annotation_scale(location = "bl",
pad_x=unit(1,"cm"), pad_y=unit(1,"cm"),
width_hint = 0.3,
text_family = cnfont,
text_cex = 2)+#字体大小
theme_classic(base_size = 30,#边框字体大小
base_family = "",
base_line_size = 0.5,
base_rect_size = 0.5)+
theme(
panel.border = element_rect(colour = "black",
fill = NA,
linewidth = 1), #加边框
axis.title.x = element_blank(),#X轴标题
axis.title.y = element_blank(),#y轴标题
#图例的相对位置，注释掉就在外面
# legend.position = c(0.89, 0.53),#图例的相对位置（x比例,y比例
legend.fill = NULL,
legend.title.col=NULL,
legend.title = (element_text(size=16)),
legend.text = (element_text(size=16)),
legend.background = (element_rect(fill = NULL)) ,
plot.background = element_rect(fill = "white", color = "white")) +
annotation_north_arrow(
location = "tl",
which_north = "false",
pad_y = unit(0.5, "cm"),
pad_x = unit(1, "cm"),
style = north_arrow_nautical(
text_family = cnfont,
text_size=20
)
)-> p
ggsave(filename = paste0("F:/A 毕业论文/original data/season_means_raster/", month, ".pdf"), plot = p, device = "pdf", width = 10, height = 9)
}
})
