简介2005-2021年中国西北地区新疆维吾尔自治区、青海省、甘肃省、内蒙古自治区和宁夏回族自治区1000米分辨率年NPP数据产品基于MOD17A3HGF.006原始产品重采样至WGS84坐标单位为g·C/m²。该年度NPP由全年8天净光合作用PSN合成反映植被总初级生产力GPP与维持呼吸MR之差。数据已完成无效值清洗可直接用于区域碳循环、植被生产力时空变化、土地退化及生态修复效果评估等研究适合生态遥感、GIS分析与GEE平台的科研人员与高校师生。压缩包共79个文件、约150.09MB包括19个逐年TIF栅格、19个TFW坐标配准文件、38个XML元数据说明及3个OVR金字塔另附NPP_mean均值栅格便于快速掌握总体水平。TIF可直接加载至ArcGIS/QGISTFW保证地理坐标精确匹配XML记录完整处理信息数据组织方式便于批量制图与质量检查。目前已有870人学习下载适合作为西北干旱区长时间序列NPP研究的参考数据也可用于验证模型模拟、开展趋势分析及空间格局对比。1. 打开这个 1000m 的 NPP 栅格之前先搞清楚这三件事从某课题组或数据分发站点拿到「西北地区_NPP_MOD17A3HGF_1000m_2005-2021.tif」第一反应往往是丢进渲染器、拉个渐变配色看个热闹。NPP净初级生产力是植被光合作用净固定的有机碳量MOD17A3HGF 是 MODIS 的年度 NPP 产品1000m 分辨率2005-2021 年逐年一景覆盖西北干旱与半干旱区。它常被用来做荒漠化评估、草地退化监测和区域碳源汇分析。但我要先泼一盆冷水这个产品里从单位到投影埋着四个暗坑直接打开第一张图的解读很可能是错的。这篇笔记就冲着「把这个 tif 变成一套可信趋势图」来写新手能照步骤跑通熟手可以跳过基础直奔避坑清单。2. 把 MOD17A3HGF 参数讲透算法来源、单位与栅格语义2.1 光能利用率模型NDVI 和气象数据如何算出一个年度的 NPPMOD17 系列不是直接用遥感反射率表出 NPP它依赖一个光能利用率模型。核心逻辑是植被每天吸收的光合有效辐射 APAR 乘以光能利用率 ε累计成总初级生产力 GPP再扣掉植被自养呼吸 R 得到 NPP。公式可以写成 GPP ε × FPAR × PARNPP ≈ GPP − R。MOD17 里 FPAR 来自 MODIS 的 NDVI/FPAR 产品PAR 来自气象再分析辐射场而 ε 不是常数它被最低温、水汽压差这些环境因子压住相当于「理想情况下能固定多少碳天气不给力时就打折」。干旱区植被对水分胁迫极度敏感这个 ε 的打折函数在西北地区的生态信号里比湿润区明显得多这也是为什么 NPP 比单纯看 NDVI 更能反映草地和荒漠生态系统的生产能力变化。MOD17A3HGF 的「HGF」是 V6 之后更新序列的标识。跟早期版本 MOD17A3 相比最关键的改动是不再依赖固定的生物群落查找表参数而是引入了动态辐射场和更新的 FPAR 算法让荒漠、草原这类低生产力、季节波动大的区域年际变化能真正反映在 NPP 数值上不至于被查表参数抹平。这也是为什么 2005-2021 这 17 年长序列拿出来做趋势分析比只看三五年的变化更有说服力。需要注意MOD17A3HGF 输出的是年累计 NPP没有逐月信号如果要做生长季物候变化得换 MOD17A2H 的 GPP 产品两者的数据流完全不一样。在西北地区用 1000m 是合理选择。这里的景观格局以荒漠、草原、绿洲过渡带为主斑块尺度通常是几公里到几十公里1000m 既能保留主要空间差异又不会像 250m 产品那样引入大量混合像元和噪声计算量也可控。做区域尺度碳收支或退化评估没必要上更高分辨率。2.2 一张表读全 MOD17A3HGF 的关键参数分辨率、单位、填充值、QC参数值说明空间分辨率1000mMODIS Sinusoidal 投影下约 1km × 1km时间粒度年累计值每年一个像元代表该年 NPP 总产量物理单位kg C/m²/yr每平方米每年固定的碳千克数尺度因子0.0001原始整数 DN × 0.0001 真值无效值/填充值原始 DN 30000换算后 3.0 判为异常需掩膜有效值范围0 ~ 约 2干旱半干旱区通常在 0 ~ 0.4局部绿洲可达 1 以上投影MODIS Sinusoidal正弦等积投影跟 WGS84 经纬度不重合波段Npp、Npp_QCNpp 是目标数据QC 用于质量筛选单位是第一个翻车点。原始 GeoTIFF 里存的大多是放大 10000 倍的整型数直接统计流域平均 NPP 会得到「几千」那是 DN 不是物理量。乘完 0.0001 得到 0.0x 到 1.x 的浮点才算对。我见过某开发者把 3000 多的均值写进报告被审稿人直接质疑数据量级离谱。填充值陷阱在文件里往往体现为 32761 或 65535如果不提前掩膜这些值会混进趋势计算一个 32761 的像元能把周边几十公里内的斜率拉偏。常见做法是把大于 30000 的 DN 全部判为 NoData这个上限比查表填 32761 更稳妥因为卫星边界和水域处有时会出现中间态异常值。QC 波段我一般不会逐 bit 深磕西北区域大、像元多多数研究只保留「质量良好且非填充」的像元就够了。如果你做连续多年像元级趋势请一定把 QC 作为第二层掩膜先按 Npp 阈值去异常再按 QC 去抖动。拿到文件后的第一步不是算趋势而是先检查逐年均值和空间分布有没有断裂。MOD17A3HGF 计算时依赖 NDVI 时间序列平滑干旱区冬季裸土反射率变化剧烈如果物候平滑不过关FPAR 会出现像元级跳变导致某一年 NPP 异常低。先把 17 年的区域均值拉一条曲线出来眼测有没有离谱的断崖再往下一步分析这能省下后面大量的排查时间。3. 用云端平台批量导出西北地区 2005-2021 年 NPP脚本与参数3.1 用规则矩形代替行政边界先能跑通再说处理的第一步是把数据从云端平台拉下来。这里用规则矩形当作研究区经纬度范围大致定在 73°E-108°E、34°N-50°N能覆盖西北主体部分。实际项目如果要求按行政边界统计把对应 shp 上传到平台 Assets 替换掉矩形即可后面的逻辑完全不变。先跑通数据处理链路再换精细边界这是我在多期项目里的习惯。var roi ee.Geometry.Rectangle([73, 34, 108, 50]); var col ee.ImageCollection(MODIS/061/MOD17A3HGF) .filterDate(2005-01-01, 2022-01-01); print(年度影像数量, col.size()); // 正常应为17这段代码先构建影像集合再打印数量。MOD17A3HGF 在云端平台每年一景时间戳一般是当年 1 月 1 日所以筛选日期范围写成 2005-01-01 到 2022-01-01结束日期带 2022 年不会误收 2022 年那景又能保证把 2021 年包含进来。如果 print 出来的 size 不是 17说明产品有缺失或时间戳偏移先解决集合问题再继续不要急着导出。3.2 年份筛选与批量导出17 个 tif 一个循环确认集合完整后用一个客户端循环逐个年份导出。代码里直接按年份重新 filterDate 并取 first()逻辑直观任务之间互不依赖即使某年失败也不影响其他年份。for (var y 2005; y 2021; y) { var start y -01-01; var end (y 1) -01-01; var img ee.ImageCollection(MODIS/061/MOD17A3HGF) .filterDate(start, end) .first(); var npp img.select(Npp).multiply(0.0001).rename(NPP); var valid npp.gt(0).and(npp.lt(3.0)); npp npp.updateMask(valid).clip(roi); Export.image.toDrive({ image: npp, description: NPP_ y, folder: NPP_2005_2021, region: roi, scale: 1000, maxPixels: 1e13 }); }参数说明multiply(0.0001) 必须放在 select 之后先把 Npp 波段取出来再转成物理量。valid 掩膜用了 0 到 3.0 的阈值换算后大于 3.0 的原始异常 DN 在这个阶段就被剔掉了。这里故意不传 crs让平台按影像自带的 MODIS Sinusoidal 投影分块导出。西北地区横跨多个 UTM 分带如果强行指定一个带号边缘区域变形会直接改变像元面积最后面积统计全靠重投影数据质量反而更差。提示如果某一年 first() 拿到空影像导出任务会在后台报错。先把第一章里的 size 确认做好再跑批量循环这是避免白等两小时的唯一办法。3.3 本地预处理三部曲重投影、去无效值、转浮点导出完成后每个 tif 是正弦投影的 1km 栅格下一步统一到 WGS84 经纬度。我一般用 GDAL 做批处理一条命令解决重投影和数据格式问题for f in NPP_*.tif; do gdalwarp -t_srs EPSG:4326 -r near -ot Float32 -dstnodata -9999 \ $f ${f%.tif}_wgs84.tif done-t_srs EPSG:4326 定义目标坐标系-r near 用最近邻重采样不插值出小数来污染物理量-ot Float32 输出 32 位浮点后续乘除运算不丢精度-dstnodata -9999 在输出边缘和无效区域写统一 NoData。如果需要等积投影做面积统计把 EPSG:4326 换成适合当地的 Albers 投影参数即可。重投影完成后再用 Python 补一道清洗因为 GDAL 参数只处理边缘的 NoData原文件里 30000 的异常 DN 经过乘法后可能变成像 3.2761 这样的残留值import numpy as np import rasterio with rasterio.open(NPP_2005_wgs84.tif) as src: data src.read(1) profile src.profile data np.where((data 0) | (data 3.0), np.nan, data) with rasterio.open(NPP_2005_clean.tif, w, **profile) as dst: dst.write(data, 1)这里把小于 0 和大于 3.0 的像元都置成 np.nan。西北干旱区有效 NPP 极少超过 1.5绿洲作物也难到 3.0超过这个数值的基本是水体残留或异常噪声清掉比留着更安全。清洗完的 clean 文件就是后续趋势分析的输入。4. 用 Theil-Sen 斜率和 Mann-Kendall 检验算 17 年趋势4.1 堆叠 17 年栅格逐像元算斜率的 Python 脚本趋势分析的最常见做法是逐像元计算斜率。先把 17 年 clean tif 读入一个三维数组import numpy as np import rasterio years list(range(2005, 2022)) stack [] for y in years: with rasterio.open(fNPP_{y}_clean.tif) as src: stack.append(src.read(1)) arr np.stack(stack, axis0) # shape: (17, rows, cols)arr 的第一维是年份后面两维是空间行列。接下来筛选有效像元17 年里至少有 15 年不是 NoData 的像元才参与趋势计算否则样本量太少算出的斜率没有意义。然后对每个有效像元调 scipy 的 theilslopesfrom scipy.stats import theilslopes valid_count np.sum(~np.isnan(arr), axis0) idx np.where(valid_count 15) slope np.full(arr.shape[1:], np.nan) intercept np.full(arr.shape[1:], np.nan) for r, c in zip(*idx): series arr[:, r, c] mask ~np.isnan(series) ys np.array(years, dtypefloat)[mask] s, b, lo, hi theilslopes(series[mask], ys[mask]) slope[r, c] s intercept[r, c] b这个双层循环纯粹为了可读性像元数量大时跑起来会很慢。优化方式很多先用常年裸地掩膜跳过 NPP 均值低于 0.01 的像元再把数组切块丢给 multiprocessing。Theil-Sen 的原理是对任意两个不同年份的像元值求斜率最后取所有斜率的中位数。它不像最小二乘那样被个别极端年份扯走对干旱区常出现的「某年因雪灾或极端干旱掉到 0.001」这类尖峰异常非常稳健。4.2 Mann-Kendall 显著性检验p0.05 才允许上色斜率大不代表趋势可信还要检验显著性。Mann-Kendall 是非参数检验方法不要求数据正态分布非常适合 NPP 这种经常出现偏态分布的生态指标import pymannkendall as mk pval np.full_like(slope, np.nan) for r, c in zip(*idx): series arr[:, r, c] mask ~np.isnan(series) ys np.array(years, dtypefloat)[mask] if len(ys) 10: continue result mk.original_test(series[mask]) pval[r, c] result.ppymannkendall 返回的 result 对象包含 trend 类别和 p 值。original_test 要求输入序列中不能有缺失值所以提前把 mask 后的年份重新排列。出图的核心逻辑是先掩掉不显著像元signif pval 0.05 trend_show np.where(signif, slope, np.nan)trend_show 里只保留通过显著性检验的像元其余都是 NaN。渲染时用 RdYlGn 配色负斜率显示红色正斜率显示绿色图例单位写成 kg C/m²/yr/yr。如果直接拿 slope 出图不加 p 值掩膜「整片西北红彤彤」的视觉效果会强烈误导判断以为全区都在退化实际上大部分像元年份间随机波动大斜率根本通不过检验。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 6)) im ax.imshow(trend_show, cmapRdYlGn, vmin-0.02, vmax0.02, extent[73, 108, 34, 50]) plt.colorbar(im, axax, labelNPP trend (kg C/m²/yr/yr)) plt.savefig(NPP_trend_2005_2021.png, dpi300, bbox_inchestight)vmin/vmax 的设定对出图影响很大这是很多第一次做趋势图的人会纠结的地方。我一般先用 np.nanpercentile(trend_show, [2, 98]) 看数据分布再根据研究目标微调。设成 ±0.02 时绿洲强趋势能显现荒漠像元如果斜率在 ±0.005 以内就保持浅色如果设成 ±0.1整张图会灰蒙蒙一片真实趋势反而被淹没。阈值没有绝对正确值但结论不能因为配色范围随意改变这一点务必写进图注。5. NPP 栅格处理避坑清单现象、原因、解决5.1 全图均值 3000 多单位陷阱现象用 gdalinfo 或 Python 打开 tif统计值几千渲染出来白花花一片完全看不出空间分异。原因读取的是原始 DN没有乘尺度因子 0.0001。解决在云端平台导出阶段就执行 multiply(0.0001)本地读取后也必须立刻检查数值量级。西北地区区域均值正常范围在 0.02 到 0.2 之间看到超过 1 的均值先查单位别急着分析。5.2 边界外一圈矩形数据clip 没生效现象按矩形导出后四角和水域有一圈非 NoData 值显示为 0 或 3.2761统计面积时被当成有效像元。原因正弦投影下 clip 的边界和输出 tif 的像元格网不对齐原始文件自带的外包矩形背景残留。解决重投影后统一再做一遍数值掩膜。gdalwarp 的 -dstnodata 只处理边缘不处理原文件里乘完尺度因子后的残留异常值必须用第一章的 rasterio 清洗逻辑再跑一次。5.3 某一年影像全是背景filterDate 拿到空图现象批量任务跑完某一年打开全是 NoData其余年份正常。原因GEE 里该年份影像的时间戳不是当年 1 月 1 日或者产品集合在特定年份有版本变体filterDate 加 first() 组合抓到了空影像。解决先在集合级别 print(col.size()) 确认数量如果缺年份改用按影像 ID 拼接年份的方式拉取对应影像绕过时间戳筛选的不确定性。批量导出前务必做这一步检查。5.4 Theil-Sen 结果一大片红色没做显著性检验现象趋势图里红色区域占比 90% 以上和同年降水数据反映的趋势完全对不上。原因干旱区大量像元 NPP 数值极低0 到 0.05 之间年份间随机波动大最小二乘或 Theil-Sen 都能算出负斜率但统计上根本不显著。解决按第四章方法先算 Mann-Kendall p 值p0.05 才上色。同时把色阶拉伸控制在 ±0.02 以内避免微小斜率被夸大。5.5 波动曲线突然断崖某年像元被 mask 掉了现象区域年均值曲线在 2009 年前后跳水检查单年栅格发现大片区域变成 NoData。原因极端干旱年份触发算法异常NPP 出现大范围负值或填充值预处理时统一阈值把大量像元删掉了。解决曲线统计时额外输出每年的有效像元覆盖率覆盖率低于 95% 时用前后年份线性插值补洞再画曲线否则曲线形态只是缺值的假象不是真实生态信号。这个覆盖率字段在第六章会继续用到。6. 把 17 年数据收进一张波动曲线区域均值的进阶画法趋势图看空间分布波动曲线看时间过程。把西北地区每年的区域均值提取出来加上标准差阴影能直观看到 NPP 年际变化和异常年份。这里有个关键技巧统计均值时只统计有效像元并同时报告有效覆盖率两个指标一起看才知道曲线是否可信。import numpy as np import rasterio means, stds, cover [], [], [] for y in years: with rasterio.open(fNPP_{y}_clean.tif) as src: data src.read(1) v data[np.isfinite(data) (data 0) (data 3.0)] means.append(np.mean(v)) stds.append(np.std(v)) cover.append(np.sum(np.isfinite(data)) / np.prod(data.shape)) import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 5)) ax.plot(years, means, color#2b6f44, linewidth2, label区域均值) ax.fill_between(years, np.array(means) - np.array(stds), np.array(means) np.array(stds), color#2b6f44, alpha0.2, label±1 标准差) ax.set_xlabel(年份) ax.set_ylabel(NPP (kg C/m²/yr)) ax.legend() plt.savefig(NPP_annual_curve.png, dpi300, bbox_inchestight)标准差阴影反映空间异质性干旱年份整体下降时阴影通常会收窄因为大量低值像元被统一压低。如果发现某年 std 突然变大多半是那一年的有效覆盖率异常需要回到预处理阶段排查。曲线画出来后可以叠加一个 5 年滑动平均线用来消除 ENSO 或极端气候事件的单年扰动让长期趋势更清楚。我第一次做西北 NPP 趋势时第一版图就是没扣单位、没做显著性检验被同组人指出「你这个结论是配色给出来的」——从那以后我把尺度因子、填充值、显著性三件事写死在处理流程最前面。这套习惯也送给你希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/20 3:06:52 · 阅读 1.2k NPP栅格预处理与趋势分析实战：MOD17A3HGF数据避坑指南
简介2005-2021年中国西北地区新疆维吾尔自治区、青海省、甘肃省、内蒙古自治区和宁夏回族自治区1000米分辨率年NPP数据产品基于MOD17A3HGF.006原始产品重采样至WGS84坐标单位为g·C/m²。该年度NPP由全年8天净光合作用PSN合成反映植被总初级生产力GPP与维持呼吸MR之差。数据已完成无效值清洗可直接用于区域碳循环、植被生产力时空变化、土地退化及生态修复效果评估等研究适合生态遥感、GIS分析与GEE平台的科研人员与高校师生。压缩包共79个文件、约150.09MB包括19个逐年TIF栅格、19个TFW坐标配准文件、38个XML元数据说明及3个OVR金字塔另附NPP_mean均值栅格便于快速掌握总体水平。TIF可直接加载至ArcGIS/QGISTFW保证地理坐标精确匹配XML记录完整处理信息数据组织方式便于批量制图与质量检查。目前已有870人学习下载适合作为西北干旱区长时间序列NPP研究的参考数据也可用于验证模型模拟、开展趋势分析及空间格局对比。1. 打开这个 1000m 的 NPP 栅格之前先搞清楚这三件事从某课题组或数据分发站点拿到「西北地区_NPP_MOD17A3HGF_1000m_2005-2021.tif」第一反应往往是丢进渲染器、拉个渐变配色看个热闹。NPP净初级生产力是植被光合作用净固定的有机碳量MOD17A3HGF 是 MODIS 的年度 NPP 产品1000m 分辨率2005-2021 年逐年一景覆盖西北干旱与半干旱区。它常被用来做荒漠化评估、草地退化监测和区域碳源汇分析。但我要先泼一盆冷水这个产品里从单位到投影埋着四个暗坑直接打开第一张图的解读很可能是错的。这篇笔记就冲着「把这个 tif 变成一套可信趋势图」来写新手能照步骤跑通熟手可以跳过基础直奔避坑清单。2. 把 MOD17A3HGF 参数讲透算法来源、单位与栅格语义2.1 光能利用率模型NDVI 和气象数据如何算出一个年度的 NPPMOD17 系列不是直接用遥感反射率表出 NPP它依赖一个光能利用率模型。核心逻辑是植被每天吸收的光合有效辐射 APAR 乘以光能利用率 ε累计成总初级生产力 GPP再扣掉植被自养呼吸 R 得到 NPP。公式可以写成 GPP ε × FPAR × PARNPP ≈ GPP − R。MOD17 里 FPAR 来自 MODIS 的 NDVI/FPAR 产品PAR 来自气象再分析辐射场而 ε 不是常数它被最低温、水汽压差这些环境因子压住相当于「理想情况下能固定多少碳天气不给力时就打折」。干旱区植被对水分胁迫极度敏感这个 ε 的打折函数在西北地区的生态信号里比湿润区明显得多这也是为什么 NPP 比单纯看 NDVI 更能反映草地和荒漠生态系统的生产能力变化。MOD17A3HGF 的「HGF」是 V6 之后更新序列的标识。跟早期版本 MOD17A3 相比最关键的改动是不再依赖固定的生物群落查找表参数而是引入了动态辐射场和更新的 FPAR 算法让荒漠、草原这类低生产力、季节波动大的区域年际变化能真正反映在 NPP 数值上不至于被查表参数抹平。这也是为什么 2005-2021 这 17 年长序列拿出来做趋势分析比只看三五年的变化更有说服力。需要注意MOD17A3HGF 输出的是年累计 NPP没有逐月信号如果要做生长季物候变化得换 MOD17A2H 的 GPP 产品两者的数据流完全不一样。在西北地区用 1000m 是合理选择。这里的景观格局以荒漠、草原、绿洲过渡带为主斑块尺度通常是几公里到几十公里1000m 既能保留主要空间差异又不会像 250m 产品那样引入大量混合像元和噪声计算量也可控。做区域尺度碳收支或退化评估没必要上更高分辨率。2.2 一张表读全 MOD17A3HGF 的关键参数分辨率、单位、填充值、QC参数值说明空间分辨率1000mMODIS Sinusoidal 投影下约 1km × 1km时间粒度年累计值每年一个像元代表该年 NPP 总产量物理单位kg C/m²/yr每平方米每年固定的碳千克数尺度因子0.0001原始整数 DN × 0.0001 真值无效值/填充值原始 DN 30000换算后 3.0 判为异常需掩膜有效值范围0 ~ 约 2干旱半干旱区通常在 0 ~ 0.4局部绿洲可达 1 以上投影MODIS Sinusoidal正弦等积投影跟 WGS84 经纬度不重合波段Npp、Npp_QCNpp 是目标数据QC 用于质量筛选单位是第一个翻车点。原始 GeoTIFF 里存的大多是放大 10000 倍的整型数直接统计流域平均 NPP 会得到「几千」那是 DN 不是物理量。乘完 0.0001 得到 0.0x 到 1.x 的浮点才算对。我见过某开发者把 3000 多的均值写进报告被审稿人直接质疑数据量级离谱。填充值陷阱在文件里往往体现为 32761 或 65535如果不提前掩膜这些值会混进趋势计算一个 32761 的像元能把周边几十公里内的斜率拉偏。常见做法是把大于 30000 的 DN 全部判为 NoData这个上限比查表填 32761 更稳妥因为卫星边界和水域处有时会出现中间态异常值。QC 波段我一般不会逐 bit 深磕西北区域大、像元多多数研究只保留「质量良好且非填充」的像元就够了。如果你做连续多年像元级趋势请一定把 QC 作为第二层掩膜先按 Npp 阈值去异常再按 QC 去抖动。拿到文件后的第一步不是算趋势而是先检查逐年均值和空间分布有没有断裂。MOD17A3HGF 计算时依赖 NDVI 时间序列平滑干旱区冬季裸土反射率变化剧烈如果物候平滑不过关FPAR 会出现像元级跳变导致某一年 NPP 异常低。先把 17 年的区域均值拉一条曲线出来眼测有没有离谱的断崖再往下一步分析这能省下后面大量的排查时间。3. 用云端平台批量导出西北地区 2005-2021 年 NPP脚本与参数3.1 用规则矩形代替行政边界先能跑通再说处理的第一步是把数据从云端平台拉下来。这里用规则矩形当作研究区经纬度范围大致定在 73°E-108°E、34°N-50°N能覆盖西北主体部分。实际项目如果要求按行政边界统计把对应 shp 上传到平台 Assets 替换掉矩形即可后面的逻辑完全不变。先跑通数据处理链路再换精细边界这是我在多期项目里的习惯。var roi ee.Geometry.Rectangle([73, 34, 108, 50]); var col ee.ImageCollection(MODIS/061/MOD17A3HGF) .filterDate(2005-01-01, 2022-01-01); print(年度影像数量, col.size()); // 正常应为17这段代码先构建影像集合再打印数量。MOD17A3HGF 在云端平台每年一景时间戳一般是当年 1 月 1 日所以筛选日期范围写成 2005-01-01 到 2022-01-01结束日期带 2022 年不会误收 2022 年那景又能保证把 2021 年包含进来。如果 print 出来的 size 不是 17说明产品有缺失或时间戳偏移先解决集合问题再继续不要急着导出。3.2 年份筛选与批量导出17 个 tif 一个循环确认集合完整后用一个客户端循环逐个年份导出。代码里直接按年份重新 filterDate 并取 first()逻辑直观任务之间互不依赖即使某年失败也不影响其他年份。for (var y 2005; y 2021; y) { var start y -01-01; var end (y 1) -01-01; var img ee.ImageCollection(MODIS/061/MOD17A3HGF) .filterDate(start, end) .first(); var npp img.select(Npp).multiply(0.0001).rename(NPP); var valid npp.gt(0).and(npp.lt(3.0)); npp npp.updateMask(valid).clip(roi); Export.image.toDrive({ image: npp, description: NPP_ y, folder: NPP_2005_2021, region: roi, scale: 1000, maxPixels: 1e13 }); }参数说明multiply(0.0001) 必须放在 select 之后先把 Npp 波段取出来再转成物理量。valid 掩膜用了 0 到 3.0 的阈值换算后大于 3.0 的原始异常 DN 在这个阶段就被剔掉了。这里故意不传 crs让平台按影像自带的 MODIS Sinusoidal 投影分块导出。西北地区横跨多个 UTM 分带如果强行指定一个带号边缘区域变形会直接改变像元面积最后面积统计全靠重投影数据质量反而更差。提示如果某一年 first() 拿到空影像导出任务会在后台报错。先把第一章里的 size 确认做好再跑批量循环这是避免白等两小时的唯一办法。3.3 本地预处理三部曲重投影、去无效值、转浮点导出完成后每个 tif 是正弦投影的 1km 栅格下一步统一到 WGS84 经纬度。我一般用 GDAL 做批处理一条命令解决重投影和数据格式问题for f in NPP_*.tif; do gdalwarp -t_srs EPSG:4326 -r near -ot Float32 -dstnodata -9999 \ $f ${f%.tif}_wgs84.tif done-t_srs EPSG:4326 定义目标坐标系-r near 用最近邻重采样不插值出小数来污染物理量-ot Float32 输出 32 位浮点后续乘除运算不丢精度-dstnodata -9999 在输出边缘和无效区域写统一 NoData。如果需要等积投影做面积统计把 EPSG:4326 换成适合当地的 Albers 投影参数即可。重投影完成后再用 Python 补一道清洗因为 GDAL 参数只处理边缘的 NoData原文件里 30000 的异常 DN 经过乘法后可能变成像 3.2761 这样的残留值import numpy as np import rasterio with rasterio.open(NPP_2005_wgs84.tif) as src: data src.read(1) profile src.profile data np.where((data 0) | (data 3.0), np.nan, data) with rasterio.open(NPP_2005_clean.tif, w, **profile) as dst: dst.write(data, 1)这里把小于 0 和大于 3.0 的像元都置成 np.nan。西北干旱区有效 NPP 极少超过 1.5绿洲作物也难到 3.0超过这个数值的基本是水体残留或异常噪声清掉比留着更安全。清洗完的 clean 文件就是后续趋势分析的输入。4. 用 Theil-Sen 斜率和 Mann-Kendall 检验算 17 年趋势4.1 堆叠 17 年栅格逐像元算斜率的 Python 脚本趋势分析的最常见做法是逐像元计算斜率。先把 17 年 clean tif 读入一个三维数组import numpy as np import rasterio years list(range(2005, 2022)) stack [] for y in years: with rasterio.open(fNPP_{y}_clean.tif) as src: stack.append(src.read(1)) arr np.stack(stack, axis0) # shape: (17, rows, cols)arr 的第一维是年份后面两维是空间行列。接下来筛选有效像元17 年里至少有 15 年不是 NoData 的像元才参与趋势计算否则样本量太少算出的斜率没有意义。然后对每个有效像元调 scipy 的 theilslopesfrom scipy.stats import theilslopes valid_count np.sum(~np.isnan(arr), axis0) idx np.where(valid_count 15) slope np.full(arr.shape[1:], np.nan) intercept np.full(arr.shape[1:], np.nan) for r, c in zip(*idx): series arr[:, r, c] mask ~np.isnan(series) ys np.array(years, dtypefloat)[mask] s, b, lo, hi theilslopes(series[mask], ys[mask]) slope[r, c] s intercept[r, c] b这个双层循环纯粹为了可读性像元数量大时跑起来会很慢。优化方式很多先用常年裸地掩膜跳过 NPP 均值低于 0.01 的像元再把数组切块丢给 multiprocessing。Theil-Sen 的原理是对任意两个不同年份的像元值求斜率最后取所有斜率的中位数。它不像最小二乘那样被个别极端年份扯走对干旱区常出现的「某年因雪灾或极端干旱掉到 0.001」这类尖峰异常非常稳健。4.2 Mann-Kendall 显著性检验p0.05 才允许上色斜率大不代表趋势可信还要检验显著性。Mann-Kendall 是非参数检验方法不要求数据正态分布非常适合 NPP 这种经常出现偏态分布的生态指标import pymannkendall as mk pval np.full_like(slope, np.nan) for r, c in zip(*idx): series arr[:, r, c] mask ~np.isnan(series) ys np.array(years, dtypefloat)[mask] if len(ys) 10: continue result mk.original_test(series[mask]) pval[r, c] result.ppymannkendall 返回的 result 对象包含 trend 类别和 p 值。original_test 要求输入序列中不能有缺失值所以提前把 mask 后的年份重新排列。出图的核心逻辑是先掩掉不显著像元signif pval 0.05 trend_show np.where(signif, slope, np.nan)trend_show 里只保留通过显著性检验的像元其余都是 NaN。渲染时用 RdYlGn 配色负斜率显示红色正斜率显示绿色图例单位写成 kg C/m²/yr/yr。如果直接拿 slope 出图不加 p 值掩膜「整片西北红彤彤」的视觉效果会强烈误导判断以为全区都在退化实际上大部分像元年份间随机波动大斜率根本通不过检验。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 6)) im ax.imshow(trend_show, cmapRdYlGn, vmin-0.02, vmax0.02, extent[73, 108, 34, 50]) plt.colorbar(im, axax, labelNPP trend (kg C/m²/yr/yr)) plt.savefig(NPP_trend_2005_2021.png, dpi300, bbox_inchestight)vmin/vmax 的设定对出图影响很大这是很多第一次做趋势图的人会纠结的地方。我一般先用 np.nanpercentile(trend_show, [2, 98]) 看数据分布再根据研究目标微调。设成 ±0.02 时绿洲强趋势能显现荒漠像元如果斜率在 ±0.005 以内就保持浅色如果设成 ±0.1整张图会灰蒙蒙一片真实趋势反而被淹没。阈值没有绝对正确值但结论不能因为配色范围随意改变这一点务必写进图注。5. NPP 栅格处理避坑清单现象、原因、解决5.1 全图均值 3000 多单位陷阱现象用 gdalinfo 或 Python 打开 tif统计值几千渲染出来白花花一片完全看不出空间分异。原因读取的是原始 DN没有乘尺度因子 0.0001。解决在云端平台导出阶段就执行 multiply(0.0001)本地读取后也必须立刻检查数值量级。西北地区区域均值正常范围在 0.02 到 0.2 之间看到超过 1 的均值先查单位别急着分析。5.2 边界外一圈矩形数据clip 没生效现象按矩形导出后四角和水域有一圈非 NoData 值显示为 0 或 3.2761统计面积时被当成有效像元。原因正弦投影下 clip 的边界和输出 tif 的像元格网不对齐原始文件自带的外包矩形背景残留。解决重投影后统一再做一遍数值掩膜。gdalwarp 的 -dstnodata 只处理边缘不处理原文件里乘完尺度因子后的残留异常值必须用第一章的 rasterio 清洗逻辑再跑一次。5.3 某一年影像全是背景filterDate 拿到空图现象批量任务跑完某一年打开全是 NoData其余年份正常。原因GEE 里该年份影像的时间戳不是当年 1 月 1 日或者产品集合在特定年份有版本变体filterDate 加 first() 组合抓到了空影像。解决先在集合级别 print(col.size()) 确认数量如果缺年份改用按影像 ID 拼接年份的方式拉取对应影像绕过时间戳筛选的不确定性。批量导出前务必做这一步检查。5.4 Theil-Sen 结果一大片红色没做显著性检验现象趋势图里红色区域占比 90% 以上和同年降水数据反映的趋势完全对不上。原因干旱区大量像元 NPP 数值极低0 到 0.05 之间年份间随机波动大最小二乘或 Theil-Sen 都能算出负斜率但统计上根本不显著。解决按第四章方法先算 Mann-Kendall p 值p0.05 才上色。同时把色阶拉伸控制在 ±0.02 以内避免微小斜率被夸大。5.5 波动曲线突然断崖某年像元被 mask 掉了现象区域年均值曲线在 2009 年前后跳水检查单年栅格发现大片区域变成 NoData。原因极端干旱年份触发算法异常NPP 出现大范围负值或填充值预处理时统一阈值把大量像元删掉了。解决曲线统计时额外输出每年的有效像元覆盖率覆盖率低于 95% 时用前后年份线性插值补洞再画曲线否则曲线形态只是缺值的假象不是真实生态信号。这个覆盖率字段在第六章会继续用到。6. 把 17 年数据收进一张波动曲线区域均值的进阶画法趋势图看空间分布波动曲线看时间过程。把西北地区每年的区域均值提取出来加上标准差阴影能直观看到 NPP 年际变化和异常年份。这里有个关键技巧统计均值时只统计有效像元并同时报告有效覆盖率两个指标一起看才知道曲线是否可信。import numpy as np import rasterio means, stds, cover [], [], [] for y in years: with rasterio.open(fNPP_{y}_clean.tif) as src: data src.read(1) v data[np.isfinite(data) (data 0) (data 3.0)] means.append(np.mean(v)) stds.append(np.std(v)) cover.append(np.sum(np.isfinite(data)) / np.prod(data.shape)) import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 5)) ax.plot(years, means, color#2b6f44, linewidth2, label区域均值) ax.fill_between(years, np.array(means) - np.array(stds), np.array(means) np.array(stds), color#2b6f44, alpha0.2, label±1 标准差) ax.set_xlabel(年份) ax.set_ylabel(NPP (kg C/m²/yr)) ax.legend() plt.savefig(NPP_annual_curve.png, dpi300, bbox_inchestight)标准差阴影反映空间异质性干旱年份整体下降时阴影通常会收窄因为大量低值像元被统一压低。如果发现某年 std 突然变大多半是那一年的有效覆盖率异常需要回到预处理阶段排查。曲线画出来后可以叠加一个 5 年滑动平均线用来消除 ENSO 或极端气候事件的单年扰动让长期趋势更清楚。我第一次做西北 NPP 趋势时第一版图就是没扣单位、没做显著性检验被同组人指出「你这个结论是配色给出来的」——从那以后我把尺度因子、填充值、显著性三件事写死在处理流程最前面。这套习惯也送给你希望帮到你。本文还有配套的精品资源点击获取