scripts/ChinaHW_cluster/s1_热浪指标年际变化.R

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

# 逐年的统计结果
library(glue)
library(sf2)
library(plyr)
library(foreach)
library(iterators)
## 考虑网格的面积加权
raster2array <- function(r) {
    as.matrix(r) %>%t() %>% flipud()
}

{
    cellsize = 0.25
    grid = make_grid(range = c(70, 140, 15, 55) + c(-1, 1, -1, 1)*cellsize,
                  cellsize, c(FALSE, FALSE)) # type = "mat"
    # grid = l$grid.origin %$% get_grid.lonlat(lon, lat) # version2
    grid@data <- data.frame(x = as.numeric(mat[,,150]))
    r = raster(grid )
    area = raster2array(area(r) / 1e3) %>% as.numeric()
    grid$area = area
    ## version2: make_grid
    # delta = cellsize
    # r@data <- data.frame(x = as.numeric(mat[,,150]))
    # spplot(r)
}

## 需要把阈值单独导出来
file = "OUTPUT/clusterId/all/HW1.nc"
# 1961-2018
dates_full = nc_date(file)
year_max = max(dates_full) %>% year()
year_min = min(dates_full) %>% year()
# year = 2018
years = year_min:year_max %>% set_names(., .)

# res = foreach(year = years, i = icount()) %do% {
# l <- ncread(file, DatePeriod = date_period)
i = 0
res = list()
for (year in years) {
    i = i + 1
    runningId(i)

    date_period <- c(glue("{year}-01-01"), glue("{year}-12-31"))
    l <- ncread(file, DatePeriod = date_period)
    dates = as.Date(date_period) %>% {seq(.[1], .[2], by = "day")}
    arr <- l$data[[1]]
    inds = array_3dTo2d(arr) %>% apply(2, function(x) !all(is.na(x))) %>% which()

    # 统计id数
    info = table(as.numeric(arr)) %>% as.data.table() %>% set_colnames(c("id", "N"))
    info$id %<>% as.numeric()

    # 统计每个id的日变化
    d_HW = llply(info$id %>% set_names(., .), HW_statistics,
        arr = arr[, , inds], dates = dates[inds], area = grid$area
        # .progress = "text"
    ) %>% melt_list("id")

    # 统计热浪指标
    res[[i]] = d_HW[, .(
        n_day = .N,
        tol_pixels = sum(num),
        tol_area = sum(area),
        mean_pixels = mean(num),
        mean_area = mean(area),
        date_begin = min(date), date_end = max(date)
    ), .(id)] %>%
        cbind(I = 1:nrow(.), .)
}

df = melt_list(res) %>% mutate(year = year(date_begin))
fwrite(df, "HW1_statistics.csv")

library(tidyverse)
df[mean_area >= 10, .N, .(year)] %>%
    ggplot(aes(year, N)) +
    geom_point() + geom_line() +
    geom_smooth()

df[mean_area >= 100, .(value = sum(tol_area)), .(year)] %>%
    ggplot(aes(year, value)) +
    geom_point() + geom_line() +
    geom_smooth() +
    labs(y = "总热浪发生面积")

## 小热浪事件不在考虑范围之内 area >= 1e4 km, 1*1deg
# 热浪
# total_area
# intensity
# severity
# duration
df2 = df %>% group_by(year)

# ngrids <- dim(mat)[1:2] %>% prod()
# nc_info()