R语言操作栅格对象触发C栈使用超限错误的排查求助

问题描述 #

在Windows 10设备上使用R 4.2.3及RStudio 2023.03处理栅格数据时，遇到Error: C stack usage 15925776 is too close to the limit错误。该栅格对象来自同事发送的.RData文件，常规栅格操作（如r+1、mean(r)、summary(r)等）均触发错误，但object.size(r)、print(r)等仅访问元数据的操作可正常执行。

每次报错伴随50条W:\rasters\trail_road_density_formal.tif: No such file or directory (GDAL error 4)警告，推测该路径为同事设备上的文件位置。已尝试调整栈大小，但无效，需无需重新获取原始.tif文件的解决方案。

环境信息 #

R version 4.2.3 (2023-03-15 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=English_Canada.utf8  LC_CTYPE=English_Canada.utf8    LC_MONETARY=English_Canada.utf8
[4] LC_NUMERIC=C                    LC_TIME=English_Canada.utf8    

attached base packages:
[1] graphics  grDevices datasets  stats     utils     methods   base     

other attached packages:
 [1] terra_1.7-23    moveHMM_1.8     CircStats_0.2-6 boot_1.3-28.1   MASS_7.3-58.2   maptools_1.1-5 
 [7] rgdal_1.6-5     sf_1.0-8        raster_3.6-20   sp_1.5-0        lubridate_1.8.0 forcats_0.5.2  
[13] stringr_1.5.0   dplyr_1.1.0     purrr_0.3.5     readr_2.1.3     tidyr_1.2.1     tibble_3.1.8   
[19] ggplot2_3.4.1   tidyverse_1.3.2

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10         lattice_0.20-45     class_7.3-21        png_0.1-7           utf8_1.2.2         
 [6] R6_2.5.1            cellranger_1.1.0    plyr_1.8.7          backports_1.4.1     reprex_2.0.2       
[11] e1071_1.7-11        httr_1.4.5          pillar_1.8.1        RgoogleMaps_1.4.5.3 rlang_1.0.6        
[16] googlesheets4_1.0.1 readxl_1.4.1        rstudioapi_0.14     geosphere_1.5-14    googledrive_2.0.0  
[21] foreign_0.8-84      munsell_0.5.0       proxy_0.4-27        broom_1.0.3         compiler_4.2.3     
[26] numDeriv_2016.8-1.1 modelr_0.1.10       pkgconfig_2.0.3     tidyselect_1.2.0    codetools_0.2-19   
[31] fansi_1.0.3         crayon_1.5.2        tzdb_0.3.0          dbplyr_2.3.1        withr_2.5.0        
[36] bitops_1.0-7        grid_4.2.3          jsonlite_1.8.2      gtable_0.3.1        lifecycle_1.0.3    
[41] DBI_1.1.3           magrittr_2.0.3      units_0.8-0         scales_1.2.1        KernSmooth_2.23-20 
[46] cli_3.4.1           stringi_1.7.8       fs_1.5.2            xml2_1.3.3          ellipsis_0.3.2     
[51] generics_0.1.3      vctrs_0.5.2         tools_4.2.3         ggmap_3.0.1         glue_1.6.2         
[56] hms_1.1.2           jpeg_0.1-9          colorspace_2.0-3    gargle_1.3.0        classInt_0.4-8     
[61] rvest_1.0.3         haven_2.5.1   

可正常执行的操作示例 #

> object.size(r)
14584 bytes

> print(r)
class      : RasterLayer 
dimensions : 5627, 5008, 28180016  (nrow, ncol, ncell)
resolution : 30.00193, 29.99805  (x, y)
extent     : 498087.4, 648337.1, 5599179, 5767978  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs 
source     : trail_road_density_formal.tif 
names      : trail_road_density_formal 
values     : 0, 36.11524  (min, max)

> str(r)
Formal class 'RasterLayer' [package "raster"] with 13 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "W:\\rasters\\trail_road_density_formal.tif"
  .. .. ..@ datanotation: chr "FLT4S"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ bandorder   : chr "BIL"
  .. .. ..@ offset      : int 0
  .. .. ..@ toptobottom : logi TRUE
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 5008
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
  .. .. ..@ values    : logi(0) 
  .. .. ..@ offset    : num 0
  .. .. ..@ gain      : num 1
  .. .. ..@ inmemory  : logi FALSE
  .. .. ..@ fromdisk  : logi TRUE
  .. .. ..@ isfactor  : logi FALSE
  .. .. ..@ attributes: list()
  .. .. ..@ haveminmax: logi TRUE
  .. .. ..@ min       : num 0
  .. .. ..@ max       : num 36.1
  .. .. ..@ band      : int 1
  .. .. ..@ unit      : chr ""
  .. .. ..@ names     : chr "trail_road_density_formal"
  ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
  .. .. ..@ type      : chr(0) 
  .. .. ..@ values    : logi(0) 
  .. .. ..@ color     : logi(0) 
  .. .. ..@ names     : logi(0) 
  .. .. ..@ colortable: logi(0) 
  ..@ title   : chr(0) 
  ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
  .. .. ..@ xmin: num 498087
  .. .. ..@ xmax: num 648337
  .. .. ..@ ymin: num 5599179
  .. .. ..@ ymax: num 5767978
  ..@ rotated : logi FALSE
  ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
  .. .. ..@ geotrans: num(0) 
  .. .. ..@ transfun:function ()  
  ..@ ncols   : int 5008
  ..@ nrows   : int 5627
  ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs"
  .. .. ..$ comment: chr "PROJCRS[\"NAD83 / UTM zone 11N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n  "| __truncated__
  ..@ history : list()
  ..@ z       : list()
  ..@ NA      : NULL
Warning message:
Not a validObject(): no slot of name "srs" for this object of class "RasterLayer"

已尝试的无效操作 #

> Cstack_info()
      size    current  direction eval_depth 
  15922790      21688          1          2 
> system("ulimit -s")
[1] 127
> system("ulimit -s unlimited")
[1] 127
> r + 1
Error: C stack usage  15927824 is too close to the limit

解决方案 #

核心原因 #

栅格对象仅保存了同事本地tif文件的引用（@file@name指向W盘路径），实际数据未存储在.RData中。当执行需要读取数据的操作时，GDAL反复尝试访问不存在的文件，最终导致栈溢出错误。

方法1：强制将栅格数据加载至内存（若可行） #

若同事保存.RData时，栅格数据已被缓存或可通过某种方式读取，使用readAll()函数将数据全部加载至内存，脱离对外部文件的依赖：

library(raster)
r <- readAll(r)
# 验证是否加载成功
print(r@data@inmemory) # 应返回TRUE

加载完成后，常规栅格操作即可正常执行。

方法2：修改栅格文件引用并重建空栅格（避免栈溢出） #

若无法读取原始数据，可通过修改栅格的文件引用，创建一个匹配元数据的空栅格，避免GDAL反复报错导致栈溢出：

library(raster)
# 创建临时tif文件
temp_tif <- tempfile(fileext = ".tif")
# 基于现有栅格元数据生成空栅格并写入临时文件
writeRaster(
  raster(extent(r), nrow=nrow(r), ncol=ncol(r), crs=crs(r)),
  temp_tif,
  overwrite = TRUE
)
# 修改栅格对象的文件路径
r@file@name <- temp_tif

此方法可消除栈溢出错误，但栅格数据为空，仅保留元信息。若需要原始数据，仍需联系同事重新发送包含数据的.RData或tif文件。

方法3：转换为terra包的SpatRaster尝试修复 #

terra包对栅格的处理更稳定，可尝试将RasterLayer转换为SpatRaster，自动处理路径问题：

library(terra)
# 转换为SpatRaster
x <- rast(r)
# 将数据移入内存
x <- mem(x)
# 转换回RasterLayer（若仍需使用raster包）
r <- raster(x)

若转换过程中能成功恢复数据，后续操作即可正常执行。

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