# Prepare tools and folder paths for the climate-data analysis.
library(terra)
library(dplyr)
library(purrr)
library(stringr)
library(readr)
library(ggplot2)
raw_dir <- "01_era5_land_original_downloads"
processed_dir <- "02_processed_data"
figures_dir <- "04_figures"
dir.create(processed_dir, showWarnings = FALSE, recursive = TRUE)
dir.create(figures_dir, showWarnings = FALSE, recursive = TRUE)Guangzhou ERA5-Land Climate Data Analysis
Supporting climate-data appendix for HLAC705 dissertation
Figure 1.1
Annual number of JJA days with Tmax ≥35°C (any-grid event), Guangzhou, 1971–2025.
# Read annual JJA Tmax files and calculate any-grid days at or above 35°C.
tmax_files <- list.files(
raw_dir,
pattern = "tmax.*\\.nc$",
full.names = TRUE
) |>
sort()
stopifnot(length(tmax_files) == 55)
process_tmax_file <- function(file) {
r <- terra::rast(file)
year <- str_extract(basename(file), "\\d{4}") |>
as.integer()
r_celsius <- r - 273.15
dates <- terra::time(r)
if (all(is.na(dates))) {
dates <- seq(
from = as.Date(paste0(year, "-06-01")),
to = as.Date(paste0(year, "-08-31")),
by = "day"
)
} else {
dates <- as.Date(dates)
}
daily_any_grid_tmax_c <- terra::global(
r_celsius,
fun = "max",
na.rm = TRUE
)[, 1] |>
as.numeric()
tibble(
year = year,
date = dates,
daily_any_grid_tmax_c = daily_any_grid_tmax_c,
tmax_ge35_any_grid = daily_any_grid_tmax_c >= 35
)
}
daily_tmax_any_grid <- map_dfr(tmax_files, process_tmax_file)
annual_tmax_ge35 <- daily_tmax_any_grid |>
group_by(year) |>
summarise(
jja_days_ge35_any_grid = sum(tmax_ge35_any_grid, na.rm = TRUE),
max_jja_any_grid_tmax_c = max(daily_any_grid_tmax_c, na.rm = TRUE),
.groups = "drop"
) |>
arrange(year) |>
mutate(
ge35_days_5yr_ma = as.numeric(
stats::filter(
jja_days_ge35_any_grid,
rep(1 / 5, 5),
sides = 2
)
)
)
annual_tmax_ge35# A tibble: 55 × 4
year jja_days_ge35_any_grid max_jja_any_grid_tmax_c ge35_days_5yr_ma
<int> <int> <dbl> <dbl>
1 1971 2 36.1 NA
2 1972 0 34.8 NA
3 1973 0 33.6 0.4
4 1974 0 34.4 0
5 1975 0 33.9 0
6 1976 0 34.1 0
7 1977 0 34.3 0.2
8 1978 0 34.7 0.4
9 1979 1 35.2 0.4
10 1980 1 36.5 0.6
# ℹ 45 more rows# Save the daily and annual Tmax results as processed CSV files.
write_csv(
daily_tmax_any_grid,
file.path(processed_dir, "daily_tmax_any_grid_guangzhou_jja_1971_2025.csv")
)
write_csv(
annual_tmax_ge35,
file.path(processed_dir, "annual_jja_days_tmax_ge35_any_grid_guangzhou_1971_2025.csv")
)# Plot annual JJA days at or above 35°C and a centred 5-year moving average.
fig_1_1 <- ggplot(annual_tmax_ge35, aes(x = year)) +
annotate(
"rect",
xmin = 1970.5, xmax = 1990.5,
ymin = -Inf, ymax = Inf,
fill = "#fff7bc", alpha = 0.1
) +
annotate(
"rect",
xmin = 1990.5, xmax = 2010.5,
ymin = -Inf, ymax = Inf,
fill = "#fec44f", alpha = 0.1
) +
annotate(
"rect",
xmin = 2010.5, xmax = 2025.5,
ymin = -Inf, ymax = Inf,
fill = "#d95f0e", alpha = 0.1
) +
annotate(
"text",
x = 1980.5, y = 11.8,
label = "1971–1990",
color = "grey40",
size = 4
) +
annotate(
"text",
x = 2000.5, y = 11.8,
label = "1991–2010",
color = "grey40",
size = 4
) +
annotate(
"text",
x = 2018, y = 11.8,
label = "2011–2025",
color = "grey40",
size = 4
) +
geom_col(
aes(
y = jja_days_ge35_any_grid,
fill = "Annual count"
),
alpha = 0.8
) +
geom_line(
data = annual_tmax_ge35 |> filter(!is.na(ge35_days_5yr_ma)),
aes(
y = ge35_days_5yr_ma,
color = "Centred 5-year moving average"
),
linewidth = 1
) +
scale_fill_manual(
values = c("Annual count" = "grey60")
) +
scale_color_manual(
values = c("Centred 5-year moving average" = "#CD661D")
) +
labs(
title = "Annual number of JJA days with Tmax ≥35°C in Guangzhou, 1971–2025",
subtitle = "Any-grid event based on daily maximum 2 m temperature",
x = "Year",
y = "Number of JJA days with Tmax ≥35°C",
caption = "JJA refers to June–August.\nA day is counted when at least one Guangzhou grid cell reaches daily maximum 2 m temperature ≥35°C."
) +
scale_x_continuous(
breaks = seq(1975, 2025, by = 5)
) +
coord_cartesian(
ylim = c(0, 12)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
plot.caption = element_text(size = 9, hjust = 0),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom",
legend.title = element_blank()
)
fig_1_1
# Save Figure 1.1 as a high-resolution PNG file.
ggsave(
filename = file.path(
figures_dir,
"figure_1_1_jja_days_tmax_ge35_any_grid_guangzhou_1971_2025.png"
),
plot = fig_1_1,
width = 8,
height = 5,
dpi = 300
)Figure 1.2
Mean JJA dewpoint (°C), Guangzhou, 1971–2025.
# Calculate annual mean JJA dewpoint for Guangzhou.
dewpoint_files <- list.files(
raw_dir,
pattern = "dewpoint.*\\.nc$",
full.names = TRUE
) |>
sort()
stopifnot(length(dewpoint_files) == 55)
process_dewpoint_file <- function(file) {
r <- terra::rast(file)
year <- str_extract(basename(file), "\\d{4}") |>
as.integer()
r_celsius <- r - 273.15
dates <- terra::time(r)
if (all(is.na(dates))) {
dates <- seq(
from = as.Date(paste0(year, "-06-01")),
to = as.Date(paste0(year, "-08-31")),
by = "day"
)
} else {
dates <- as.Date(dates)
}
daily_mean_dewpoint_c <- terra::global(
r_celsius,
fun = "mean",
na.rm = TRUE
)[, 1] |>
as.numeric()
tibble(
year = year,
date = dates,
daily_mean_dewpoint_c = daily_mean_dewpoint_c
)
}
daily_dewpoint <- map_dfr(
dewpoint_files,
process_dewpoint_file
)
annual_dewpoint <- daily_dewpoint |>
group_by(year) |>
summarise(
mean_jja_dewpoint_c = mean(daily_mean_dewpoint_c, na.rm = TRUE),
.groups = "drop"
) |>
arrange(year) |>
mutate(
dewpoint_5yr_ma = as.numeric(
stats::filter(
mean_jja_dewpoint_c,
rep(1 / 5, 5),
sides = 2
)
)
)
annual_dewpoint# A tibble: 55 × 3
year mean_jja_dewpoint_c dewpoint_5yr_ma
<int> <dbl> <dbl>
1 1971 24.2 NA
2 1972 24.1 NA
3 1973 24.2 24.1
4 1974 23.9 24.1
5 1975 24.3 24.2
6 1976 23.9 24.2
7 1977 24.6 24.3
8 1978 24.5 24.4
9 1979 24.3 24.5
10 1980 24.6 24.4
# ℹ 45 more rows# Save the daily and annual dewpoint results as processed CSV files.
write_csv(
daily_dewpoint,
file.path(processed_dir, "daily_mean_dewpoint_guangzhou_jja_1971_2025.csv")
)
write_csv(
annual_dewpoint,
file.path(processed_dir, "annual_mean_jja_dewpoint_guangzhou_1971_2025.csv")
)# Plot mean JJA dewpoint and a centred 5-year moving average.
fig_1_2 <- ggplot(annual_dewpoint, aes(x = year)) +
annotate(
"rect",
xmin = 1971.5, xmax = 1990.5,
ymin = -Inf, ymax = Inf,
fill = "#deebf7", alpha = 0.2
) +
annotate(
"rect",
xmin = 1990.5, xmax = 2010.5,
ymin = -Inf, ymax = Inf,
fill = "#9ecae1", alpha = 0.2
) +
annotate(
"rect",
xmin = 2010.5, xmax = 2025.5,
ymin = -Inf, ymax = Inf,
fill = "#3182bd", alpha = 0.2
) +
annotate(
"text",
x = 1980.5, y = 25.8,
label = "1971–1990",
color = "grey40",
size = 4
) +
annotate(
"text",
x = 2000.5, y = 25.8,
label = "1991–2010",
color = "grey40",
size = 4
) +
annotate(
"text",
x = 2018, y = 25.8,
label = "2011–2025",
color = "grey40",
size = 4
) +
geom_hline(
yintercept = c(24, 25),
linetype = "dashed",
color = "grey40",
linewidth = 0.4
) +
geom_line(
aes(
y = mean_jja_dewpoint_c,
color = "Annual mean"
),
linewidth = 0.5,
alpha = 0.7
) +
geom_point(
aes(
y = mean_jja_dewpoint_c,
color = "Annual mean"
),
size = 1,
alpha = 0.8
) +
geom_line(
data = annual_dewpoint |> filter(!is.na(dewpoint_5yr_ma)),
aes(
y = dewpoint_5yr_ma,
color = "Centred 5-year moving average"
),
linewidth = 1.1
) +
scale_color_manual(
values = c(
"Annual mean" = "grey40",
"Centred 5-year moving average" = "#212B58"
)
) +
labs(
title = "Mean JJA dewpoint in Guangzhou, 1971–2025",
subtitle = "Based on daily mean 2 m dewpoint temperature",
x = "Year",
y = "Mean JJA dewpoint (°C)",
caption = "JJA refers to June–August.\nDashed lines mark 24°C and 25°C reference levels for humid-heat conditions."
) +
scale_x_continuous(
breaks = seq(1975, 2025, by = 5)
) +
coord_cartesian(
ylim = c(23, 26)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
plot.caption = element_text(size = 9, hjust = 0),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom",
legend.title = element_blank()
)
fig_1_2
# Save Figure 1.2 as a high-resolution PNG file.
ggsave(
filename = file.path(
figures_dir,
"figure_1_2_mean_jja_dewpoint_guangzhou_1971_2025.png"
),
plot = fig_1_2,
width = 8,
height = 5,
dpi = 300
)Figure 4.1
Average number of JJA wet days and heavy-precipitation days, Guangzhou, 2010–2025.
# Calculate annual JJA wet days and heavy-precipitation days for Guangzhou.
precip_files <- list.files(
raw_dir,
pattern = "total_precipitation.*\\.nc$",
full.names = TRUE
) |>
sort()
process_precip_file <- function(file) {
r <- terra::rast(file)
year <- str_extract(basename(file), "\\d{4}") |>
as.integer()
month <- str_extract(basename(file), "_\\d{2}\\.nc$") |>
str_remove_all("_|\\.nc") |>
as.integer()
# ERA5-Land total precipitation is stored in metres, so convert to millimetres.
r_mm <- r * 1000
dates <- terra::time(r)
if (all(is.na(dates))) {
first_day <- as.Date(sprintf("%04d-%02d-01", year, month))
last_day <- seq(first_day, by = "month", length.out = 2)[2] - 1
dates <- seq(first_day, last_day, by = "day")
} else {
dates <- as.Date(dates)
}
daily_mean_precip_mm <- terra::global(
r_mm,
fun = "mean",
na.rm = TRUE
)[, 1] |>
as.numeric()
tibble(
year = year,
month = month,
date = dates,
daily_mean_precip_mm = daily_mean_precip_mm
)
}
daily_precipitation <- map_dfr(
precip_files,
process_precip_file
) |>
filter(
year >= 2010,
year <= 2025,
month %in% c(6, 7, 8)
)
annual_precipitation_days <- daily_precipitation |>
group_by(year) |>
summarise(
wet_days_ge1mm = sum(daily_mean_precip_mm >= 1, na.rm = TRUE),
heavy_precip_days_ge10mm = sum(daily_mean_precip_mm >= 10, na.rm = TRUE),
jja_days = n(),
.groups = "drop"
)
average_precipitation_days <- tibble(
precipitation_category = c(
"Wet days\n≥1 mm",
"Heavy-precipitation days\n≥10 mm"
),
average_days = c(
mean(annual_precipitation_days$wet_days_ge1mm, na.rm = TRUE),
mean(annual_precipitation_days$heavy_precip_days_ge10mm, na.rm = TRUE)
)
) |>
mutate(
percentage_of_jja = average_days / 92 * 100
)
average_precipitation_days# A tibble: 2 × 3
precipitation_category average_days percentage_of_jja
<chr> <dbl> <dbl>
1 "Wet days\n≥1 mm" 77.9 84.7
2 "Heavy-precipitation days\n≥10 mm" 31.2 33.9# Save the daily, annual and average precipitation results as processed CSV files.
write_csv(
daily_precipitation,
file.path(processed_dir, "daily_mean_precipitation_guangzhou_jja_2010_2025.csv")
)
write_csv(
annual_precipitation_days,
file.path(processed_dir, "annual_jja_precipitation_days_guangzhou_2010_2025.csv")
)
write_csv(
average_precipitation_days,
file.path(processed_dir, "average_jja_precipitation_days_guangzhou_2010_2025.csv")
)# Plot average JJA wet days and heavy-precipitation days.
plot_precipitation_days <- average_precipitation_days |>
mutate(
precipitation_category = factor(
precipitation_category,
levels = c(
"Wet days\n≥1 mm",
"Heavy-precipitation days\n≥10 mm"
)
)
)
fig_4_1 <- ggplot(
plot_precipitation_days,
aes(
x = precipitation_category,
y = average_days,
fill = precipitation_category
)
) +
geom_col(
width = 0.6,
alpha = 0.9
) +
geom_text(
aes(
label = paste0(
round(average_days, 1),
" days\n(",
round(percentage_of_jja, 1),
"%)"
)
),
vjust = -0.4,
size = 4
) +
scale_fill_manual(
values = c(
"Wet days\n≥1 mm" = "#9ecae1",
"Heavy-precipitation days\n≥10 mm" = "#212B58"
)
) +
labs(
title = "Average JJA wet and heavy-precipitation days in Guangzhou, 2010–2025",
subtitle = "Based on area-mean daily precipitation",
x = NULL,
y = "Average number of JJA days",
caption = "JJA refers to June–August.\nPercentages refer to the share of all JJA days."
) +
scale_y_continuous(
limits = c(0, 92),
breaks = seq(0, 90, by = 15)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold"),
plot.caption = element_text(size = 9, hjust = 0),
axis.text.x = element_text(size = 11),
legend.position = "none"
)
fig_4_1
# Save Figure 4.1 as a high-resolution PNG file.
ggsave(
filename = file.path(
figures_dir,
"figure_4_1_average_jja_precipitation_days_guangzhou_2010_2025.png"
),
plot = fig_4_1,
width = 8,
height = 5,
dpi = 300
)