R中rasterToPoints失效时，如何从.tif栅格提取像元值与坐标

问题背景 #

从WorldPop的澳大利亚2020年人口约束调整.tif文件中提取像元值及对应x、y坐标时，遇到以下问题：

• 使用raster包的rasterToPoints()函数时，R会话直接卡住，无法生成空间点数据框（spdf）
• 尝试getValues()或readAll()提取时，因栅格过大报错：Error: cannot allocate vector of size 17.8 Gb

环境信息 #

• R版本：4.2.0（Windows10 x64）
• 系统内存：总内存31.781 GiB，空闲26.164 GiB
• 注意：R 4.2.0中memory.limit()已弃用，无法通过该命令调整内存

尝试过的代码 #

library(raster)
# 读取栅格文件
Raster <- raster("C:/file path/aus_ppp_2020_UNadj_constrained.tif")
# 尝试转换为空间点数据框
Raster <- rasterToPoints(Raster, spatial = TRUE)

内存与系统信息代码及输出：

sessionInfo()
# R version 4.2.0 (2022-04-22 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 22000)

library(memuse)
memuse::Sys.meminfo()
# Totalram:  31.781 GiB 
# Freeram:   26.164 GiB 

memory.limit() 
# Warning: 'memory.limit()' is no longer supported[1] Inf

需求 #

寻求以下三类解决方案：

1. 让rasterToPoints()函数能正常处理该栅格
2. 栅格分块处理后合并结果的方法
3. 提取像元值与坐标的替代方案

解决方案 #

方案1：优化rasterToPoints()参数 #

默认rasterToPoints()会提取所有像元（包括NA值），大幅增加内存占用。通过fun参数过滤无效像元，只保留有效数据：

library(raster)
Raster <- raster("C:/file path/aus_ppp_2020_UNadj_constrained.tif")
# 仅提取人口值大于0的像元（根据数据实际情况调整过滤条件）
spdf <- rasterToPoints(Raster, spatial = TRUE, fun = function(x) x > 0)

若数据中NA占比高，该方法能显著减少处理的像元数量，避免内存溢出或会话卡住。

方案2：分块处理栅格并合并结果 #

将大栅格切割为多个小瓦片，逐个处理后合并最终结果：

library(raster)
library(sp)

Raster <- raster("C:/file path/aus_ppp_2020_UNadj_constrained.tif")
n_rows <- nrow(Raster)
n_cols <- ncol(Raster)

# 设置分块大小（根据内存情况调整，比如每块1000行）
block_size <- 1000
n_blocks <- ceiling(n_rows / block_size)

# 初始化空的空间点数据框
final_spdf <- SpatialPointsDataFrame()

for (i in 1:n_blocks) {
  start_row <- (i - 1) * block_size + 1
  end_row <- min(i * block_size, n_rows)
  
  # 读取当前块的栅格数据
  block <- raster(Raster, row=start_row, nrow=end_row)
  # 转换为点数据框（可添加fun过滤NA）
  block_spdf <- rasterToPoints(block, spatial = TRUE, fun = function(x) x > 0)
  
  # 合并结果
  if (i == 1) {
    final_spdf <- block_spdf
  } else {
    final_spdf <- rbind(final_spdf, block_spdf)
  }
  
  # 清理临时变量释放内存
  rm(block, block_spdf)
  gc()
}

注意：若仍报错，可调小block_size（如500）；每次循环后调用gc()强制垃圾回收，释放内存。

方案3：替代工具/包处理 #

方法A：使用terra包（raster升级版本）

terra内存效率更高，处理大栅格更友好，内部自动分块：

library(terra)
# 读取栅格
r <- rast("C:/file path/aus_ppp_2020_UNadj_constrained.tif")
# 转换为点数据框
points_df <- as.points(r, values=TRUE)
# 转为sp格式空间点数据框（若需要）
library(sp)
spdf <- as(points_df, "SpatialPointsDataFrame")

方法B：GDAL命令行工具提取（绕开R内存限制）

用GDAL直接将栅格转为CSV，再导入R：

1. 执行GDAL命令（需提前安装GDAL）：

gdal_translate -of XYZ "C:/file path/aus_ppp_2020_UNadj_constrained.tif" "output_points.csv"

2. R中读取CSV并转换为空间点数据框：

points_df <- read.csv("output_points.csv", header=FALSE, col.names=c("x", "y", "population"))
library(sp)
coordinates(points_df) <- ~x+y
# 需匹配原栅格的坐标系，示例为WGS84
proj4string(points_df) <- CRS("+proj=longlat +datum=WGS84")

内容的提问来源于stack exchange，提问作者Mikin