基于大型地理空间数据源计算人口加权热日数

原代码的核心性能瓶颈是嵌套循环逐个遍历栅格单元，每次循环还执行栅格转多边形、局部mask、全局求和等冗余操作，对于大型国家的栅格数据来说，这种方式时间复杂度极高，效率极低。以下是基于terra包矢量运算的优化方案，能将运算速度提升几个数量级：

优化思路 #

1. 移除所有循环操作，利用terra的栅格矢量运算特性直接批量计算
2. 消除冗余的栅格转多边形、局部mask操作，直接使用已裁剪对齐的栅格数据计算
3. 确保人口栅格与热日栅格的分辨率/范围完全对齐（若原数据未对齐，先做重采样）

优化后的完整代码 #

# Load input data (hot days – geospatial file – 1 deg. resolution)
hot_days <- rast('file_path.nc')
pop <- rast('file_path.tif')

# Import shape country 
shapefile <- st_read('country_path.shp')     
   
# Crop, reproject and mask raster by shapefile 
v <- vect(shapefile)        
clip1 <- crop(hot_days, ext(shapefile), snap='out') 
clip2 <- crop(pop, ext(shapefile))
mask1 <- mask(clip1, v) 
mask2 <- mask(clip2, v)           

# 关键：确保人口栅格与热日栅格分辨率/范围完全对齐（若原数据不一致则执行）
mask2 <- resample(mask2, mask1, method = "bilinear") # 可根据需求选择重采样方法

# 直接计算加权总和与总人口总和
sum_pop_hd <- global(mask1 * mask2, fun = "sum", na.rm = TRUE)[[1]]
sum_pop <- global(mask2, fun = "sum", na.rm = TRUE)[[1]]

# Population-weighted number of hot days - final estimate
pwhd <- sum_pop_hd / sum_pop

优化点说明 #

• 移除嵌套循环：terra的栅格乘法mask1 * mask2会自动对每个对应单元格执行运算，底层用C++实现，比R循环快数十倍甚至上百倍
• 消除冗余操作：原代码中创建grid_raster、循环内转多边形再mask的操作完全多余，因为mask2已经是裁剪后的栅格，直接使用即可
• 栅格对齐处理：通过resample确保两个栅格的单元格完全对应，避免因分辨率不一致导致的计算错误
• 批量求和：global函数一次性计算全局总和，替代循环内的逐个单元格求和

内容的提问来源于stack exchange，提问作者Mikaël