利用众包气温数据和LCZ4r分析城市热岛效应
Max Anjos
April 08, 2026
Source:vignettes/local_func_modeling_crows.Rmd
local_func_modeling_crows.Rmd介绍
城市热岛 (UHI) 是城市中的一个严重环境问题,城市地区的气温高于农村地区。本案例研究演示了如何使用来自公民气象站 (CWS) 的众包气温 (Tair) 数据和 R 中的 LCZ4r 软件包来分析 UHI。我们重点关注 2023 年 5 月的每小时 Tair 数据,由德国柏林的 Netatmo CWS 收集,可通过 R 中的 CrowdQCplus 包.
# 安装并加载必要的软件包
if (!require("pacman")) install.packages("pacman")
install.packages("sf", type = "binary")
pacman::p_load(remotes, dplyr, tmap, ggplot2, devtools, data.table)
library(dplyr) # 用于数据处理
library(sf) # 用于向量数据处理
library(tmap) # 用于交互式地图可视化
library(ggplot2) # 用于数据可视化
library(data.table)
# 安装 CrowdQCplus 以进行 CWS 数据质量控制
if (!require("CrowdQCplus")) {
remotes::install_github("dafenner/CrowdQCplus")
}
library(CrowdQCplus)为什么要众包数据?
众包气象站提供前所未有的城市温度数据空间覆盖: - 高密度:全市有数百或数千个车站 - 实时数据:每小时或每小时一次的测量 - 具有成本效益:利用现有的公民基础设施 - 补充:填补官方监测网络的空白
然而,仔细的质量控制对于确保数据可靠性至关重要。
数据集准备
加载柏林 LCZ 地图
# 使用 LCZ 生成器平台获取柏林的 LCZ 地图
lcz_map <- LCZ4r::lcz_get_map_generator(ID = "8576bde60bfe774e335190f2e8fdd125dd9f4299")
# 可视化LCZ图
LCZ4r::lcz_plot_map(lcz_map)
柏林 LCZ 地图显示当地气候带的空间分布。该地图用作分析众包温度数据的空间框架。
CWS 数据的质量控制
为了确保数据的可靠性,我们采用了由 芬纳等人。 (2021) 使用 CrowdQCplus 包。此过程可以过滤掉统计上不可信的读数,从而提高数据准确性。
QC流程步骤:
- 数据填充:确保时间序列一致
- 输入验证:检查数据结构
- 质量控制:应用统计过滤器来识别异常值
- 输出统计数据:总结 QC 结果
# 执行数据填充和输入检查
data_crows <- cqcp_padding(cqcp_cws_data)
data_check <- cqcp_check_input(data_crows)
# 应用质量控制
if (data_check) {
data_qc <- cqcp_qcCWS(data_crows) # QC process
}
# 查看质量控制输出
head(data_qc)
# 显示 QC 输出统计信息
n_data_qc <- cqcp_output_statistics(data_qc)
print(n_data_qc)了解 QC 输出:
应用 QC 后,会生成多个列: - ta_int:通过所有质量控制步骤的温度值 - qc_flag:指示数据质量状态的质量控制标志 - qa_flag:用于额外验证的质量保证标志
对于气温插值,我们使用 ta_int
列,代表最可靠的温度估计。
LCZ 的昼夜周期分析
我们检查了不同 LCZ 特定日期(2023 年 5 月 1 日)气温的昼夜循环,以了解城市形态如何影响每日温度模式。
# 分析LCZ各类别中的日温度变化
cws_ts <- lcz_ts(lcz_map,
data_frame = data_qc,
var = "ta_int",
station_id = "p_id",
time.freq = "hour",
extract.method = "two.step",
year = 2023, month = 5, day = 1,
by = "daylight")
# 显示时间序列图
cws_ts
2023 年 5 月 1 日,柏林不同 LCZ 等级的昼夜温度循环。请注意,与植被和开放区域相比,紧凑的建筑区域全天保持更高的温度。
城市热岛分析
UHI 强度时间序列
我们使用 LCZ 方法计算了 2023 年 5 月每小时的 UHI 强度,该方法根据 LCZ 类别自动选择城市和农村温度。
# 计算2023年5月的城市热岛强度
cws_uhi <- lcz_uhi_intensity(lcz_map,
data_frame = data_qc,
var = "ta_int",
station_id = "p_id",
time.freq = "hour",
extract.method = "two.step",
year = 2023, month = 5,
method = "LCZ")
# 显示城市热岛强度时间序列
cws_uhi
2023年5月柏林的UHI强度时间序列。观测到的最高UHI强度为5月13日16:00 3.7°C,表明城市变暖效应显着。
城市热岛空间测绘
我们使用基于 LCZ 的空间插值对 2023 年 5 月 13 日 16:00(强度最大的时间)的 UHI 进行了建模。
温度图
# 绘制城市热岛效应高峰期的气温图
my_interp_map <- lcz_interp_map(
lcz_map,
data_frame = data_qc,
var = "ta_int",
station_id = "p_id",
sp.res = 100,
tp.res = "hour",
year = 2023, month = 5, day = 13, hour = 16
)
# 使用标题和标签自定义图表
lcz_plot_interp(
my_interp_map,
title = "Thermal Field - CWS Data",
subtitle = "Berlin - May 13, 2023 at 16:00",
caption = "Source: LCZ4r, 2024",
fill = "Temperature [°C]"
)
2023年5月13日16:00城市热岛高峰事件期间气温空间分布。该地图揭示了不同的热模式,紧凑的城市地区温度较高,植被和周边地区温度较低。
热异常图
# 针对同一事件生成热异常图
my_anomaly_map <- lcz_anomaly_map(
lcz_map,
data_frame = data_qc,
var = "ta_int",
station_id = "p_id",
sp.res = 100,
tp.res = "hour",
year = 2023, month = 5, day = 13, hour = 16
)
# 自定义异常图
lcz_plot_interp(
my_anomaly_map,
title = "Thermal Anomaly Field - CWS Data",
subtitle = "Berlin - May 13, 2023 at 16:00",
caption = "Source: LCZ4r, 2024",
fill = "Temperature Anomaly [°C]",
palette = "bl_yl_rd"
)
显示与城市平均水平偏差的热异常图。红色区域表示温度高于平均水平(热点),而蓝色区域表示温度低于平均水平(冷点)。
UHI地图精度评估
我们使用留一交叉验证 (LOOCV) 和计算指标(如 RMSE、MAE 和 sMAPE)评估插值精度。
交叉验证
# 评估插值精度
df_eval <- lcz_interp_eval(
lcz_map,
data_frame = data_qc,
var = "ta_int",
station_id = "p_id",
year = 2023, month = 5, day = 13,
LOOCV = TRUE,
extract.method = "two.step",
sp.res = 500,
tp.res = "hour",
vg.model = "Sph",
LCZinterp = TRUE
)LCZ 类别的评估指标
# 按LCZ类计算评估指标
df_eval_metrics <- df_eval %>%
group_by(lcz) %>%
summarise(
n_obs = n(), # Number of observations
rmse = sqrt(mean((observed - predicted)^2)), # RMSE
mae = mean(abs(observed - predicted)), # MAE
smape = mean(2 * abs(observed - predicted) /
(abs(observed) + abs(predicted)) * 100) # sMAPE
) %>%
arrange(rmse)
# 显示指标
df_eval_metrics相关性分析
# 相关图,包含回归方程和 R²
p1 <- ggplot(df_eval, aes(x = observed, y = predicted)) +
geom_point(alpha = 0.3, color = "#1D9E75", size = 1) +
geom_smooth(method = "lm", color = "red", se = FALSE, linewidth = 1) +
stat_density2d(aes(fill = after_stat(level)), geom = "polygon", alpha = 0.3) +
scale_fill_viridis_c(option = "viridis", name = "Density") +
stat_poly_eq(
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
formula = y ~ x,
parse = TRUE,
label.x.npc = "left",
label.y.npc = "top",
size = 4
) +
labs(
title = "Observed vs. Predicted Temperatures",
x = "Observed Temperature [°C]",
y = "Predicted Temperature [°C]"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
axis.title = element_text(face = "bold")
)
# 添加边际直方图
p1_with_marginals <- ggMarginal(p1,
type = "histogram",
fill = "#1D9E75",
bins = 30,
alpha = 0.6)
# 打印图表
p1_with_marginals
众包数据观测到的温度和预测的温度之间的相关性。高 R² 值表明模型性能良好,证明了 QC 后众包数据的可靠性。
按 LCZ 类别进行残差分析
# 按LCZ类别划分的残差图
p2 <- ggplot(df_eval, aes(x = observed, y = residual)) +
geom_point(alpha = 0.3, color = "darkgreen", size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 0.8) +
geom_smooth(method = "loess", color = "blue", se = FALSE, linewidth = 1) +
stat_density2d(aes(fill = after_stat(level)), geom = "polygon", alpha = 0.3) +
scale_fill_viridis_c(option = "plasma", direction = -1, name = "Density") +
facet_wrap(~ lcz, scales = "free", ncol = 3) +
labs(
title = "Residuals by LCZ Class",
x = "Observed Temperature [°C]",
y = "Residuals [°C]"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
axis.title = element_text(face = "bold"),
strip.text = element_text(size = 12, face = "bold"),
legend.position = "right",
panel.grid.major = element_line(color = "gray90", linetype = "dotted"),
panel.grid.minor = element_blank()
)
# 打印图表
p2
用于众包数据评估的 LCZ 类残差。零附近的随机散布表明模型在不同的城市形态中表现良好。
主要发现和讨论
保存结果
# 将评估指标保存到 CSV 文件
write.csv(df_eval_metrics,
file = "cws_uhi_metrics_may2023.csv",
row.names = FALSE)
# 保存完整评估数据框
write.csv(df_eval,
file = "cws_uhi_evaluation_may2023.csv",
row.names = FALSE)
# 保存地块
ggsave("cws_correlation_plot.png", p1_with_marginals, width = 10, height = 8, dpi = 300)
ggsave("cws_residuals_by_lcz.png", p2, width = 12, height = 10, dpi = 300)