★ FEATURED ARTICLE 1. 需求场景与项目思路1.1 为什么要批量提取栅格统计值做遥感或空间数据分析的朋友应该都遇到过这种场景拿到一个文件夹里面铺了几十上百幅栅格影像可能是不同年份的NDVI、不同波段的反射率也可能是模型输出的连续网格数据。老板或合作方开口就问这些数据的整体范围是多少波动大不大平均水平和离散程度怎么样这时候如果你还在一幅图一幅图地打开软件去看属性、记数字效率低不说还容易抄错。更关键的是当文件夹里有几百个文件时手点已经完全不现实了。这个需求本质上是批量栅格描述性统计要提取的核心指标就是四个最小值min、最大值max、标准差sd和平均值mean。这四个指标加起来基本能刻画一个栅格数据集的分布特征min和max给出值域范围mean代表整体水平sd反映空间异质性。理解每个栅格的这四个值是快速质控数据、筛选异常图层、对比不同时期影像变化的第一步。R语言在这个场景里非常合适。一方面raster和terra这套栅格处理栈成熟稳定读取、统计、批量循环都方便另一方面R的data.frame和write.csv让结果整理和导出特别顺畅。我建议用terra包作为主力它是raster的下一代版本内存管理和速度都有明显提升而且函数命名更统一。对于从没接触过terra的朋友这篇里的代码同样能帮你快速上手。1.2 项目目标与适用对象这个项目的目标很简单写一段R脚本输入一个文件夹路径自动读取里面所有栅格文件依次计算每个栅格的min、max、sd和mean把结果汇总成一张表输出。既能一次性跑完整个文件夹也能按特定后缀筛选文件还考虑到NoData值、CRS不一致、堆栈文件等实际数据里的坑。适合谁来参考首先是刚接触R栅格处理的初学者可以照着代码逐行跑通理解批量处理的套路其次是经常处理多幅影像的从业者可以把脚本存成模板改改路径就能复用到不同项目里最后是那些文件量特别大、之前靠手工统计的朋友这篇里的方法能帮你把几分钟的重复劳动压缩成几秒钟。代码不复杂但背后涉及的批量操作逻辑、内存管理思路和栅格统计原理值得好好捋一遍。2. 环境准备与核心工具选型2.1 R语言与必要的包开始写代码之前先把环境准备好。R的安装就不展开了注意两点一是建议用最新版R4.x二是RStudio顺手装上写脚本、看变量、调试都方便。针对栅格处理需要安装的关键包是terra它同时承担了数据读取、统计计算和输出三大任务。install.packages(terra)如果是在公司内网或者镜像源受限的环境可以指定镜像安装install.packages(terra, repos https://mirrors.tuna.tsinghua.edu.cn/CRAN/)除了terra强烈建议装上dplyr或data.table来整理输出结果。不过我们的统计结果量级很小文件数行 × 指标数列用R自带的data.frame就完全够了没必要额外引入依赖保持脚本轻量更好。如果以后要扩展到上千个文件、还要合并更多属性再考虑data.table也不迟。加载包和设置工作目录的固定动作library(terra) # 设置工作目录建议用绝对路径避免相对路径带来的混乱 setwd(E:/R_projects/raster_stats)2.2 为什么选择terra而不是raster很多老教程还在用raster包但2023年起raster包基本进入维护模式新功能都集中在terra上。我自己是两个包都用过明显感觉到terra在三个地方更顺手处理速度快。terra基于C的底层架构重写了栅格IO和计算逻辑对大文件的读取和聚合操作比raster快一个量级。我有一个900MB的DEM文件raster读一遍要十几秒terra压缩到三秒上下。做批量统计时文件一多这个差距就很致命了。内存管理更智能。栅格数据动辄几百MB甚至几个GBraster处理大文件时会频繁地整块读入内存卡顿甚至崩掉是常事。terra默认的分块读取机制更激进能把内存占用控制在一个合理范围。当然这跟数据本身也有关系下面会说怎么进一步优化。函数设计更统一。raster包里有cellStats、calc、extract一堆函数各有各的用法terra里直接就是global()、app()、extract()参数和返回值更规范。批量统计这种活用terra写出来的代码明显更短。2.3 栅格文件夹结构说明为了后面演示方便我先说明一下数据组织方式。假设我们的文件夹长这样E:/R_projects/raster_stats/data/ ├── ndvi_2024_01.tif ├── ndvi_2024_02.tif ├── ndvi_2024_03.tif └── ndvi_2024_04.tif这些是月度NDVI栅格覆盖同一研究区域分辨率、范围理论上一致但实际情况中也可能有偏差。批量统计的脚本要考虑这种不一致性不能假设所有文件完全对齐。如果你的数据不是tif格式比如是.img、.ascterra同样支持只需在匹配文件后缀时做调整。3. 分步实现从单个栅格到批量统计3.1 第一步读取单个栅格并验证项目一开始先别急着写循环先把单幅栅格跑通。这样能把读取数据—计算指标—输出结果每个环节的问题单独暴露再进入批量化时心理就有底了。# 读取单个栅格 r - rast(data/ndvi_2024_01.tif) print(r)rast()是terra里读取栅格的核心函数不管是单波段还是多波段的GeoTIFF都能读。print出来的信息很关键能看到分辨率、范围、CRS、文件路径、波段数。我第一次接触terra时习惯忽略这步直接算后来发现很多报错都源于数据和预期的偏差好好看一眼print输出的信息能省一大半排查时间。计算四个核心指标在terra里一行搞定# 计算统计值 stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) print(stats)global()函数的逻辑是对整个栅格的所有像元应用指定函数返回一个标量或向量。这里的na.rm TRUE必须加上不加的话碰到NoData值结果就直接变成NA后面批量统计会全线崩盘。输出长这样min max sd mean lyr1 -0.213 0.876 0.231045 0.452887注意sd和mean的列名terra生成的默认列名就是函数名后面批量拼接时直接用。如果只想算单值比如只看平均值也可以写成global(r, mean, na.rm TRUE)不过标题要求四个指标一起返回所以用向量传给fun参数更简洁。3.2 第二步实现文件夹批量读取单幅验证通过后进入批量环节。第一步要解决的是怎么把文件夹里的tif文件全部找出来。基础做法是用list.files()# 获取所有tif文件 file_list - list.files(data, pattern \\.tif$, full.names TRUE) print(file_list)逐行拆解一下第一个参数是文件夹路径相对路径的写法是相对于当前工作目录。如果你设置了setwd直接写文件夹名就行。pattern \\.tif$是一个正则表达式\\.匹配点号$表示结尾。这样能精确匹配所有以.tif结尾的文件不会误伤.tiff或其他格式。如果数据是.img改成\\.img$即可。full.names TRUE是重中之重。不设这个参数返回的只是文件名比如ndvi_2024_01.tif后面传给rast()时会因为找不到文件而报错。设成TRUE后返回完整路径data/ndvi_2024_01.tif才可以直接用于读取。如果你还想对文件排序别依赖list.files()默认的返回顺序这个顺序在不同操作系统上不一样。稳妥做法是加一句sort()file_list - sort(list.files(data, pattern \\.tif$, full.names TRUE))这一步很重要因为后面结果表会和文件列表一一对应顺序错乱会导致文件A的统计结果写在文件B那一行这种隐蔽错误很难发现。3.3 第三步写循环逐个计算并汇总有了文件列表接下来就是经典的初始化结果对象 → 循环处理 → 汇总绑定三步走。先初始化一个空的结果表用data.frame预先声明列名# 初始化结果数据框 results - data.frame( file character(), min numeric(), max numeric(), sd numeric(), mean numeric(), stringsAsFactors FALSE )预声明列名有两个好处一是保证列类型正确character列不会被误判成factor二是后面每行绑定进来时结构始终一致不会出现第一次绑定后某一列变成list的情况。写循环# 循环处理每个文件 for (f in file_list) { # 读取栅格 r - rast(f) # 计算统计值 stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) # 提取文件名不含路径和后缀作为标识 file_short - basename(f) # 获取含后缀文件名 file_core - tools::file_path_sans_ext(file_short) # 去掉后缀 # 组装当前文件的结果行 row - data.frame( file file_core, min stats$min, max stats$max, sd stats$sd, mean stats$mean, stringsAsFactors FALSE ) # 绑定到总结果表 results - rbind(results, row) }代码逻辑非常直观每次循环读入一个栅格算出四个统计量组装成一行rbind进结果表。这里有个小技巧tools::file_path_sans_ext()可以干净地去掉后缀比用gsub去匹配更不容易出错尤其是文件名里可能有多余的点号时。循环结束后results里已经有所有文件的统计结果直接打印看看print(results)3.4 第四步结果排序与导出统计结果默认按文件列表的顺序排列如果文件列表没排过序结果就很随机。按文件名排序用order()results - results[order(results$file), ]这步是针对时间序列数据的常用操作。比如文件名叫ndvi_2024_01.tif、ndvi_2024_02.tif这种字符串排序和自然时间排序一致。但如果跨越年份比如有ndvi_2023_12.tif和ndvi_2024_01.tif字符串排序结果仍然一致因为年份在月份前面暂时不需要额外处理。如果将来要在循环里直接保证顺序也可以提前对file_list做自定义排序函数。导出到CSV方便后续Excel或Python继续处理write.csv(results, output/raster_stats_results.csv, row.names FALSE)建议先创建output文件夹再写文件不然首次运行时会报无法打开文件之类的错误。一个更稳的写法if (!dir.exists(output)) dir.create(output) write.csv(results, output/raster_stats_results.csv, row.names FALSE)到这一步核心功能已经全部实现。整个脚本去掉注释不到20行非常轻量却能解决一个实际工作中的高频痛点。4. 进阶优化处理真实数据中的各种意外4.1 文件后缀不统一怎么办实际情况里文件夹中的栅格可能既有.tif又有.img或者文件名混乱不堪。这时可以改用更宽泛的匹配方式把所有R能读取的栅格格式都纳入# 匹配常见栅格后缀 file_list - list.files( data, pattern \\.(tif|tiff|img|asc|dat)$, full.names TRUE, ignore.case TRUE )ignore.case TRUE是为了防止有人存了.TIF这种大写后缀。正则里的|是或的逻辑\\.(tif|tiff|img|asc|dat)$表示匹配所有以这些后缀结尾的文件。用这个方式后只要terra能读的格式都能批量处理。但也要小心万一文件夹里有辅助文件比如.tif.aux.xml这类用上面的正则会漏掉因为结尾不是.tif。如果出现读取失败或文件数量对不上的情况先用print(file_list)检查匹配结果再做针对性调整。4.2 NoData值与栅格统计准确性这是个很多人踩过的坑。栅格数据里经常用-9999、NA或0来表示无数据区域如果不去掉这些值统计结果就会严重失真。terra读取GeoTIFF时能自动识别文件里标注的NoData值但如果原始文件的NoData定义不规范或者你想手动指定可以在rast()读取后单独设置r - rast(f) # 如果NoData值没有被正确识别手动指定 NAflag(r) - -9999重要的是global()里的na.rm TRUE必须配合正确的NAflag才能生效。如果文件的NoData值是-9999但没有被识别成NA那么na.rm TRUE处理的是真正的NA-9999仍然会参与计算把均值拉低一大截。这种错误非常隐蔽因为你看到的平均值是有值的、看起来很合理的数字实际上已经被污染了。批量处理时怎么避免两种思路一种是循环里检查每个文件的NAflag发现异常就修正for (f in file_list) { r - rast(f) if (any(is.na(NAflag(r))) || (is.na(NAflag(r)[1]))) { # 某些没有定义NAflag的栅格默认NAflag为NA # 这时按实际数据的NoData值手动设置 } }另一种更稳妥统计前先做一次去极值处理把小于某个合理阈值的值直接设置成NA。比如NDVI的合理范围是-1到1之间超过这个范围的基本都是NoData或异常值r[r -1 | r 1] - NA这种方式适合已知数据合理取值范围的场景不只解决NoData问题还能剔除传感器异常值。缺点是需要你对数据的物理意义有足够了解阈值设错反而会引入更大误差。4.3 多波段栅格如何处理前面演示的都是单波段栅格。如果文件夹里的tif是四波段的卫星影像global()返回的是每个波段各一行统计结果。这时再按一个文件一行去组装数据就会出问题。处理多波段栅格的逻辑很简单拿到每个波段的统计结果后按行展开r - rast(f) # 假设有4个波段 stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) print(stats)输出会是min max sd mean lyr1 -0.213 0.876 0.231045 0.452887 lyr2 -0.101 0.654 0.187655 0.312876 lyr3 0.002 0.789 0.198654 0.423765 lyr4 0.111 0.923 0.201234 0.534567这时如果想要每个波段一行的结果表可以改成两层循环外层是文件内层是波段。为了标识清楚文件名后面加上波段编号for (f in file_list) { r - rast(f) stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) # 为每个波段生成一行 for (i in 1:nrow(stats)) { row - data.frame( file paste0(tools::file_path_sans_ext(basename(f)), _band, i), min stats$min[i], max stats$max[i], sd stats$sd[i], mean stats$mean[i] ) results - rbind(results, row) } }如果波段很多比如高光谱数据有几十上百个波段这样逐行rbind效率会变低更高效的做法是先把每个文件的波段统计结果存成一个list循环结束后一次性do.call(rbind, list)归并。4.4 大规模数据时的内存优化文件特别多或栅格特别大时循环里的每个rast(f)都会占用内存处理完的文件要及时释放。R的垃圾回收虽然自动运行但大对象占用的内存未必能及时还给系统。手动加一句gc()在循环末尾是常用的做法for (f in file_list) { r - rast(f) stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) # ... 组装结果 rm(r) # 删除大对象 gc() # 请求垃圾回收 }但说实话对大多数场景几百个tif、普通分辨率terra已经处理得很好不需要特别担心。真正需要注意的反倒是不要一次性把所有栅格都读入内存再批量统计。有些新手会写# 反面教材一次性读入所有文件再统计 all_rasters - lapply(file_list, rast) # 内存爆炸的根源这种做法在文件少时可以跑但文件一多内存直接见底甚至系统卡死。正确的打开方式是边读边统计、统计完就释放也就是前面的循环写法。terraform的分块机制已经帮我们做了底层优化我们要做的是不在业务逻辑层制造额外的内存压力。5. 常见问题与排查技巧实录5.1 读取文件时报错file does not exist这个报错八成是路径问题。list.files()返回的路径和当前工作目录不一致或者full.names FALSE导致读取时找不到文件。排查思路# 检查当前工作目录 getwd() # 检查文件列表是否包含完整路径 print(file_list) # 如果full.names FALSE手动拼接路径 file_list - paste0(data/, file_list)还有一种情况是Windows系统下路径分隔符的问题。R在Windows里接受/和\\两种分隔符但如果你从别的地方复制了C:\Users\name\data这种路径反斜杠会被解释成转义字符导致路径错误。建议统一用正斜杠/。5.2 统计结果全部为NA如果global()返回的统计值全是NA大概率是NoData值问题。分两步排查第一步打印栅格对象看NAflagr - rast(f) print(NAflag(r))如果NAflag显示为NA说明terra可以正常处理文件内定义的NoData。如果显示不是NA但明显不对比如显示为0但数据的真实NoData是-9999就需要手动设置NAflag(r) - -9999第二步直接看数据的值域分布# 检查像元值的独特性 unique_values - unique(values(r)) print(head(unique_values, 20))这样能看到数据里是否有-9999这类异常值判断是不是NoData没被识别。解决后重新统计。5.3 CRS坐标系不一致导致的问题统计min、max、sd、mean只关心像元值本身不涉及投影变换所以CRS不一致通常不会导致这个脚本报错。但在两种情况下要额外留意一是栅格范围extent差异很大有的文件覆盖全域有的只覆盖一个小块这时直接对比mean值没有意义需要加上范围大小作为辅助变量一起输出。方法是读取时记录像元数量或面积ncell_stats - ncell(r) # 加进results里作为一列二是后续如果要对栅格做叠加分析CRS不一致会直接报错。提前统一CRS可以这样# 将所有栅格重投影到同一CRS r_reprojected - project(r, EPSG:4326)但重投影会插值改变原始像元值尤其是分类数据所以提取统计值这个场景里不建议做投影变换除非你有明确的对比需求。5.4 速度慢到无法接受批量处理上百个几十MB的栅格循环可能要跑几分钟这是正常现象。但如果慢得离谱比如一个文件要十几秒建议检查一是global()是否在做全像元扫描。理论上统计min、max、sd、mean必须扫描所有像元没有捷径。但如果你用terra的global()底层已经是C实现比我早期用raster::cellStats()快很多。二是数据是否特别大。如果栅格分辨率是0.5米、覆盖整个县单文件就有几个GB读取需要时间。这时可以先用terra::aggregate()做降采样预览但会牺牲精度不推荐用于正式统计。三是后台有没有其他程序占用内存。R在处理大文件时内存不足会触发swap速度直线下降。关掉不必要的程序或者换一台大内存机器跑批量任务。5.5 结果写进CSV后再读入出现乱码或格式变化write.csv()写出的文件默认使用系统编码在Windows上通常是GBK在macOS或Linux上是UTF-8。用Excel打开GBK文件一般没问题但如果用Python或其他工具读可能遇到乱码。推荐显式指定编码write.csv(results, output/raster_stats_results.csv, row.names FALSE, fileEncoding UTF-8)还有一种情况是统计值本来保留了十几位小数CSV里看起来很长Excel打开后可能被截断显示。如果想要更规整的结果在写CSV之前先roundresults$min - round(results$min, 4) results$max - round(results$max, 4) results$sd - round(results$sd, 4) results$mean - round(results$mean, 4)保留四位小数对大多数应用场景足够也可以根据自己的精度需求调整。6. 把脚本固化成可复用函数6.1 封装成通用函数做了几次类似的项目后我建议把这段逻辑封装成函数以后一行代码调用省去每次重写循环的功夫。函数里把文件格式、NoData值、输出路径都设成参数灵活调整batch_raster_stats - function(folder_path, pattern \\.tif$, output_file NULL, na_value NULL, round_digits 4) { library(terra) # 获取文件列表 file_list - sort(list.files(folder_path, pattern pattern, full.names TRUE)) if (length(file_list) 0) { stop(No files found in the specified folder.) } # 初始化结果 results - data.frame( file character(), min numeric(), max numeric(), sd numeric(), mean numeric(), stringsAsFactors FALSE ) # 循环处理 for (f in file_list) { r - rast(f) # 设置NoData值 if (!is.null(na_value)) { NAflag(r) - na_value } # 计算统计值 stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) # 组装一行 file_core - tools::file_path_sans_ext(basename(f)) row - data.frame( file file_core, min stats$min, max stats$max, sd stats$sd, mean stats$mean, stringsAsFactors FALSE ) results - rbind(results, row) } # 排序并四舍五入 results - results[order(results$file), ] if (!is.null(round_digits)) { results[, 2:5] - round(results[, 2:5], round_digits) } # 导出 if (!is.null(output_file)) { if (!dir.exists(dirname(output_file))) { dir.create(dirname(output_file), recursive TRUE) } write.csv(results, output_file, row.names FALSE, fileEncoding UTF-8) } return(results) }调用只需要一行res - batch_raster_stats( folder_path data, pattern \\.tif$, output_file output/stats.csv, na_value -9999, round_digits 4 )如果以后要处理.img文件直接改pattern参数即可。如果某个文件夹的NoData值是-9999但另一个不是每次调用时传入对应值就行不用改函数本身。6.2 扩展同时处理多个文件夹有些项目的输入数据分布在多个子文件夹里比如按年份分了data/2023/、data/2024/。可以用list.files(recursive TRUE)一键递归拉取全部子文件夹里的文件file_list - list.files( data, pattern \\.tif$, full.names TRUE, recursive TRUE )这样找出来的文件路径里包含子文件夹名但basename()只保留文件名跨年份同名文件比如两年的数据都叫ndvi_01.tif会被后面的排序和展示搞混。解决方法是在文件名前列里带上相对路径# 把file改成包含相对路径的标识 file_core - sub(^data/, , tools::file_path_sans_ext(f))或者在结果表中增加一列folder从路径里提取子文件夹名folder_name - basename(dirname(f)) row - data.frame( folder folder_name, file file_core, ... )这样导出后可以直接用数据透视表pivot table按文件夹分组看统计结果非常方便。7. 实际项目演练以月度NDVI数据为例为了让你看到这套方法的完整应用效果我用一个模拟场景走一遍全流程。假设data/文件夹里有2024年1月到4月的月度NDVI栅格每幅图覆盖同样区域范围是0到1之间NDVI的取值范围空间分辨率250米。目标是计算每月的min、max、sd、mean并输出一张汇总表。完整脚本library(terra) # 数据路径 data_dir - data output_dir - output if (!dir.exists(output_dir)) dir.create(output_dir) # 文件列表 file_list - sort(list.files(data_dir, pattern \\.tif$, full.names TRUE)) # 初始化结果表 results - data.frame( month character(), min numeric(), max numeric(), sd numeric(), mean numeric(), stringsAsFactors FALSE ) # 批量计算 for (f in file_list) { # 读取栅格 r - rast(f) # 手动指定NoData NAflag(r) - -9999 # 空间聚合前先做必要的裁剪如果范围不统一按公共范围裁剪 # 这一步可选需要以第一幅图为基准 # r - crop(r, rast(file_list[1])) # 计算统计值 stats - global(r, fun c(min, max, sd, mean), na.rm TRUE) # 提取文件名中的月份假设为 ndvi_2024_01.tif 格式 month_label - sub(ndvi_(.*)\\.tif, \\1, basename(f)) # 组装行 row - data.frame( month month_label, min round(stats$min, 4), max round(stats$max, 4), sd round(stats$sd, 4), mean round(stats$mean, 4), stringsAsFactors FALSE ) # 汇总 results - rbind(results, row) } # 按月份排序 results - results[order(results$month), ] # 输出 print(results) write.csv(results, output/ndvi_monthly_stats.csv, row.names FALSE, fileEncoding UTF-8)跑完后你会得到一张类似这样的表monthminmaxsdmean2024_010.0210.8760.2310.4532024_020.0140.8990.2450.4872024_030.0450.9340.2180.5422024_040.0670.9580.1970.601从这张表能读出什么mean从1月的0.453上升到4月的0.601说明植被覆盖在生长季逐步增加sd整体稳定在0.2左右说明空间异质性没有发生剧烈变化max在0.9上下浮动说明始终存在高植被覆盖区域。这些判断如果只看单个栅格很难形成整体印象但汇总成一张表后趋势一目了然。8. 经验技巧与日常使用心得8.1 先小范围测试再全量运行写批量脚本最忌讳的就是直接在几百个文件上跑。我习惯的做法是先拿3到5个文件测试确认输出结果合理后再替换成完整文件列表。测试文件的选取要有代表性最好包含一个边界案例——比如范围最小的、NoData最多的、文件名带特殊字符的。这样可以一次性暴露大部分潜在问题。测试时的技巧是把file_list切片# 只取前5个文件做测试 test_list - file_list[1:5]跑通后再把test_list换回file_list。如果测试文件的结果看起来合理全量运行基本不会出大问题。8.2 记录运行日志如果数据量大、运行时间长建议在脚本里打日志记录每个文件的处理状态。这样即使中途报错也能快速定位到具体是哪个文件出了问题for (i in seq_along(file_list)) { f - file_list[i] message(Processing , i, /, length(file_list), : , basename(f)) # ... 原有逻辑 }message()会把信息输出到R的控制台又不影响结果对象的生成。如果日志要多还可以配合sink()把输出保存到txt文件里跑完再慢慢查。8.3 边学边用的拓展思路掌握了基础的批量栅格统计后这个项目还可以向多个方向延展按区域掩膜统计。如果只关心某个行政区或研究区的统计值可以先读入一个面状矢量文件shapefile或GeoJSON用terra::mask()或crop()裁出研究区域再统计。代码只增加一两行但结果的专业价值高很多尤其是做区域对比分析时。时间序列变化检测。如果文件是逐时期的影像算完每个时期的mean后可以继续做趋势分析。比如用lm()拟合时间序列的斜率快速识别植被退化或城市扩张区域。统计指标从单幅图的描述统计升级成跨时相的格局分析。与Python的衔接。R算完统计结果并导出CSV后可以无缝交给Python的pandas、matplotlib做可视化和后续分析。我在实际项目里经常是R做栅格预处理和统计Python画图出表各取所长。回到标题本身——提取文件夹内栅格数据的最小值、最大值、标准差和平均值这看起来是个很小的需求但背后涉及的路径处理、批量循环、NoData管理、结果汇总是R处理空间数据的通用基本功。把这套逻辑吃透了以后遇到类似的批量任务比如批量裁剪、批量重投影、批量计算面积都能举一反三。我自己第一次完整跑通这段代码后处理栅格数据的效率明显上了一个台阶希望这篇分享也能帮上你的忙。
阅读完成 · 觉得有帮助？
