最近我整理一批长时间序列的栅格数据时遇到了一个特别常见的需求同一个文件夹里放了几十个GeoTIFF文件每个文件代表不同时间或不同区域的地表参数而我需要把每个文件的最小值、最大值、标准差和平均值批量提取出来再做后续的整理和分析。说实话这个需求看起来简单真操作起来却有不少容易踩的坑比如文件里带NoData值、坐标系和范围不一致、文件名顺序混乱、文件太大导致内存不足等等。我之前也用过ArcGIS里的栅格统计功能数据量少的时候还能应付文件一多就麻烦了一个个右键查看图层属性再把结果抄到表格里不仅效率低还容易抄错。后来我直接改用R写了一个批量提取脚本只需要把文件夹路径改一下几分钟就能把整个目录跑完输出一张包含所有统计值的表格。这篇文章就把我这一套思路和代码完整记录下来主要围绕R语言和栅格数据的处理展开适合刚入门R的人也适合已经有一定基础但想提高数据处理效率的从业者。1. 这个任务到底在解决什么问题在开始写代码之前我非常建议先想清楚一个事我们要提取的这几个统计值实际应用里到底用来做什么为什么要写脚本而不是直接在GIS软件里看1.1 栅格统计值的使用场景栅格数据的“最小值、最大值、平均值、标准差”听起来像是很基础的描述性统计但不同领域都有非常实际的用途。比如做遥感反演的时候我经常要比较不同时期NDVI栅格的整体变化。直接看平均值能快速判断这一段时期植被覆盖是上升还是下降看标准差能反映空间差异是否变大如果标准差越来越高意味着研究区内土地覆盖状况在分化。再看最小值和最大值可以快速判断有没有异常像元比如小于0或大于1的值如果出现在NDVI里多半是数据预处理环节出了问题需要回头查掩膜或云检测结果。再比如做气象插值研究拿到一批逐月的降水栅格经常要先批量统计各月份的空间平均值和极值用来判断气候异常或者为后续的时空建模提供输入。手动在ArcGIS里拉一个月的图层属性容易但连续做120个月还要记录成表格那就不是工作量的问题而是人很容易出错。除了这些在精度验证、模型校准、图层批量检查这些环节里也都会用到四个统计指标。本质上这是一个“批量数据体检”的过程先把每个文件的内部情况摸清楚才能放心地进入下一步分析。1.2 为什么用R来做这件事用R处理这个需求优势在于可复制、可批量、结果可进一步分析。如果在QGIS或ArcGIS里手动操作流程大致是逐个打开栅格、查看元数据、打开统计信息、复制结果、粘到Excel。第一次做可能觉得还好但只要你需要换一批数据或者换一个文件夹就得重新操作一遍。而R脚本不一样路径一换回车一敲自动跑完整个过程还能给别人复现。更关键的一点是R处理完统计之后结果可以直接进数据框接着和站点数据、区域统计数据、时间序列信息合并甚至马上画图。我以前经常遇到一个场景跑完统计值之后还要把每个文件对应的年份、季节、实验组别加进去再做一个趋势分析。如果统计结果是在Excel里做的下一步又得重新读入R中间多了一步文件来回转的流程还可能出现格式不对齐的问题。所以我的结论是如果只是偶尔看一个两个栅格的统计值直接在GIS软件里点几下也没问题但凡你手头有超过5个栅格文件需要做同样的统计R是更高效、更稳妥的选择。2. 栅格处理工具选型raster包还是terra包在R里处理栅格数据早期最流行的是raster包这几年越来越多的人转向terra包。我现在的建议很明确新项目直接学terra除非你要维护老代码否则不要从raster开始。2.1 raster包和terra包的对比raster包曾经是R语言栅格分析的事实标准网上大量教程都在用raster包和它的rgdal依赖链。但问题在于raster包底层依赖的rgdal、sp等包已经停止维护官方也已经宣布进入“维护模式”在R的更新版本里安装越来越麻烦。terra包是同一个作者基于C重写的新一代栅格分析工具包你可以理解成raster的“完全重做版”。它继承了raster的核心设计思路但在速度上有了很大提升尤其是处理大型栅格的时候terra底层用S4类系统支持更多的栅格格式和sf矢量框架的配合也更好。我做过一个小对比同样是读取一个约500MB的GeoTIFF并计算全图统计值terra的耗时大概是raster的三分之一到一半。虽然单次操作差别看起来不大但在批量处理几十个文件时节省的时间就很可观了。2.2 为什么不建议用“自带”的summary或cellStats早期有些R用户遇到这个需求会想到直接对栅格对象执行summary()或者使用raster包里的cellStats()函数。这两个方法不是不行但都有一个共同的问题输出信息比较模糊而且批量处理时需要额外写循环代码来整理结果。terra包里的global()函数是我最常用的它可以直接指定一个或多个统计函数返回一个整齐的数据框。你用这个函数处理单张栅格感觉只是比summary()多了一步但放在批量循环里数据结构的一致性优势就特别明显。可能有人会问“那我直接读取栅格到内存然后转成向量再用R自带的min()、max()、mean()、sd()不就行了吗”理论上也可以但实际执行时一个几百MB的栅格转成向量后内存占用可能飙升到几个GB很容易造成R会话卡死或崩溃。terra的global()函数在底层是按块读取和计算的不会一次把全部像元载入内存这才是处理大数据时更安全的方式。3. 数据准备与核心代码实现接下来是整篇博文最核心的部分完整实操代码。我会讲得比较细每一步都说说为什么这么写方便你直接拿去改路径复用。3.1 安装并加载必要R包首先要确保R环境已经装好。如果还没安装terra用下面的命令install.packages(terra)安装完成后加载library(terra)如果你使用的是RStudio可以顺手在右上角检查R版本是否在4.1以上。terra包目前对R版本有一定要求太老的版本会出现安装失败或函数不兼容的情况。如果网络环境不太顺畅安装时还可以指定镜像源比如使用清华或中科大镜像install.packages(terra, repos https://mirrors.tuna.tsinghua.edu.cn/CRAN/)这里我额外提醒一句如果你同时装过raster包加载terra之后在同一个R会话里尽量不要混用两类栅格对象。因为terra的SpatRaster和raster的RasterLayer不是同一个类函数之间经常无法直接互相转换一旦混用很容易出现“找不到方法”或者“invalid class”的报错。3.2 设置文件夹路径并列出所有栅格文件首先把路径赋值给一个变量这样后面管理起来方便。folder_path - D:/R_work/raster_data接下来用list.files()列出文件夹内所有tif文件。我建议在pattern参数里明确写明要匹配的文件后缀files - list.files(folder_path, pattern \\.tif$, full.names TRUE) files这里的\\.tif$是正则表达式意思是匹配以.tif结尾的文件名。之所以这样写是因为文件夹里可能还有其他非栅格文件比如.tfw世界文件、.csv表格、.txt说明文档如果不加过滤后面读取时很可能报错。full.names TRUE也特别重要。如果不写这个参数list.files()返回的只是文件名不带路径后面rast()读取时R会在当前工作目录里找文件容易找不到。把完整路径返回出来才能确保程序在任何工作目录下都能正确定位文件。如果你还希望递归读取子文件夹里的所有tif文件可以把参数改成files - list.files(folder_path, pattern \\.tif$, full.names TRUE, recursive TRUE)不过我要提醒一下开启recursive TRUE后如果某个子目录里存放了临时文件或者不需要分析的备份tif它们也会被一起读进来。所以我更推荐的做法是先用上面的代码列出文件然后人工扫一眼files向量确认数量对不对、路径对不对再继续往下走。3.3 读取第一个栅格并查看基本信息批量循环之前我习惯先读取第一个文件确认数据能不能正常打开、坐标系和分辨率是否符合预期。r1 - rast(files[1]) r1打印出来的信息会包含像元数量、分辨率、范围、坐标系等。比如可能看到class : SpatRaster dimensions : 1000, 2000, 1 (nrow, ncol, nlyr) resolution : 0.01, 0.01 (x, y) extent : 100, 120, 30, 40 (xmin, xmax, ymin, ymax) coord. ref. : WGS 84 source : file.tif看到nlyr为1说明这是一个单波段栅格。如果显示nlyr大于1比如是3说明这个文件里包含多个波段后续计算时就需要决定是分别统计还是只取某一个波段。这个细节很容易被忽略但一旦忽略统计结果就会变得很不可控。3.4 用global函数计算单个栅格的统计值terra包中提取统计值的核心函数是global()。它的基础用法是global(r1, fun min, na.rm TRUE)这段代码会返回一个数据框行是图层列是对应统计值。如果你想同时计算多项统计值可以只调用一次global(r1, fun c(min, max, mean, sd), na.rm TRUE)我看到过很多人写代码时忘记加na.rm TRUE结果统计值全是NA。这里必须强调一下栅格数据中经常存在NoData值有些是用NA表示有些是用特定值如-9999表示如果在计算统计值时不做任何处理默认情况下很多函数会把NA传递给结果最后你拿到的就是一堆空值。不过“去掉NA”和“把-9999也当NA”是两回事。terra的global()默认只认识NA如果你的文件里用-9999表示无效值需要先做一步替换否则-9999会被当真值参与统计r1 - subst(r1, -9999, NA)这一步看似简单却是整个流程里最容易出问题的地方。我曾经接过一批高程数据文件里用的是-32768表示无效一开始没处理就直接统计标准差被拉得特别大后来排查了半天才发现是无效值混进了统计。4. 批量提取所有栅格的四项统计值单张栅格搞定后批量处理就是水到渠成的事核心思路是先准备一个空表然后循环读取每个文件计算统计值填进表里最后保存。4.1 拿到全部统计值后应处理的多波段问题批量处理时很多人会直接用rast(files[i])读取文件然后交给global()计算。对于单波段文件这个流程没问题。但如果是多波段文件比如一个tif里包含多个时相的栅格global()返回的结果会有多行直接把结果塞进一个以文件为单位的结果表就会产生行数不匹配的问题。我通常的处理方式是先判断一下波段数r_i - rast(files[i]) if (nlyr(r_i) 1) { r_i - r_i[[1]] # 取第一个波段或根据需要选择 }如果是做时间序列分析拍脑袋取第一个波段并不严谨。我会在读取之前就明确好自己的想法是需要逐波段统计还是只需要某一波段。如果只是对“每个文件”汇总一个值那最好确保所有文件的波段结构一致并在循环里统一取同一波段。4.2 完整批量循环代码下面这段代码是我实际项目中简化后的版本结构清晰可以直接套用到自己的数据上library(terra) folder_path - D:/R_work/raster_data files - list.files(folder_path, pattern \\.tif$, full.names TRUE) results - data.frame( file basename(files), min NA, max NA, mean NA, sd NA, stringsAsFactors FALSE ) for (i in seq_along(files)) { r_i - rast(files[i]) if (nlyr(r_i) 1) { r_i - r_i[[1]] } r_i - subst(r_i, -9999, NA) stat_names - c(min, max, mean, sd) stats_vec - sapply(stat_names, function(f) { global(r_i, fun f, na.rm TRUE)[1, 1] }) results[i, c(min, max, mean, sd)] - stats_vec cat(已完成, i, /, length(files), :, basename(files[i]), \n) } results这里我用了sapply()逐个计算四个统计值是因为在实际操作中global()一次性计算多个统计值时有些函数组合在某些版本的terra里会出现结果列名不够直观的情况。逐个计算虽然多写一行但结果更可控。循环里加了一行cat()输出进度。可能你会觉得这只是锦上添花但我在处理上百个大文件时最怕的就是运行到一半看起来像卡死了其实还在跑。输出一个当前文件名至少能知道程序卡在哪个文件上方便排查。4.3 用lapply更优雅地替代for循环有R基础的人可能会问这种逐文件批量操作是不是用lapply()更符合R的习惯确实可以。下面这种写法更简洁process_one - function(file_path) { r_i - rast(file_path) if (nlyr(r_i) 1) r_i - r_i[[1]] r_i - subst(r_i, -9999, NA) data.frame( file basename(file_path), min global(r_i, min, na.rm TRUE)[1, 1], max global(r_i, max, na.rm TRUE)[1, 1], mean global(r_i, mean, na.rm TRUE)[1, 1], sd global(r_i, sd, na.rm TRUE)[1, 1] ) } results_list - lapply(files, process_one) results - do.call(rbind, results_list)两种写法的效果是一样的。for循环在代码里更容易加进度提示和断点调试新手看明白的概率更高lapply适合已经熟练的人写起来短也确实更符合R的向量化风格。我自己早期写函数时比较喜欢lapply后来数据量大、需要加各种容错处理和日志输出时还是回到了for循环。这个没有绝对的优劣选自己习惯的就行。4.4 结果保存到本地CSV文件运行完批量计算后结果已经是一个data.frame了我通常习惯先看一眼前几行确认数据再写盘head(results) write.csv(results, raster_summary_stats.csv, row.names FALSE)这里有一个很实际的经验如果你用RStudio或Notepad打开CSV文件结果通常正常但如果直接用Excel双击打开中文文件名或中文路径可能会导致乱码。因为Excel在Windows下默认使用本地编码而R的write.csv()默认写出UTF-8编码。要避免这个问题可以在写出时指定带BOM的UTF-8编码con - file(raster_summary_stats.csv, encoding UTF-8) write.csv(results, con, row.names FALSE)或者直接输出一个没有中文列名的CSV这样绝大多数环境都不会乱码。另一个更省事的方法是保存成R自带格式saveRDS(results, raster_summary_stats.rds)下次用readRDS()一行读回来格式不会变形适合继续在R里做分析用。5. 实际操作中的常见问题与排查技巧说了这么多真正跑数据的时候还是会遇到一堆问题。我把这几年最常踩的坑集中整理一下方便你对照排查。5.1 统计值出现NA或异常值如果某个文件的统计值突然变成NA优先检查三件事文件路径是否读取成功、是否使用了错误的无效值替换、该区域是否真的全部是NoData。有些极端情况下比如裁剪后的栅格在某个范围内全是NoDataglobal()计算mean时即使设置了na.rm TRUE也可能返回NA因为有效像元数量为0。这种结果是有意义的代表这个文件在指定区域内没有有效数据。我建议在结果表里同时记录一个有效像元数防止被NA误导valid_n - global(r_i, notNA, na.rm TRUE)[1, 1]terra里的global()也可以计算notNA这种函数或者你自己统计一下非NA数量。这样一来当你看到均值为NA时就能知道这个栅格是否“空得彻底”。5.2 最小值和最大值看起来不对劲如果标准差没有问题但最小值和最大值远超出研究区域的合理范围那大概率是无效值没清理干净。前面提到的-9999、-999、-32768只是几种常见情况还有的影像会把背景值设为0但0可能是水体或裸地也可能真是无效值这需要结合具体数据含义来判断。我现在的习惯是在批量计算前先把每个文件的直方图信息打印出来看看数值分布是否连续。用terra快速获取分位数也能做到quantile(r_i, probs c(0.001, 0.01, 0.99, 0.999), na.rm TRUE)如果0.1%分位数附近突然出现一个非常突兀的数值十有八九是无效值混进去了。5.3 栅格之间分辨率和范围不一致批量处理时不同tif文件的分辨率、范围很可能不一样。比如某个文件夹里的tif有的覆盖整个研究区有的只覆盖一个县分辨率一个是10米一个是30米。如果你只是分别统计每个文件自己的内部统计值那不一样也不影响计算。但你如果需要把这几张栅格放在一起比较平均值那就要非常小心因为不同分辨率下像元数量不同简单平均并不等价于真正意义上的空间平均。范围不一致则会影响“同一区域”的比较两个文件如果范围不同平均值就代表的是不同区域的数值不能直接放到同一个时间序列里去解读。如果你需要统一到同一网格后再比较可以在R里先做重采样或提取公共区域。但这个过程涉及插值算法选择会比单纯统计复杂很多。建议先问清楚自己的业务需求到底只是看每个文件自己的状态还是要做跨文件严格比较。5.4 文件名排序和正则匹配问题list.files()返回的文件顺序在Windows和Linux下可能不一样而且如果没有零填充文件名排序还会出现类似“1、10、11、2”的序列。比如文件名是2020.tif、2021.tif这样还好如果是raster_1.tif到raster_100.tif单纯按字符串排序会得到raster_1、raster_10、raster_100这种乱序。解决办法是在读取前把文件名里的数字编号提取出来再按数字排序nums - as.numeric(sub(.*?([0-9])\\.tif$, \\1, files)) order_idx - order(nums) files - files[order_idx]如果文件名里自带年份、月份、日期这一步就特别实用。我通常会顺手把日期信息提取出来加进结果表里后面做时间序列分析就不用重新匹配了。5.5 大数据量处理时的内存和耗时控制处理非常大的栅格文件时最容易出现内存爆炸。虽然terra的global()已经按块读取但如果你在读取后进行subst()等操作可能还是会增加内存负担。建议开头就设置好terra的全局内存参数terraOptions(memfrac 0.6, progress 10)memfrac表示R可用的最大内存比例一般设在0.5到0.8之间progress 10表示每处理10%时输出一次进度信息。这样设置之后脚本在长时间运行时能更明确地看到进度也不会因为一次调用吃掉太多内存。如果文件实在太大比如单个几十GB还要考虑分块读取。terra里的readStart()、readValues()可以手动控制分块但一般做统计值用不到这么细global()内部已经优化过。6. 把统计结果接进后续分析流程批量提取完四个统计值之后很多人的流程就结束了但我还想多说一点这些统计值在数据分析流程里的价值往往比你看到的表格更大。6.1 根据文件名提取时间信息如果文件名规范比如NDVI_2020_07_15.tif可以直接在R里用正则表达式提取时间并生成year和month列time_info - regmatches(results$file, regexpr([0-9]{4}_[0-9]{2}, results$file)) results$year - as.numeric(substr(time_info, 1, 4)) results$month - as.numeric(substr(time_info, 6, 7))有了时间和统计值之后画出时间序列就很方便plot(results$year, results$mean, type b, xlab 年份, ylab 年均值, main NDVI年平均值变化趋势)这一步我非常推荐因为只看表格很难看出趋势画图后数据异常和长期变化一眼就能发现。6.2 按实验分组做后续比较如果这批栅格数据对应不同的实验处理比如对照区、实验区那么把处理名称从文件路径中提取出来结合统计值做分组比较是下一步常用操作。你可以使用tapply()或者tidyverse风格的group_by()函数汇总。这实际上体现了在R里做完整流程的好处数据提取、数据整理、统计检验、可视化全部能在同一个工具里完成不需要在不同软件之间频繁切换。6.3 验证统计结果的准确性批量脚本跑完之后验证非常关键特别是在做正式分析之前。我通常会在所有文件里抽2到3个文件和ArcGIS/ QGIS里的栅格属性统计结果对比一下。对比时特别要注意“标准差”的算法差异。ArcGIS里某些工具计算的标准差可能使用总体标准差即除以n而R基础包sd()函数默认使用样本标准差即除以n-1。当像元数量特别多时两者的差异很小几乎可以忽略但在算法验证时还可能出现个位数的差异让不熟悉的人误以为代码算错了。terra的global()里在计算sd时使用的统计公式也可以查一下帮助文档如果业务上有严格要求最好统一口径。我个人遇到过的情况是栅格像元数量动辄几百万甚至几千万样本标准差和总体标准差的差异通常在小数点后第三位才体现出来所以大多数分析场景下不用过度纠结。但如果涉及论文里的方法细节还是建议把统一标准写清楚。最后再分享一个实用习惯批量提取统计值这种事看起来是简单的土办法但用得好了能帮你在数据检查阶段就发现很多潜在问题。我现在拿到一批新栅格数据时的标准流程是先跑一遍批量统计生成结果表然后立刻画一张平均值的排序图或者时间序列图。如果哪张图的均值突然跳变我会第一时间回去看原始文件而不是等分析做到一半才发现数据有问题。这也是我写这个R脚本最想强调的一点自动化提取统计值不只是为了省事更是为了让数据质量检查变得可重复、可追踪。你不需要每次都用肉眼去翻一个图层表格和图会让你在几分钟内找到异常。希望这份记录和代码对你的工作有帮助也欢迎你在实际使用中根据自己的文件命名和业务需求做调整。