R语言 热浪识别（Heatwave）（ERA5数据、逐像元）
R语言初学者，欢迎大家交流学习
本代码基于R语言
逐像元识别重点使用的是 terra 包中的 app 函数
内部过程与全局计算相似，下面是全局计算的函数
#########################################
热浪识别重点使用的是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)
##################################################################
#######################################################################
#逐像元计算
#！！！！！！！！！！！！！！！！！！！浮动阈值（p=90）
heatDetect<- function(x,y){
if(length(na.omit(x))<10)
return(c(NA,NA,NA))
#数据准备
x <- x-273.15
tempdata <- as.data.frame(x)
tempdata$t <- y
names(tempdata) <- c("temp","t")
#热浪识别
heat_ts <- heatwaveR::ts2clm(tempdata,
climatologyPeriod = c("1963-03-08","1993-03-08"),
pctile = 90 ,
roundClm = 4 )
heatwave <- heatwaveR::detect_event(heat_ts,
maxGap = 2,
minDuration = 3,
categories = TRUE)
#输出准备
event_no <- heatwave$event_no[length(heatwave$event_no)]
duration_sum <- sum(heatwave$duration)
intensity_cumulative_mean <- mean(heatwave$intensity_cumulative,na.rm=T)
out <- c(event_no,duration_sum,intensity_cumulative_mean)
return(out)
}
heatdetect_app <- app(era_tem_2m,heatDetect,y=dataname,cores=4)
#命名lyres
names(heatdetect_app) <- c("event_no","duration_sum","intensity_cumulative_mean")
#展示
plot(heatdetect_app)
filename="E:/Desktop/Heatwave/heatdetect_app.tif"
writeRaster(heatdetect_app,filename = filename,overwrite=TRUE)
#### 固定阈值
EX_heatDetect<- function(x,y){
if(length(na.omit(x))<10)
return(c(NA,NA,NA))
#数据准备
x <- x-273.15
tempdata <- as.data.frame(x)
tempdata$t <- y
names(tempdata) <- c("temp","t")
#热浪识别
EX_heatwave <- heatwaveR::exceedance(tempdata,threshold = 35,minDuration = 3,maxGap = 2)
#输出准备
event_no <- EX_heatwave$exceedance$exceedance_no[length(EX_heatwave$exceedance$exceedance_no)]
duration_sum <- sum(EX_heatwave$exceedance$duration)
intensity_cumulative_mean <- mean(EX_heatwave$exceedance$intensity_cumulative,na.rm=T)
out <- c(event_no,duration_sum,intensity_cumulative_mean)
return(out)
}
heatDetect_EX_app <- app(era_tem_2m,EX_heatDetect,y = dataname,cores=4)
names(heatDetect_EX_app) <- c("event_no","duration_sum","intensity_cumulative_mean")
plot(heatDetect_EX_app)
