scripts/ChinaHW_cluster/temp/Figure1_HW_characteristics_variation.R
In CUG-hydro/heatwave: What the Package Does (One Line, Title Case)
#待解决问题 增加slope pvalue ylab大小 调整-输出长度
library(ncdf4)
library(glue)
library(data.table)
library(ggplot2)
library(lubridate)
library(magrittr)
library(foreach)
library(iterators)
library(egg)
library(ggpmisc)
library(dplyr)
library(raster)
theme_new <- function(base_size = 30,
base_family = "",
base_line_size = base_size / 170,
base_rect_size = base_size / 170){
theme_void(base_size = base_size,
base_family = base_family,
base_line_size = base_line_size) %+replace%
theme(
axis.title.x = element_text(vjust = -0.5),
plot.title = element_text(
color = rgb(25, 43, 65, maxColorValue = 255),
face = "bold",
hjust = 0),
axis.title = element_text(
size = base_size,
color = '#2f2f2f'),
axis.text = element_text(
color = '#121b30',
size = rel(0.78)),
axis.text.y = element_text(hjust = 1),
panel.grid.major = element_line(
color = '#d0d0d0',
linetype = "dashed",
size = 0.5),
axis.title.y = element_blank(),
panel.border = element_rect(colour = "#5e6164", fill = NA, size = 1.5),
plot.margin = unit(c(1,1,1,1), "cm"),
panel.spacing.y = unit(1, 'cm'),
panel.spacing.x = unit(1, 'cm'),
axis.ticks.length = unit(.3, "cm"),
axis.ticks=element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size = rel(0.95))
)
}
label = data.frame(variable = c('D', 'I', 'S',
'area_mean', 'area_max', 'area_sum'),
tag = c('(a)~HWD', '(b)~HWI', '(c)~HWS',
'(d)~HWA[mean]', '(e)~HWA[max]', '(f)~HWA[sum]'),
stringsAsFactors = TRUE)
ylab_mean = c(`D` = 'Annual~Mean~Duration~(day)',
`I` = 'Annual~Mean~Intensity~(degree*C)',
`S` = 'Annual~Mean~Severity~(degree*C%.%day)',
`area_mean` = 'Annual~Mean~Area~(10^4~km^{2})',
`area_max` = 'Annual~Mean~Area~(10^4~km^{2})',
`area_sum` = 'Annual~Mean~Area~(10^4~km^{2})')
ylab_max = c(`D` = 'Annual~Max~Duration~(day)',
`I` = 'Annual~Max~Intensity~(degree*C)',
`S` = 'Annual~Max~Severity~(degree*C%.%day)',
`area_mean` = 'Annual~Max~Area~(10^6~km^{2})',
`area_max` = 'Annual~Max~Area~(10^6~km^{2})',
`area_sum` = 'Annual~Max~Area~(10^6~km^{2})')
lm_trend <- function(dat) {
fit <- lm(value ~ year, data = dat)
info <- summary(fit)
slope <- sprintf("slope == %.2f * ~decade ^ {-1}",
round(info$coefficients[2, 1] * 10, 2))
if(info$coefficients[2, 4] < 0.001) {
pvalue <- '"p-value < 0.001"'
} else {
pvalue <- sprintf('"p-value" == %.3f', round(info$coefficients[2, 4], 3))
}
data.frame(slope = slope, pvalue = pvalue,
y_t = 'if'(round(info$coefficients[2, 1] * 10, 2) >= 0, 0.11, 0.9),
y_p = 'if'(round(info$coefficients[2, 1] * 10, 2) >= 0, 0.03, 0.80),
stringsAsFactors = FALSE)
}
## HI
file = "OUTPUT/clusterId/HI/HW1.nc"
nc = nc_open(file)
date <- ncvar_get(nc, 'time') %>% as_date(origin = '1970-01-01')
a <- raster(file) %>% area
load('OUTPUT/HI/tmin_arr_HI.rda')
load('OUTPUT/HI/tmax_arr_HI.rda')
load('OUTPUT/HI/tmean_arr_HI.rda')
library(nctools)
lon = ncread(file, "lon")
lat = ncread(file, "lon")
dims <- ncdim_def_lonlat(lon, lat, date)
ncwrite(list(HI_Tmax = tmax_arr_HI), "OUTPUT/Observed_HI_tmax-196101_201704.nc", dims = dims)
foreach(i = c(1:13), icount()) %do% {
if(i %in% c(1, 2, 5, 6)) {
raw <- tmean_arr_HI
} else if (i %in% c(3, 7, 8, 10, 11)) {
raw <- tmax_arr_HI
} else if (i %in% c(4, 9, 12, 13)) {
raw <- tmin_arr_HI
}
dim(raw) <- 283 * 163 * 20574
load(paste0('OUTPUT/diff/HI/HW', i, '.rda'))
cmd = glue("diff <- HW{i}")
eval(parse(text = cmd))
dim(diff) <- 283 * 163 * 20574
rm(list = paste0('HW', i))
gc()
clusterId <- nc_open(paste0('OUTPUT/clusterId/HI/HW', i, '.nc')) %>%
ncvar_get('clusterId')
dim(clusterId) <- 283 * 163 * 20574
dt <- data.table(raw, clusterId, diff,
date = rep(date, each = 283 * 163), area = values(a))[!is.na(clusterId)]
#plan1 yearly average
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = mean(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = mean(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = mean(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e4), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = mean(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e4), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = mean(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e4), .(clusterId, year(date))
][, .(area_sum = mean(area_sum)), year]$area_sum
draw <- draw[year <= 2016] %>% melt(id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
scale_x_continuous(breaks = c(1961, 1972, 1983, 1994, 2005, 2016)) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_mean, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/HI/HW', i, '_yearly_mean.png'),
plot, width = 17, height = 14)
#plan2 yearly max
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = max(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = max(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = max(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e6), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = max(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e6), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = max(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e6), .(clusterId, year(date))
][, .(area_sum = max(area_sum)), year]$area_sum
draw <- draw[year <= 2016] %>% melt(id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
scale_x_continuous(breaks = c(1961, 1972, 1983, 1994, 2005, 2016)) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_max, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/HI/HW', i, '_yearly_max.png'),
plot, width = 17, height = 14)
}
rm(tmax_arr_HI, tmin_arr_HI, tmean_arr_HI)
gc()
## nonHI
file = "OUTPUT/clusterId/nonHI/HW1.nc"
date <- nc_open(file) %>% ncvar_get('time') %>% as_date(origin = '1970-01-01')
file = "data-raw/CN05.1_Tm_1961_2018_daily_025x025.nc"
raw_mean <- nc_open(file) %>%
ncvar_get('tm')
dim(raw_mean) <- 283 * 163 * 21184
file = "data-raw/CN05.1_Tmax_1961_2018_daily_025x025.nc"
raw_max <- nc_open(file) %>%
ncvar_get('tmax')
dim(raw_max) <- 283 * 163 * 21184
file = "data-raw/CN05.1_Tmin_1961_2018_daily_025x025.nc"
raw_min <- nc_open(file) %>%
ncvar_get('tmin')
dim(raw_min) <- 283 * 163 * 21184
foreach(i = c(1:13), icount()) %do% {
if(i %in% c(1, 2, 5, 6)) {
raw <- raw_mean
} else if (i %in% c(3, 7, 8, 10, 11)) {
raw <- raw_max
} else if (i %in% c(4, 9, 12, 13)) {
raw <- raw_min
}
dim(raw) <- 283 * 163 * 21184
load(paste0('OUTPUT/diff/nonHI/HW', i, '.nc'))
cmd = glue("diff <- HW{i}")
eval(parse(text = cmd))
dim(diff) <- 283 * 163 * 21184
rm(list = paste0('HW', i))
gc()
clusterId <- nc_open(paste0('OUTPUT/clusterId/nonHI/HW', i, '.nc')) %>%
ncvar_get('clusterId')
dim(clusterId) <- 283 * 163 * 21184
dt <- data.table(raw, clusterId, diff,
date = rep(date, each = 283 * 163), area = values(a))[!is.na(clusterId)]
#plan1 yearly average
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = mean(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = mean(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = mean(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e4), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = mean(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e4), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = mean(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e4), .(clusterId, year(date))
][, .(area_sum = mean(area_sum)), year]$area_sum
draw <- melt(draw, id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
xlim(1960, 2018) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_mean, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/nonHI/HW', i, '_yearly_mean.png'),
plot, width = 17, height = 14)
#plan2 yearly max
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = max(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = max(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = max(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e6), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = max(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e6), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = max(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e6), .(clusterId, year(date))
][, .(area_sum = max(area_sum)), year]$area_sum
draw <- melt(draw, id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
xlim(1960, 2018) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_max, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/nonHI/HW', i, '_yearly_max.png'),
plot, width = 17, height = 14)
}
CUG-hydro/heatwave documentation built on Dec. 17, 2021, 1:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#待解决问题 增加slope pvalue ylab大小 调整-输出长度
library(ncdf4)
library(glue)
library(data.table)
library(ggplot2)
library(lubridate)
library(magrittr)
library(foreach)
library(iterators)
library(egg)
library(ggpmisc)
library(dplyr)
library(raster)
theme_new <- function(base_size = 30,
base_family = "",
base_line_size = base_size / 170,
base_rect_size = base_size / 170){
theme_void(base_size = base_size,
base_family = base_family,
base_line_size = base_line_size) %+replace%
theme(
axis.title.x = element_text(vjust = -0.5),
plot.title = element_text(
color = rgb(25, 43, 65, maxColorValue = 255),
face = "bold",
hjust = 0),
axis.title = element_text(
size = base_size,
color = '#2f2f2f'),
axis.text = element_text(
color = '#121b30',
size = rel(0.78)),
axis.text.y = element_text(hjust = 1),
panel.grid.major = element_line(
color = '#d0d0d0',
linetype = "dashed",
size = 0.5),
axis.title.y = element_blank(),
panel.border = element_rect(colour = "#5e6164", fill = NA, size = 1.5),
plot.margin = unit(c(1,1,1,1), "cm"),
panel.spacing.y = unit(1, 'cm'),
panel.spacing.x = unit(1, 'cm'),
axis.ticks.length = unit(.3, "cm"),
axis.ticks=element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size = rel(0.95))
)
}
label = data.frame(variable = c('D', 'I', 'S',
'area_mean', 'area_max', 'area_sum'),
tag = c('(a)~HWD', '(b)~HWI', '(c)~HWS',
'(d)~HWA[mean]', '(e)~HWA[max]', '(f)~HWA[sum]'),
stringsAsFactors = TRUE)
ylab_mean = c(`D` = 'Annual~Mean~Duration~(day)',
`I` = 'Annual~Mean~Intensity~(degree*C)',
`S` = 'Annual~Mean~Severity~(degree*C%.%day)',
`area_mean` = 'Annual~Mean~Area~(10^4~km^{2})',
`area_max` = 'Annual~Mean~Area~(10^4~km^{2})',
`area_sum` = 'Annual~Mean~Area~(10^4~km^{2})')
ylab_max = c(`D` = 'Annual~Max~Duration~(day)',
`I` = 'Annual~Max~Intensity~(degree*C)',
`S` = 'Annual~Max~Severity~(degree*C%.%day)',
`area_mean` = 'Annual~Max~Area~(10^6~km^{2})',
`area_max` = 'Annual~Max~Area~(10^6~km^{2})',
`area_sum` = 'Annual~Max~Area~(10^6~km^{2})')
lm_trend <- function(dat) {
fit <- lm(value ~ year, data = dat)
info <- summary(fit)
slope <- sprintf("slope == %.2f * ~decade ^ {-1}",
round(info$coefficients[2, 1] * 10, 2))
if(info$coefficients[2, 4] < 0.001) {
pvalue <- '"p-value < 0.001"'
} else {
pvalue <- sprintf('"p-value" == %.3f', round(info$coefficients[2, 4], 3))
}
data.frame(slope = slope, pvalue = pvalue,
y_t = 'if'(round(info$coefficients[2, 1] * 10, 2) >= 0, 0.11, 0.9),
y_p = 'if'(round(info$coefficients[2, 1] * 10, 2) >= 0, 0.03, 0.80),
stringsAsFactors = FALSE)
}
## HI
file = "OUTPUT/clusterId/HI/HW1.nc"
nc = nc_open(file)
date <- ncvar_get(nc, 'time') %>% as_date(origin = '1970-01-01')
a <- raster(file) %>% area
load('OUTPUT/HI/tmin_arr_HI.rda')
load('OUTPUT/HI/tmax_arr_HI.rda')
load('OUTPUT/HI/tmean_arr_HI.rda')
library(nctools)
lon = ncread(file, "lon")
lat = ncread(file, "lon")
dims <- ncdim_def_lonlat(lon, lat, date)
ncwrite(list(HI_Tmax = tmax_arr_HI), "OUTPUT/Observed_HI_tmax-196101_201704.nc", dims = dims)
foreach(i = c(1:13), icount()) %do% {
if(i %in% c(1, 2, 5, 6)) {
raw <- tmean_arr_HI
} else if (i %in% c(3, 7, 8, 10, 11)) {
raw <- tmax_arr_HI
} else if (i %in% c(4, 9, 12, 13)) {
raw <- tmin_arr_HI
}
dim(raw) <- 283 * 163 * 20574
load(paste0('OUTPUT/diff/HI/HW', i, '.rda'))
cmd = glue("diff <- HW{i}")
eval(parse(text = cmd))
dim(diff) <- 283 * 163 * 20574
rm(list = paste0('HW', i))
gc()
clusterId <- nc_open(paste0('OUTPUT/clusterId/HI/HW', i, '.nc')) %>%
ncvar_get('clusterId')
dim(clusterId) <- 283 * 163 * 20574
dt <- data.table(raw, clusterId, diff,
date = rep(date, each = 283 * 163), area = values(a))[!is.na(clusterId)]
#plan1 yearly average
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = mean(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = mean(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = mean(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e4), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = mean(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e4), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = mean(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e4), .(clusterId, year(date))
][, .(area_sum = mean(area_sum)), year]$area_sum
draw <- draw[year <= 2016] %>% melt(id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
scale_x_continuous(breaks = c(1961, 1972, 1983, 1994, 2005, 2016)) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_mean, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/HI/HW', i, '_yearly_mean.png'),
plot, width = 17, height = 14)
#plan2 yearly max
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = max(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = max(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = max(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e6), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = max(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e6), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = max(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e6), .(clusterId, year(date))
][, .(area_sum = max(area_sum)), year]$area_sum
draw <- draw[year <= 2016] %>% melt(id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
scale_x_continuous(breaks = c(1961, 1972, 1983, 1994, 2005, 2016)) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_max, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/HI/HW', i, '_yearly_max.png'),
plot, width = 17, height = 14)
}
rm(tmax_arr_HI, tmin_arr_HI, tmean_arr_HI)
gc()
## nonHI
file = "OUTPUT/clusterId/nonHI/HW1.nc"
date <- nc_open(file) %>% ncvar_get('time') %>% as_date(origin = '1970-01-01')
file = "data-raw/CN05.1_Tm_1961_2018_daily_025x025.nc"
raw_mean <- nc_open(file) %>%
ncvar_get('tm')
dim(raw_mean) <- 283 * 163 * 21184
file = "data-raw/CN05.1_Tmax_1961_2018_daily_025x025.nc"
raw_max <- nc_open(file) %>%
ncvar_get('tmax')
dim(raw_max) <- 283 * 163 * 21184
file = "data-raw/CN05.1_Tmin_1961_2018_daily_025x025.nc"
raw_min <- nc_open(file) %>%
ncvar_get('tmin')
dim(raw_min) <- 283 * 163 * 21184
foreach(i = c(1:13), icount()) %do% {
if(i %in% c(1, 2, 5, 6)) {
raw <- raw_mean
} else if (i %in% c(3, 7, 8, 10, 11)) {
raw <- raw_max
} else if (i %in% c(4, 9, 12, 13)) {
raw <- raw_min
}
dim(raw) <- 283 * 163 * 21184
load(paste0('OUTPUT/diff/nonHI/HW', i, '.nc'))
cmd = glue("diff <- HW{i}")
eval(parse(text = cmd))
dim(diff) <- 283 * 163 * 21184
rm(list = paste0('HW', i))
gc()
clusterId <- nc_open(paste0('OUTPUT/clusterId/nonHI/HW', i, '.nc')) %>%
ncvar_get('clusterId')
dim(clusterId) <- 283 * 163 * 21184
dt <- data.table(raw, clusterId, diff,
date = rep(date, each = 283 * 163), area = values(a))[!is.na(clusterId)]
#plan1 yearly average
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = mean(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = mean(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = mean(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e4), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = mean(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e4), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = mean(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e4), .(clusterId, year(date))
][, .(area_sum = mean(area_sum)), year]$area_sum
draw <- melt(draw, id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
xlim(1960, 2018) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_mean, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/nonHI/HW', i, '_yearly_mean.png'),
plot, width = 17, height = 14)
#plan2 yearly max
draw <- data.table(year = unique(year(dt$date)))
draw$D <- dt[, .(D = uniqueN(date)), .(clusterId, year(date))
][, .(D = max(D)), year]$D
draw$I <- dt[, .(I = weighted.mean(diff, area)), .(clusterId, date)
][, .(I = mean(I)), .(clusterId, year(date))
][, .(I = max(I)), year]$I
draw$S <- dt[, .(S = weighted.mean(diff, area)), .(clusterId, date)
][, .(S = sum(S)), .(clusterId, year(date))
][, .(S = max(S)), year]$S
draw$area_mean <- dt[, .(area_mean = sum(area) / 1e6), .(clusterId, date)
][, .(area_mean = mean(area_mean)), .(clusterId, year(date))
][, .(area_mean = max(area_mean)), year]$area_mean
draw$area_max <- dt[, .(area_max = sum(area) / 1e6), .(clusterId, date)
][, .(area_max = max(area_max)), .(clusterId, year(date))
][, .(area_max = max(area_max)), year]$area_max
draw$area_sum <- dt[, .(area_sum = sum(area) / 1e6), .(clusterId, year(date))
][, .(area_sum = max(area_sum)), year]$area_sum
draw <- melt(draw, id.vars = 'year')
labels <- draw %>%
group_by(variable) %>%
do(lm_trend(.))
plot <- ggplot(draw, aes(year, value)) +
geom_smooth(fill = '#d9ecf9', colour="#014ed6", size = 3, alpha = 0.5, method = 'lm') +
geom_line(colour = '#3b435b', alpha = 0.7, size = 1.3) +
theme_new(base_size = 22, base_family = 'Arial') +
xlim(1960, 2018) +
xlab('Year') +
geom_text_npc(data = labels, aes(label = slope, npcy = y_t), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = labels, aes(label = pvalue, npcy = y_p), parse = TRUE,
npcx = 0.62, size = 6, color = '#014ed6', hjust = 0) +
geom_text_npc(data = label, aes(label = tag), parse = TRUE, color = '#0d0e28',
npcx = 0.05, npcy = 0.95, size = 7) +
facet_wrap(~variable, scales="free", nrow = 3,
label = as_labeller(ylab_max, label_parsed), strip.position = 'left')
ggsave(paste0('Figures/nonHI/HW', i, '_yearly_max.png'),
plot, width = 17, height = 14)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.