R語言初學(xué)者，歡迎大家交流學(xué)習(xí)

本代碼基于R語言?

重點(diǎn)的包為heatwaveR包

使用的氣溫數(shù)據(jù)為?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)

#加載矢量邊界數(shù)據(jù)

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")

#提取復(fù)合圖層中的溫度層

era_tem_2m <- rast()

##溫度數(shù)據(jù)

for (i in 1:length(dataPath)) {

? #lyrs=2是2米大氣溫度

? heatdata <- rast(dataPath[i],lyrs=2) %>%?

? ? ? ? ? ? ? crop(.,ili_vect) %>%?

? ? ? ? ? ? ? mask(.,ili_vect)

? #為每層數(shù)據(jù)命名

? names(heatdata) <- substr(dataPath[i],46,53)

? #將每層數(shù)據(jù)添加到數(shù)據(jù)集中

? add(era_tem_2m) <- heatdata

}

era_tem_2m

#保存結(jié)果

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(.)

#準(zhǔn)備熱浪識別數(shù)據(jù)

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)