In CUG-hydro/heatwave: What the Package Does (One Line, Title Case)

目前13种热浪指标可以分为如下几类：

热浪指标定义

按照输入变量可以划分为

• 基于Tmin的：夜间热浪
• 基于Tmax的：日间热浪
• 基于Tavg的：综合热浪

按照定义方法可以划分为

• 相对阈值法

  • At least k consecutive days >= T1

• 滑动平均相对阈值法

  • At least k consecutive days >= T1
  • 15-day moving sampling by year

• 双变量阈值法

  • Everyday > T2
  • Average > T1
  • At least k consecutive days > T1

• 绝对阈值法

  • At least k consecutive days ≥ T1

| name | index | description | |------ | --------------| --------------------------- | | HW1 | Tavg_90 | K = 2 T1 = 90th percentile | | HW2 | Tavg_95 | K = 2 T1 = 95th percentile | | HW3 | Tmax_95 | K = 2 T1 = 95th percentile | | HW4 | Tmin_95 | K = 2 T1 = 95th percentile | | HW5 | Tavg_98 | K = 2 T1 = 98th percentile | | HW6 | Tavg_99 | K = 2 T1 = 99th percentile | | HW7 | Tmax_cma | K = 3 T1 = 35℃ | | HW8 | Tmax_mov_95 | K = 3 T1 = 95th percentile | | HW9 | Tmin_mov_95 | K = 3 T1 = 95th percentile | | HW10 | Tmax_2con_975 | K = 3 T1 = 97.5th percentile T2 = 81th percentile | | HW11 | Tmax_2con_90 | K = 3 T1 = 90th percentile T2 = 75th percentile | | HW12 | Tmin_2con_975 | K = 3 T1 = 97.5th percentile T2 = 81th percentile | | HW13 | Tmin_2con_90 | K = 3 T1 = 90th percentile T2 = 75th percentile |

source('scripts/main_pkgs.R')

TRS = c(75, 81, 90, 95, 97.5, 98, 99, 99.9, 99.99)

热浪的变化趋势

首先把体感温度heat index转存为nc格式数据。

file = "OUTPUT/clusterId/HI/HW1.nc"
date <- nc_date(file)
lon = ncread(file, "lon")
lat = ncread(file, "lat")
dims <- ncdim_def_lonlat(lon, lat, date)

load("OUTPUT/HI/tmax_arr_HI.rda")
ncwrite(list(HI_Tmax = tmax_arr_HI)  , "OUTPUT/Observed_HI_tmax-196101_201704.nc", dims = dims)

load("OUTPUT/HI/tmin_arr_HI.rda")
ncwrite(list(HI_Tmin = tmin_arr_HI), "OUTPUT/Observed_HI_tmin-196101_201704.nc", dims = dims)

load("OUTPUT/HI/tmean_arr_HI.rda")
ncwrite(list(HI_Tmean = tmean_arr_HI), "OUTPUT/Observed_HI_tmean-196101_201704.nc", dims = dims)
# InitCluster(4)
# llply_par()

files_HI = dir("OUTPUT/", "*.nc", full.names = TRUE)
files_noHI = dir("data-raw", "CN05.1_T", full.names = TRUE)

# mark NA values
file = files_HI[1]
lat = ncread(file, "lat")
lon = ncread(file, "lon")

dates = nc_date(file)
inds = which(year(dates) <= 2016)
dates = dates[inds]
arr <- ncread(files[1])$data[[1]][,,inds]

if (0) {
    arr_year = apply_3d(arr, by = year(dates))
    arr_multiYearMean = apply_3d(arr_year)
    inds_zh = which.notna(arr_multiYearMean)
    # arr = arr[inds_zh]
    # 仅保留有值的部分
    # load("OUTPUT/HI/tmax_arr_HI.rda"

    ## 保存中国数据范围至data
    # if  
    grid_d025 = get_grid.lonlat(lon, lat) 
    grid_d025$id = NA
    grid_d025$id[inds_zh] = 1
    spplot(grid_d025)
    use_data(grid_d025, overwrite = TRUE)
}
# save(grid_d025, inds_d025)
I_grid = grid_d025$id %>% which.notna() 

devtools::load_all()
files_HI = dir("OUTPUT/", "*.nc", full.names = TRUE)
files_noHI = dir("data-raw", "CN05.1_T", full.names = TRUE)

I_grid = grid_d025$id %>% which.notna() 
files = c(files_HI, files_noHI)

InitCluster(6, kill = FALSE)
l_TRS = foreach(file = files) %dopar% {
    l <- read_nc(file)
    cal_threshold(l$arr, l$dates, NULL, I_grid = I_grid)    
}
names <- c("HI_Tmax", "HI_Tavg", "HI_Tmin", "noHI_Tavg", "noHI_Tmax", "noHI_Tmin")
l_TRS = set_names(l_TRS, names)
l_TRS = l_TRS[c("HI_Tavg", "HI_Tmax", "HI_Tmin", "noHI_Tavg", "noHI_Tmax", "noHI_Tmin")]
# check l_TRS
save(l_TRS, file = "OUTPUT/Observed_HI&noHI_TRS.rda")
load("OUTPUT/Observed_HI&noHI_TRS.rda")
map(l_TRS, ~.[,150,1]) %>% as.data.table() %>% na.omit()

dates = seq(as.Date('2010-01-01'), as.Date('2013-12-31'), "day")
set.seed(1)
x = rnorm(length(dates)) %>% as.matrix()

CUG-hydro/heatwave documentation built on Dec. 17, 2021, 1:53 p.m.