一个做草地遥感或生态模型的人拿到标题里这种“NPP 草原中国土木基1981-1990 年R1”的栅格数据时第一反应通常不是直接跑分析而是先弄清楚这个“R1”到底代表什么以及这十年的草原净初级生产力数据拿出来之后能算出哪些真正有说服力的结论。这个数据集我实际处理过里头坑不少单位、投影、波段含义任何一个环节理解偏了后面全盘皆错。这篇就把我自己的处理流程、代码和踩坑记录整理出来给同样要碰这套数据的人一个参考。1. 数据集本身1981-1990年中国草原NPP到底是什么1.1 NPP不是“长得多快”而是“固定了多少碳”净初级生产力Net Primary ProductivityNPP这个概念外行听起来像“草原长得多茂盛”但做生态的人心里清楚它真正衡量的是绿色植物在单位时间和单位面积内通过光合作用固定的有机碳总量减去自身呼吸消耗之后剩下的那部分。换句话说NPP回答的是“这片草地一年下来到底往生态系统里净存了多少碳”。单位通常是g C/m²/年也就是每平方米每年固定的碳克数。这个数值对草原生态研究之所以关键是因为草原不像森林那样有大量地上生物量积累它的碳主要存在地下根系和土壤里而且草原对气候波动极其敏感。1981到1990这十年恰好是中国北方草原区一个比较特殊的阶段既有气候的年度波动也有放牧压力变化这十年的NPP数据可以作为后续三十年变化的基线。我在实际项目中常拿它当“基准期”用来对比2000年之后MODIS时代草原NPP的升降幅度。1.2 R1版本一个栅格波段还是一个数据发布版这个“R1”的歧义我见了不止一次。在不同的数据源里它可能是两个完全不同的东西。第一种情况R1指的是“Release 1”也就是第一版发布数据。很多长时序植被生产力数据集会分版本发布比如GLASS、GIMMS、或一些区域再分析产品第一版和第二版之间的算法差异、输入气象数据不同会导致同一区域同一年的NPP数值有明显差别。如果你手上同时有R1和R2两个版本先别急着拿R2填充R1的缺失值两个版本之间的数值可比性需要先做散点回归检验。第二种情况R1是栅格文件里的第一个波段。遥感数据经常把多年结果打包成一个多波段文件比如一个GeoTIFF里面有10个波段分别对应1981到1990年的年均NPP那么波段1就是1981年波段10就是1990年。标题里的“R1”如果出现在文件名后缀多半是波段标记。我处理过一套数据文件名是“NPP_China_Grassland_1981-1990_R1.tif”里面其实只有一个波段R1纯粹是标识符表示这是该数据集的第一条记录。判断方法是直接读文件元数据。用R或Python打开栅格后看nbands如果只有一个波段但文件名有R1那它大概率是版本标识如果有十个波段那R1多半指第一个波段。千万别靠猜我见过有人把十波段数据错当成单一年份用统计结果偏到没法解释。2. 数据拿到手之后从文件到可用栅格的处理流程2.1 先认清楚文件格式和投影这类中国区域历史NPP数据常见格式有GeoTIFF、IMGERDAS、NetCDF偶尔也有以ASCII Grid形式分发的。第一步永远是读元数据搞清楚三件事投影坐标系、像元大小、无效值设定。中国区域的栅格数据投影可能是WGS84经纬度、Albers等面积投影、Lambert conformal conic或者UTM分带。不同投影下像元代表的实际地面面积差异很大。在WGS84经纬度投影下一个0.1度的像元在北方草原区大约北纬40度到50度对应的东西向距离大约8.5到6.5公里左右面积从70多平方公里到90多平方公里不等不能直接用像元面积乘统一值算总量。如果要做区域汇总建议先投影转换到Albers等面积投影中国常用中央经线105°E标准纬线25°N和47°N确保每个像元面积一致。检查无效值也很关键。很多NPP产品把水体、沙漠、裸岩、城市区域设为无效值编码可能是-9999、-3.402823e38、255或者NaN。如果一开始没把无效值设成NA后面求均值的时候这些极端负数会把结果拉得完全失真。我在第一次处理一套ASCII Grid数据时忘了看头文件里的NODATA_value结果全区平均NPP算出来是负的排查了半天才发现是沙漠区的大片-9999没过滤。2.2 单位换算和异常值过滤别偷懒NPP产品的单位五花八门。有的是g C/m²/年有的是kg C/m²/年还有的是g C/m²/月甚至有的模型输出是g DM/m²/年干物质不是碳。干物质转碳要乘0.45到0.5的系数月值转年值要把十二个月加起来而不是取平均。这些换算不做对数值完全没意义。我习惯在数据处理流程一开始就把单位统一成g C/m²/年并且做一次频率直方图检查。草地NPP在中国的典型范围干旱半干旱区大约在50到400 g C/m²/年典型草原和草甸草原可以到200到500极端湿润的草甸局部能超过600但超过1000的像元在自然草地里非常罕见。如果直方图上出现大量几千甚至上万的值先别急着欢呼“草原碳汇巨大”多半是单位看错了或者有异常像元混进来。实际处理时我会用分位数截断而不是硬阈值。比如先把小于0的值设为NA然后计算全区的第99.5百分位数把它作为可疑偏高值检查阈值对比对应位置的原始数据确认真实性。有些模型在灌溉农田上会输出极高NPP如果研究区里混了农田像元不排除是真实值需要结合土地覆盖数据掩膜。3. 用R语言做区域统计和时间趋势一套能直接抄的代码3.1 区域平均NPP提取拿到单年或多年NPP栅格后最常规的操作是提取某个区域省界、流域、草原类型区的平均NPP。我用R的terra包老一点的版本是raster包来做代码核心思路就三步读栅格、设无效值为NA、用区域矢量做zonal统计。下面这段代码以“十年单波段、R1标识”的数据为例直接提取内蒙古草原区的平均NPPlibrary(terra) # 读取1981-1990十年平均NPP栅格假设已经是g C/m2/yr r - rast(NPP_China_Grassland_1981-1990_R1.tif) # 查看基本信息 print(r) # 如果有无效值例如-9999需要替换为NA r[r -9999] - NA # 读取中国草原区矢量边界自己准备例如内蒙古典型草原区 grass - vect(typical_steppe.shp) # 确保栅格和矢量坐标系一致 if (crs(r) ! crs(grass)) { r - project(r, crs(grass)) } # 提取区域平均NPP region_vals - extract(r, grass, fun mean, na.rm TRUE, weights TRUE) print(region_vals)注意这里我加了weights TRUE这对于经纬度投影下的栅格很重要相当于按像元实际面积加权避免高纬度像元被等权平均。如果栅格本身已经是等面积投影weights加不加差别不大但加上没坏处。如果数据集是十年十个波段的单文件提取每年区域均值的代码稍微改一下r_multi - rast(NPP_China_Grassland_1981-1990_10bands.tif) n - nlyr(r_multi) year_means - numeric(n) for (i in 1:n) { layer_i - r_multi[[i]] layer_i[layer_i -9999] - NA year_means[i] - terra::extract(layer_i, grass, fun mean, na.rm TRUE, weights TRUE)[[1]] } names(year_means) - 1981:1990 print(year_means)这段代码我在不同数据集上跑过核心就是循环取波段、过滤无效值、加权统计。慢是慢一点但结果可靠。3.2 十年趋势Mann-Kendall检验与Sen斜率十年NPP数据光看平均值没意思能支撑文章的关键是趋势这十年草原生产力是在上升、下降还是基本稳定。样本数量只有10个用普通线性回归也可以但栅格逐像元做回归时异常值和数据非正态性会有影响。我偏向用Mann-Kendall趋势检验配合Sens slope估计这也是水文气象学界做长时序趋势检测的标配组合。Mann-Kendall检验不要求数据正态分布对异常值不敏感适合环境数据。Sens slope是取两两配对点斜率的中位数比最小二乘回归更稳。对单点时间序列做M-K检验R里有现成的包比如trend包library(trend) # 假设vec是某个像元1981-1990的NPP序列 vec - c(187.2, 201.5, 176.3, 208.9, 195.1, 213.6, 188.4, 224.7, 209.3, 231.8) mk_test - mk.test(vec) print(mk_test) # 计算Sen斜率 sen - sens.slope(vec) print(sen)mk.test输出里的pvalue和S统计量是核心p小于0.05认为趋势显著。sens.slope返回的估计值就是年均变化速率单位是g C/m²/年/年。比如Sen斜率为2.3意味着这十年平均每年NPP增加2.3 g C/m²。如果要全栅格逐像元算趋势就是对该像元的十个时间值做一次mk.test循环或使用apply。但要注意全栅格几十万个像元跑下来速度会比较慢我建议在提取到区域统计值之后先对区域均值序列做趋势分析再考虑逐像元分析避免一上来就跑全图。# 对区域均值序列做趋势检验 region_year - c(192.1, 198.4, 185.2, 204.7, 202.3, 211.8, 195.6, 208.2, 216.4, 223.1) mk.test(region_year) sens.slope(region_year)结果通常显示内蒙古典型草原区在1981到1990年间NPP呈小幅波动上升但很多区域p值不显著这与那个年代的降水波动大有关。真正显著的趋势在2000年后更常见。4. 常见问题与排查实录4.1 栅格范围对不齐怎么办历史NPP数据最头疼的问题之一是和现势遥感数据范围对不齐。1981-1990的数据经常是老坐标系比如北京54或西安80投影而2000年后的数据是WGS84或CGCS2000叠加在一起错位几百米到几公里很常见。我的处理习惯是统一投影到一个基准然后做重采样。具体用terra的project函数指定目标分辨率和重采样方法。对NPP这种连续变量用bilinear双线性插值比nearest neighbor更合理但注意重采样会平滑掉一部分极值。如果后续要做两个时期栅格的逐像元差值重采样方法不一致会导致虚假的边界变化务必保持两个栅格采用同一套处理参数。r_base - rast(MODIS_NPP_2001_base.tif) r_1980s - rast(NPP_China_Grassland_1981-1990_R1.tif) # 统一到同一CRS和网格 r_1980s_re - project(r_1980s, crs(r_base), res res(r_base), method bilinear)对齐之后还要检查像元的起始坐标是否完全一致。有时候两个栅格的分辨率相同但原点差了半个像元看起来对齐了实际没对齐。检查方法就是把两个栅格的extent打印出来对比origin要一致。4.2 为什么算出来的NPP数值偏大或偏小数值明显偏离文献范围通常逃不出下面几个原因一是单位换算少了系数。干物质和碳之间的转换很多人记成0.5实际要看植被类型草本植物含碳率大约在0.45左右灌木略高。如果源数据本身就是干物质直接当碳用结果会高估一倍多。二是时间聚合方式错了。如果源数据给的是每月NPP全年值应该是1月到12月求和而不是平均。拿月平均当成年均值NPP会低了一个数量级。三是掩膜范围没控制好。如果研究区矢量包含了农田或者城市边缘那统计出来的“草原NPP”其实混合了非草地像元。我在处理当雄和锡林郭勒的数据时都遇到过类似问题最后都是靠严格掩膜纯草地像元解决的用同时期的土地覆盖产品筛选出草地类别再和NPP栅格相乘。四是无效值没有彻底滤除。这个前面说过-9999和NaN这种要分清。有些GeoTIFF把无效值存在内部掩膜里直接读数值看不到需要额外读mask图层。我给一个排查清单现象可能原因排查方法全区均值是负数无效值未过滤查最小值分布看是否大量-9999数值比文献高1倍干物质没乘碳系数查数据说明单位是DM还是C数值比文献低一个量级月值当年均值确认是求和还是平均空间分布有大片0值掩膜范围不对叠加土地覆盖类型检查和同区域文献趋势相反投影或重采样偏移检查像元和边界是否对齐4.3 R1和其他版本混用会踩什么坑混用不同版本的数据需要格外小心。有一年我为了补一个1987年的缺失数据直接用另一个数据集的R2版本对应年份填充结果发现拼接后的时间序列在1987年出现一个明显的跳变区域的NPP从190突然跳到245然后又回落。最初怀疑是气候异常后来检查才发现是两套数据集的算法不一致R1用的是旧版光能利用率模型R2改进了水分胁迫系数导致湿润年份的模拟值系统性偏高。从那以后我给自己定了个规矩任何时间序列拼接都要先在重叠年份做双数据集对比。如果两个版本在同一个研究区同一年的像元值散点接近1:1线才允许互相填充如果系统性偏移超过10%就得做偏差校正或者干脆放弃拼接改用单一版本保证一致性。还有一种容易忽略的情况R1可能不是年份标识而是区域标识。比如“China North R1”和“China South R2”分别表示北方区第一块和南方区第二块。这时候如果你默认R1就是第一年整个区域都会张冠李戴。拿到文件先看配套的README或者xml元数据别急着开跑。4.4 几个我踩过的坑和现在的处理习惯第一个坑用terra读老IMG格式时坐标参考信息丢失。很多老IMG文件内嵌的投影信息不完整读取后crs显示为空。这时候不能直接叠加矢量必须手动给栅格设置投影。我在处理一批80年代内蒙古NPP数据时就是因为没设置投影后续所有空间统计全错返工了两天。设置投影的代码很简单crs(r_old) - EPSG:4326 # 如果确认数据是经纬度关键是你要确认真实投影。最稳妥的办法是找到原始数据的投影字符串从元数据文件里复制过来而不是凭印象设。第二个坑边界区的像元归属。用矢量边界做zonal统计时默认情况下部分像素可能落在边界外导致统计偏差。我的习惯是用terra的extract函数加weightsTRUE做加权平均或者先栅格化边界再裁剪出完整像元。边界影响对草原区来说通常不大但如果研究区很小比如只研究某一块典型草原边界一个像元的差异就会影响结果。第三个坑R语言里内存不够。全中国范围0.05度的NPP栅格十个波段叠起来内存占用不小。老笔记本跑起来很容易卡死。我现在的做法是分块处理用terra的aggregate先把分辨率降低或者用blocks分块循环读取统计。实在不行就转成NetCDF用ncdf4包分块读。第四个坑把“十年平均”直接当成“每一年”。很多数据集的标题看着像1981-1990实际上给的是这十年的平均态而不是十个年份分别的值。如果只有单波段那只能代表十年平均NPP不适合做趋势分析。要区分这种情况最简单的方法就是看波段数。单波段就没法算Mann-Kendall只能作为基准期均值和其他时期做差值。我现在的处理习惯拿到任何NPP数据后第一件事是写一个三行代码的快速体检print(r) summary(values(r)) hist(values(r), breaks 100)看一眼最小值、最大值、直方图形状心里就有数了。最大值几万或者最小值负几千的数据基本都有问题先解决再谈分析。这个习惯帮我避开了至少一半的数据坑。5. 我个人在这些年实践中的体会做草原NPP数据分析最大的感受是数据质量把控比分析算法本身更重要。算法再花哨输入数据的单位、投影、无效值任何一个环节出错输出全是垃圾。1981-1990这段历史时期的NPP数据算得上是国内草地生产力研究的一个基准数据源把它和2000年后的遥感NPP产品做对比能看出近几十年来草原生产力的变化趋势对生态修复评估和碳收支估算都有实际参考价值。但前提条件就是前面说的那些版本含义搞清楚、单位统一、无效值滤干净、投影对齐。如果你手头正好拿到一批类似命名的数据建议先花半小时读元数据再用上面提到的体检代码看一眼分布最后才进入正式分析。别看这半小时好像耽误了时间实际是在帮你避免后面几天的返工。这套流程我用下来无论是写论文还是做工程项目都省心得多。