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

本代碼基于R語言?

逐像元識別重點使用的是 terra 包中的 app 函數(shù)

內(nèi)部過程與全局計算相似，下面是全局計算的函數(shù)

R語言 熱浪（Heatwave）識別（ERA5數(shù)據(jù)）

#########################################

熱浪識別重點使用的是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)

##################################################################

#######################################################################

#逐像元計算

#?。。。。。。。。。。。。。。。。。?！浮動閾值（p=90）

heatDetect<- function(x,y){

? if(length(na.omit(x))<10)

? ? return(c(NA,NA,NA))

??

? #數(shù)據(jù)準備

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

??

? #數(shù)據(jù)準備

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