使用大型地理空间数据源计算人口加权的炎热天数

Question

我正在计算一组国家多年来的人口加权炎热天数。我必须根据网格单元中的人口对网格单元中的每个炎热日值进行加权，以计算人口加权炎热日数的最终估计值。我有一个有效的代码，但与大国打交道确实需要很长时间。这是由嵌套循环引起的，我在其中有效地单独循环遍历每个网格单元。因此，该代码的时间效率不高，我正在寻找任何建议来更有效地编码。

我无法添加任何数据，因为这些数据涉及大型地理空间文件。

# 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)           

# Prepare grid raster to loop through
grid_rows <- nrow(mask1)  # Number of rows in the grid 
grid_cols <- ncol(mask1)  # Number of columns in the grid        

# Specify cell size 
cell_size <- res(mask1)[[1]]        

# Calculate the extent of the grid based on the raster's extent 
grid_extent <- ext(mask1)        

grid_raster <- rast(nrows = grid_rows, ncols = grid_cols, nlyrs = 1, crs = 'ESRI:54009', extent = grid_extent, resolution = cell_size)        

# Base value for calculating population-weighted number of hot days 
sum_pop_hd <- 0
sum_pop <- 0        

# Loop through each cell       
for (row in 1:grid_rows) {
        for (col in 1:grid_cols) {
          hd_cell <- mask1[row, col][[1]]                      
          if (is.na(hd_cell) == FALSE) {                         
          # Extract the cell as a SpatRaster             
          cell <- grid_raster[row, col, drop = FALSE]             
          cell <- terra::as.polygons(cell)                          

          # Produce values per cell             
          pop_cell <- terra::mask(mask2, cell, touches = FALSE)             
          pop_cell_sum <- terra::global(pop_cell, fun = "sum", na.rm = T)             
          pop_cell_sum <- pop_cell_sum[[1]]                         

            if (is.na(pop_cell_sum) == TRUE) {               
            sum_pop <- sum_pop # basically not change in values               
            sum_pop_hd <- sum_pop_hd             
            } else if (is.na(pop_cell_sum) == FALSE | ((pop_cell_sum == 0) == TRUE)) {               
            sum_pop <- sum_pop + pop_cell_sum               
            weighted <- pop_cell_sum * hd_cell               
            sum_pop_hd <- sum_pop_hd + weighted            
            }                         
          }                            
}              
}        

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

Answer 1

没有数据，很难猜测你的输出。然而，所有这些循环似乎都是多余的并且没有必要（对于光栅操作，循环通常是不必要的）。 我将使用一些虚拟数据和示例代码来说明我的想法。 首先，我将跳过 shapefile 边界上的裁剪和遮罩部分，因为此操作可以为每个国家/地区 shp 执行一次，并且不会更改主要工作流程。 其次，我假设人口和炎热天气栅格位于相同的 crs 中，具有相同的范围和分辨率。 第三，我会将一个栅格网格转换为矢量多边形，并使用该网格来提取炎热天气的加权平均值。

library(terra)
library(exactextractr)
library(sf)
pop <- rast(nrows=10, ncols=10)
values(pop) <- round(runif(100),2)*100
pop_vect <- as.polygons(pop,dissolve=F)
hd <- rast(nrows=10, ncols=10)
values(hd) <- round(runif(100),2)*10

wm <- exact_extract(hd,sf::st_as_sf(pop_vect),fun="weighted_mean",weights=pop)

这只是一个示例，但我认为您可以对其进行调整以获得您需要的内容。如果您可以提供有关栅格的更多信息，那么将很容易提供帮助。例如，

hot_days

和

pop

的分辨率、层数、范围、crs 等