R语言 热浪（Heatwave）识别（ERA5数据）
R语言初学者，欢迎大家交流学习
本代码基于R语言
重点的包为heatwaveR包
使用的气温数据为 ERA5-Land Daily Aggregated
下载自GEE
https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_LAND_DAILY_RAW
下面是代码
#清空变量
rm(list = ls())
gc()
#安装包
#install.packages("terra")
#install.packages("heatwaveR")
#install.packages("magrittr")
#install.packages("tidyverse")
#install.packages("tidyr")
#install.packages("lubridate")
#install.packages("dplyr")
#加载包
library(magrittr)
library(terra)
library(heatwaveR)
library(tidyverse,warn.conflicts = FALSE)
library(ggplot2)
library(tidyr)
library(scales)
library(lubridate)
library(trend)
library(dplyr)
#加载矢量边界数据
iliPath <- "E:/Desktop/ili_riverbasin_1984/ili_riverbasin_1984.shp"
ili_vect <- vect(iliPath)
#文件夹路径
path = "E:/Desktop/ERA5_Daily_ili _1"
#提取所有文件路径
dataPath <- list.files(path, recursive = TRUE, pattern = ".tif$",
full.names = TRUE, include.dirs = TRUE)
#提取文件时间
dataname <- substr(dataPath,46,53) %>%
as.character() %>%
as.Date(.,"%Y%m%d")
#提取复合图层中的温度层
era_tem_2m <- rast()
##温度数据
for (i in 1:length(dataPath)) {
#lyrs=2是2米大气温度
heatdata <- rast(dataPath[i],lyrs=2) %>%
crop(.,ili_vect) %>%
mask(.,ili_vect)
#为每层数据命名
names(heatdata) <- substr(dataPath[i],46,53)
#将每层数据添加到数据集中
add(era_tem_2m) <- heatdata
}
era_tem_2m
#保存结果
filename <- "E:/Desktop/ERA5_air_tem_2m.tif"
writeRaster(era_tem_2m, filename = filename , overwrite=TRUE)
##################################################################
#全局计算
##计算平均值
tempData <- global(era_tem_2m,"mean",na.rm=T,cores=4) %>%
as.data.frame(.)
#准备热浪识别数据
tempData$date <- dataname
tempData$mean <- tempData$mean - 273.15
names(tempData) <- c("temp","t")
#写文件
filename <- "E:/Desktop/ili_global_mean.csv"
write.csv(tempData,file = filename)
#绘图
ggplot(data = tempData, mapping = aes(x=t,y=temp, group = 1))+
xlab("year")+geom_line(size=0.5)+
scale_x_date(breaks = "8 years",date_labels="%Y")
#######################
#detect_heatwave
# Make a climatology from a daily time series
heatwave <- ts2clm(tempData, x=t,y=temp,
climatologyPeriod = c("1963-03-08","1993-03-08"),
pctile = 90 ,roundClm = 4) %>%
# Detect heatwaves
detect_event(., x=t, y=temp,
minDuration = 3,
maxGap = 2,
categories = TRUE)
heatwave
heatwave$event_no[length(heatwave$event_no)]
sum(heatwave$duration)
mean(heatwave$duration)
max(heatwave$intensity_cumulative)
# Create a line plot of heatwaves .
event_line(heatwave,
min_duration = 5,
metric = "intensity_cumulative",
category = TRUE,
spread = 150,
start_date = "2000-01-01",
end_date = "2020-01-01")
event_line(heatwave,
min_duration = 5,
metric = "intensity_mean",
category = TRUE,
spread = 150,
start_date = "2000-01-01",
end_date = "2020-01-01")
#Detect consecutive days in exceedance of a given threshold.
EX_heatwave <- exceedance(data = tempData,threshold = 35,minDuration = 3,maxGap = 2)
EX_heatwave
filename <- "E:/Desktop/ili_global_EX_heatwave.csv"
write.csv(EX_heatwave$threshold,file = filename)
