ggplot2绘制时间序列变化图

疫情时间序列变化

涉及主要方法：

• list.files遍历文件
• read_excel读取EXCEL文件
• for循环读取文件夹内所有EXCEL
• ymd时间提取
• rbind/cbind数据框行列的合并
• ggplot2制图
• labs标注X轴、Y轴标题
• theme_bw使用空白模版可视化效果
• theme标题、坐标轴、图例等文字格式、文字大小、位置设置
• geom_smooth平滑曲线绘制
• geom_point散点绘制

library(lubridate)
library(ggplot2)
library(grid)
library(dplyr)
library(sf)
library(knitr)
library(tidyverse)


a = list.files("Blocks")                                             #list.files命令将input文件夹下所有文件名输入a
dir = paste("./Blocks/",a,sep="")                                      #用paste命令构建路径变量dir
n = length(dir)                                                                 #读取dir长度，也就是文件夹下的文件个数
mergedata2 = read_excel(path =  dir[1],sheet = 1, col_names = T)   #读入第一个文件内容（可以不用先读一个，但是为了简单，省去定义data.frame的时间，我选择先读入一个文件。
for (i in 2:n){
  new.data = read_excel(path =  dir[i],sheet = 1, col_names = T)
  mergedata2 = rbind(mergedata2,new.data)
}
Daydate <- ymd(mergedata2$数据时间)
BJCOVID <- cbind(mergedata2, Daydate)
BJTOWNCOVID <- BJCOVID[BJCOVID$town != "北京", ]
BJTotal <- BJCOVID[BJCOVID$town == "北京", ]
BJNewConform <- BJTotal[,c(1:7,13)]
colnames(BJNewConform)[7] <- "Conform"
ConTYPE<- "新增确诊"
BJConform <- cbind(BJNewConform, ConTYPE)
BJTotalCon <- BJTotal[,c(1:6,8,13)]
colnames(BJTotalCon)[7] <- "Conform"
ConTYPE<- "累计确诊"
BJConform <- rbind(BJConform,cbind(BJTotalCon, ConTYPE))

BJConform <- read.csv("北京确诊.csv", header = T)
BJConform$Daydate <- ymd(BJConform$Daydate)
#北京新增累计曲线图
p1 <- ggplot(BJConform, aes(Daydate, Conform, color=ConTYPE))+
  labs(x="时间", y="病例数", title="北京市新发疫情病例数变化图")+
  theme_bw()+
  theme(title = element_text(size = 22),
        axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18), 
        legend.position = "bottom",
        legend.title = element_blank())+     #y轴标签大小
  geom_smooth()+
  geom_point(size=2)
p1

p1

#报告疫情的街道数量统计与制图
nTownNewCon <-  BJTOWNCOVID%>% 
  filter(!is.na(新增确诊))%>%
  group_by(Daydate)%>%
  summarise(n())
ConTYPE<- "新增确诊"
nTownNewCon <- cbind(nTownNewCon, ConTYPE)
nTownTotCon <-  BJTOWNCOVID%>% 
  filter(!is.na(累计确诊))%>%
  group_by(Daydate)%>%
  summarise(n())
ConTYPE<- "累计确诊"
nTownTotCon <- cbind(nTownTotCon, ConTYPE)
nTown <- rbind(nTownNewCon, nTownTotCon)
#街道数
p2 <- ggplot(nTown, aes(Daydate, `n()`, color=ConTYPE))+
  labs(x="时间", y="报告疫情街道（乡镇）数", title="报告新增/累计确诊病例的街道数量变化图")+
  theme(title = element_text(size = 22),
        axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18), 
        legend.position = "bottom",
        legend.title = element_blank())+     #y轴标签大小
  geom_smooth()+
  geom_point(size=2)
p2

p2

需要注意的是，geom_smooth默认使用的是"loess"回归方法，是一种非参回归方法，无法给出回归方程。

loess（locally weighted regression）是一种用于局部回归分析的非参数方法，它主要是把样本划分成一个个小区间，对区间中的样本进行多项式拟合，不断重复这个过程得到在不同区间的加权回归曲线，最后再把这些回归曲线的中心连在一起合成完整的回归曲线.

NDVI时间序列拟合回归

对NDVI时间序列进行线性回归，并给出回归方程、R2和p值。

• ggplot2绘图
• ggpmisc进行回归方程标注
• stat_poly_eq做标注
• formula使用了一个变量，指定回归方程形式

library(ggplot2)
library(ggpmisc)
#年度NDVI曲线
NDVI.formula <- y~x
p2 <- ggplot(data = Year_Max_NDVI, aes(Date, avgNDVI))+
  labs(x="Year", y="NDVI")+
  theme_bw()+
  geom_smooth(method = "lm", formula = NDVI.formula) +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), stat(p.value.label), sep = "*\", \"*")),
               formula = NDVI.formula, parse = TRUE)+
  geom_point()
p2

NDVI线性回归

1. https://ggplot2.tidyverse.org/reference/geom_smooth.html
2. https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/loess
3. https://blog.csdn.net/jesseyule/article/details/95245350