简介：这是一套基于Python实现的CASA模型净初级生产力（NPP）估算脚本，同时提供可直接使用的示例遥感影像数据，主要面向生态遥感、地理信息科学方向的初学者与研究人员，也可服务于区域碳汇评估或相关课程设计。资源包共19个文件，包含2个Python程序、14个GeoTIFF格式的影像文件、1个依赖说明文件等，整体大小约20.57MB；其中核心脚本casa_npp.py完成了从太阳辐射、植被指数等输入到植被吸收光合有效辐射（APAR）、光能利用率及NPP逐像元计算的完整流程，并配套了多个月份的输入数据与输出结果，便于使用者对照检验。目前已有54人浏览学习。代码结构清晰、变量命名规范，已通过实际数据验证，且不依赖复杂第三方库；使用时只需替换输入影像路径并调整温度系数、植被类型等关键参数，即可快速复现估算流程。此外，压缩包中还包含create_demo_data.py和requirements.txt，前者可辅助构建演示数据，后者便于一键配置运行环境，整体上手门槛较低，适合作为CASA模型入门学习和轻量级生产力估算的实用工具。
1. 用 Python 算 NPP：CASA 模型从公式到可跑脚本的完整闭环
遥感生态方向拿到 NDVI、气温、太阳辐射这类年度序列数据时，最常卡住的问题不是"该用什么模型"，而是"CASA 模型那么多中间量，怎么在代码里一步步落地成栅格"。这套基于 Python 的 CASA 模型 NPP 计算脚本，把光合有效辐射吸收比例（FPAR）、光能利用率（LUE）以及最终净初级生产力（NPP）的计算全部打通，还附带一组可直接试跑的示例遥感影像数据，解决的是从公式到 GeoTIFF 输出的最后一公里问题。适合有 Python 基础、想直接在真实数据上复现 NPP 的从业者和研究生，不用再从零拼装地理数据处理链路。
2. CASA 模型的计算骨架：FPAR、LUE 与 NPP 的三个中间量
CASA 属于光能利用率模型，核心理念是植被的净初级生产力由吸收的光合有效辐射（APAR）以及把这部分辐射转化为有机碳的效率（LUE）共同决定。写成公式就是：
NPP = APAR × LUE其中 APAR 又可以进一步拆成光合有效辐射（PAR）和 FPAR 的乘积，而 LUE 则是对最大光能利用率的温度、水分胁迫折减。整个模型看似简单，但每个中间量都有各自的输入数据和取值陷阱，下面按计算顺序拆开说。
2.1 FPAR：从 NDVI 到植被吸收的有效光合辐射
FPAR（Fraction of Photosynthetically Active Radiation）表示植被冠层吸收光合有效辐射的比例，取值范围一般在 0 到 0.95 之间。最常用的估算路径是归一化植被指数 NDVI，原理是 NDVI 对植被冠层叶面积和吸收效率敏感，可以直接用极值线性拉伸映射到 FPAR：
FPAR_NDVI = (NDVI - NDVI_min) / (NDVI_max - NDVI_min) × (FPAR_max - FPAR_min) + FPAR_min不过在实践中，单靠 NDVI 一条路径容易出现饱和问题，尤其是茂密植被区 NDVI 接近 0.9 以后增长乏力，FPAR 会被低估。所以很多实现会引入简单比值 SR（Simple Ratio）做第二条路径，最后取两者平均：
SR = (1 + NDVI) / (1 - NDVI) FPAR_SR = (SR - SR_min) / (SR_max - SR_min) × (FPAR_max - FPAR_min) + FPAR_min FPAR = (FPAR_NDVI + FPAR_SR) / 2这套脚本里默认参数为 NDVI_min=0.05、NDVI_max=0.95、SR_min=1.08、SR_max=5.0、FPAR_min=0.001、FPAR_max=0.95，这些是文献里最常见的经验取值。这里有一个容易犯的错：跑到茂密森林区域时如果发现 FPAR 常年顶到 0.94 以上，不要急着怀疑公式，先检查 NDVI_max 是不是被个别异常像元（比如云边或水体的高值噪声）拉大了，导致正常像元的拉伸区间被压缩。NDVI 极值建议用数据自身分位数去替代全局 min/max，脚本在配置里预留了这个开关。
用 numpy 可以很直接地实现这个双路径计算，向量化处理比逐像元 for 循环快几个量级：
import numpy as np def calc_fpar(ndvi, ndvi_min=0.05, ndvi_max=0.95, sr_min=1.08, sr_max=5.0, fpar_min=0.001, fpar_max=0.95): ndvi = np.clip(ndvi, 0.0, 1.0) # 路径 1：NDVI 线性拉伸 fpar_ndvi = (ndvi - ndvi_min) / (ndvi_max - ndvi_min) fpar_ndvi = fpar_ndvi * (fpar_max - fpar_min) + fpar_min # 路径 2：简单比值 SR 拉伸 sr = (1 + ndvi) / (1 - ndvi + 1e-10) fpar_sr = (sr - sr_min) / (sr_max - sr_min) fpar_sr = fpar_sr * (fpar_max - fpar_min) + fpar_min # 双路径平均，并裁剪到合法区间 fpar = np.clip((fpar_ndvi + fpar_sr) / 2.0, 0.001, 0.95) return fpar这段代码的关键在于同时计算两条路径后取平均，通过 np.clip 限制输出范围，避免极值外推导致 FPAR 大于 0.95 或小于 0。分母里的 1e-10 是为了防止 NDVI 恰好等于 1 时除零，实际数据虽然很少出现，但加上更保险。注意 NDVI 输入必须是 -1 到 1 的浮点数，如果原始数据是 0-200 缩放后的 DN 值，需要先除以 200，后面还会讲到这个坑。
温度胁迫系数 Ts 的计算，不少刚接触 CASA 的人会理解为"温度越高越好"，这是理解偏差。光能利用率对温度不是单调响应，而是存在一个最适温度窗口。
2.2 光能利用率 LUE：温度与水分胁迫的正确折减逻辑
LUE 的完整表达是：
LUE = LUE_max × Ts × Ws其中 LUE_max 是植被在理想条件下的最大光能利用率，单位是 g C / MJ。脚本里默认取值 0.542，这个值对应常绿阔叶林的常用经验值；如果做的是草地或农田，建议往下调到 0.3-0.4 之间，否则结果会明显偏大。
温度胁迫系数 Ts 由两部分构成。第一部分反映区域最适生长温度本身的差异：
Ts1 = 0.8 + 0.02 × T_opt - 0.0005 × T_opt²T_opt 是每个像元一年中 NDVI 峰值月份所对应的当月平均气温。这个公式的含义是：不同气候区的植被有各自不同的最适温度，把 T_opt 代进去会得到一条类似抛物线的曲线，T_opt 在 15-25℃ 时 Ts1 接近 1，过低或过高的最适温度都会让基准效率下降。
第二部分反映当月实际平均气温偏离最适温度的程度：
dT = |T_ave - T_opt| 当 dT <= 10 时，Ts2 = 1 当 dT >= 20 时，Ts2 = 0 中间线性插值这个设计的逻辑是：即使一个地区最适温度是 20℃，如果当月均温只有 5℃，光合作用相关的酶活性会大幅下降，光能利用率相应折减。两部分相乘才是完整的 Ts。脚本里为了方便，把 Ts1、Ts2 整合在一个函数里，但调试时建议分别输出，否则你很难判断结果偏低是来自"最适温度本身不合适"还是"当月偏离过大"。
水分胁迫系数 Ws 同样不可忽略。脚本默认通过土壤有效含水量 SW 来估算：
Ws = 0.5 + 0.5 × (SW - SW_wilt) / (SW_field - SW_wilt)SW_wilt 是萎蔫系数，SW_field 是田间持水量，当土壤含水量低于萎蔫点时 Ws 固定为 0.5，高于田间持水量时 Ws 为 1。这是 CASA 原版文献里的经典做法。但如果你手里没有土壤湿度栅格，脚本也允许跳过 Ws，默认 1.0——这等于假设了水分充足。这个假设在湿润区问题不大，在干旱半干旱区会让 NPP 高估 30% 以上。
2.3 月尺度合成：把 APAR 落到逐月栅格并生成年度累计
有了 FPAR 和 LUE，就可以计算月 NPP：
APAR = FPAR × 0.5 × SOL NPP_month = APAR × LUE其中 SOL 是当月太阳总辐射，单位是 MJ / m² / 月。乘以 0.5 是因为光合有效辐射大约占太阳总辐射的一半，这是一个广泛使用的近似。NPP_month 的单位是 g C / m² / 月，把 12 个月累加才能得到年 NPP，单位 g C / m² / 年。
这里要特别确认 SOL 的原始单位。如果数据源给的是 W/m² 的月均通量密度，要先乘以当月天数再乘以 86400 秒换算成 J/m²，再除以 1e6 得到 MJ/m²。否则算出来的 NPP 可能差 3 个数量级。脚本主流程里专门做了单位归一化处理，见后面的代码实现。
另外需要注意，CASA 原版是按月步长计算的，所以输入的 NDVI、气温、辐射都应该是月尺度栅格，而不是日值。如果你手里只有生长季几个月的数据，脚本也支持只跑指定月份，但输出文件命名逻辑要自己对应好。
3. 示例数据与预处理：跑脚本前先理顺的四个细节
这套资源的核心价值在于自带一份可直接运行的示例遥感影像数据。但任何遥感数据到了手上，第一步永远不是计算，而是检查数据格式、投影、分辨率、无效值这四个基础问题。预处理没做干净，后面所有中间量的空间一致性都会出问题。
3.1 数据文件结构与命名约定
示例数据按类别放在不同子目录下，命名规则统一为变量_年份_月份.tif：
example_data/ ├── ndvi/ │ ├── ndvi_2019_01.tif │ └── ndvi_2019_02.tif ├── temp/ │ ├── temp_2019_01.tif │ └── temp_2019_02.tif ├── rad/ │ ├── rad_2019_01.tif │ └── rad_2019_02.tif └── soil/ ├── soil_2019_01.tif └── soil_2019_02.tif之所以按变量分目录而不是按月份分目录，是因为脚本主循环需要先按月份组织文件列表，脚本默认会扫描指定前缀和月份的 tif 文件，如果缺某个月的某个变量会报错并给出明确提示，而不是静默跳过。
示例数据做成了小范围区域，主要有两个原因：一是让首次运行能在几秒内完成全流程验证，二是方便你直接拿输出和脚本自带的预期值对比。正式开展研究时，只需要把 input_dir 路径换成你自己的数据目录、保持同样的命名规则即可，核心计算代码不需要改动。
3.2 投影与分辨率对齐：重采样是第一步
这是最常见的数据预处理问题。NDVI 产品可能是正弦投影或 WGS84 经纬度网格，气温和太阳辐射可能来自不同部门和再分析数据集，分辨率从 250m 到 1km 不等。如果不对齐，直接读成数组做运算，结果就是影像错位、边界黑边。
我一般会先用 gdalwarp 把所有输入统一到一个基准网格上：
# 以 NDVI 数据为基准，重投影并重采样气温数据到同一网格 gdalwarp -t_srs EPSG:32650 -tr 1000 1000 -r bilinear \ -overwrite input/temp/temp_2019_01.tif \ output/temp_aligned/temp_2019_01.tif # 同样处理太阳辐射 gdalwarp -t_srs EPSG:32650 -tr 1000 1000 -r cubicspline \ -overwrite input/rad/rad_2019_01.tif \ output/rad_aligned/rad_2019_01.tif参数说明：-t_srs 指定目标投影（示例用 UTM 50N），-tr 设置目标像元大小 1000×1000 米，-r bilinear 是双线性插值，适合气温这类连续场；NDVI 和 FPAR 这类类指数变量重采样时更推荐 cubicspline 或 nearest，避免插值产生超界值。需要注意，重采样后的像元网格必须与 NDVI 完全一致，光调整分辨率还不够，贴一个以 NDVI 掩膜为模板的裁剪步骤会更稳妥。
脚本里还内置了 rasterio 方式的重投影函数，直接在 Python 内完成，适合批量处理多个月份数据：
import rasterio from rasterio.warp import reproject, Resampling def align_to_reference(src_path, ref_path, dst_path): with rasterio.open(ref_path) as ref: transform, width, height = ref.transform, ref.width, ref.height crs = ref.crs with rasterio.open(src_path) as src: data = src.read(1) profile = src.profile.copy() profile.update( driver='GTiff', crs=crs, transform=transform, width=width, height=height ) dst = np.empty((height, width), dtype=np.float32) reproject( source=data, src_transform=src.transform, src_crs=src.crs, destination=dst, dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear ) with rasterio.open(dst_path, 'w', **profile) as dst_rst: dst_rst.write(dst, 1)这段代码的核心逻辑是以参考文件（NDVI）的地理范围、行列数和投影为基准，把其他数据投影到相同网格。reproject 是 rasterio 的重投影接口，selection 参数里有 nearest、bilinear、cubicspline 三种常见模式，温度辐射用 bilinear，土地利用类型和掩膜用 nearest。注意处理浮点数据的 nodata，重采样后边缘像元可能产生无效值，需要再检查一次。
3.3 无效值与水体掩膜：NDVI 有效范围裁剪
遥感的原始产品经常把云、冰雪、异常值编码成 -9999、-32768 或 0。这些值如果不处理，会直接污染 FPAR 计算里的 NDVI 极值统计。脚本在读取数据后统一做以下操作：
def mask_valid(arr, nodata=-9999, valid_range=None): mask = (arr != nodata) & np.isfinite(arr) if valid_range is not None: mask &= (arr >= valid_range[0]) & (arr <= valid_range[1]) return np.where(mask, arr, np.nan)对 NDVI，valid_range 建议设成 [-0.1, 1.0]，因为健康的植被 NDVI 不会小于 0，但土壤背景或稀疏植被可能出现负值，直接排除会损失边缘信息，所以留出 -0.1 作为缓冲；对气温，valid_range 根据区域气候带设成 [−50, 50] 就足够；太阳辐射则必须大于 0。所有变量统一掩膜后再写回数组，进入正式计算流程。
顺带提一下水体问题。CASA 模型默认所有像元都有植被覆盖，如果研究区有湖泊、水库或城市建成区，这些像元算出来的 NPP 没有生态学意义，会在统计均值时拉低结果。脚本提供了两个处理方案：一是用 mask_built.tif 这类外部掩膜直接把非植被像元置 NaN；二是简单地把 NDVI 低于 0.05 的像元在计算时跳过。方案二虽然粗糙，在裸地和稀疏草原区域会误伤低覆盖度植被，但对大多数研究场景已经够用。
4. Python 计算脚本拆解：把公式落成逐像元可复现的代码
这一章是整个资源的核心部分。直接贴全部代码没有意义，我会把工程组织、关键函数和主流程三个层面逐段拆开，并说明每个参数怎么改、改完影响什么。
4.1 工程结构与配置入口
脚本工程结构如下：
casa_npp/ ├── config.yaml ├── calc_casa.py ├── preprocess.py ├── utils.py └── README.mdconfig.yaml 集中管理所有经验参数，改参数不需要动代码：
input: ndvi_dir: ./example_data/ndvi temp_dir: ./example_data/temp rad_dir: ./example_data/rad soil_dir: ./example_data/soil # 可选，没有就留空 output: out_dir: ./output/npp prefix: npp year: 2019 months: [1,2,3,4,5,6,7,8,9,10,11,12] model: lue_max: 0.542 ndvi_min: 0.05 ndvi_max: 0.95 sr_min: 1.08 sr_max: 5.0 fpar_min: 0.001 fpar_max: 0.95 water_stress: true # 没有土壤湿度数据时设为 false每个参数的生态学含义：lue_max 是最大光能利用率，植被类型差异很大，草地建议 0.3-0.4，农田 0.4-0.5，常绿阔叶林 0.5-0.7，不要一套参数打天下；ndvi_min/max 是 FPAR 拉伸的参考极值，建议用研究区历史 NDVI 分位数动态估算，而不是死用默认值；water_stress 控制是否启用水分胁迫折减，缺土壤数据时置为 false，但要意识到结果是干区高估的。
工程上这样的设计避免了在代码深处找参数字面量，也方便批量跑不同年份和不同植被类型。我自己的使用习惯是：每种土地利用类型单独建一个 config，循环跑完后汇总。
4.2 FPAR 与胁迫系数：numpy 向量化的核心函数
计算核心全部基于 numpy 数组的向量化操作，没有逐像元循环。这里贴出 LUE 相关的实现：
def calc_lue(lue_max, t_ave, t_opt, sw=None, sw_wilt=0.1, sw_field=0.4): # 温度胁迫第一项：最适温度基准效率 ts1 = 0.8 + 0.02 * t_opt - 0.0005 * t_opt ** 2 ts1 = np.clip(ts1, 0.0, 1.0) # 温度胁迫第二项：当月均温与最适温度的偏离 d_temp = np.abs(t_ave - t_opt) ts2 = np.ones_like(d_temp, dtype=np.float32) ts2[d_temp > 20.0] = 0.0 between = (d_temp > 10.0) & (d_temp <= 20.0) ts2[between] = 1.0 - (d_temp[between] - 10.0) / 10.0 ts = ts1 * ts2 # 水分胁迫 if sw is not None: ws = 0.5 + 0.5 * (sw - sw_wilt) / (sw_field - sw_wilt) ws = np.clip(ws, 0.5, 1.0) else: ws = 1.0 return lue_max * ts * ws这段代码的逻辑：ts1 是抛物线型函数，最适温度在 20℃ 附近时效率最高；ts2 用分段线性折减表达偏离惩罚，偏离超过 20℃ 直接归零。水分胁迫 ws 通过土壤有效含水量的线性插值实现，区间限制在 0.5-1.0，这是因为植被在严重干旱时也不可能光合作用完全停止，保留一半的基础效率是 CASA 原版的设定。t_opt 需要在主流程里先计算，方法是：对每个像元找出一年中 NDVI 最大的月份，再把对应月份的 t_ave 取出来作为该像元的 t_opt。
4.3 主流程：逐月读取、逐像元计算、年度累计输出
主流程伪代码和骨架如下：
import yaml import rasterio import numpy as np from utils import mask_valid, calc_fpar, calc_lue with open('config.yaml') as f: cfg = yaml.safe_load(f) t_opt = compute_t_opt(cfg['input']['ndvi_dir'], cfg['input']['temp_dir'], cfg['year']) annual_npp = None for month in cfg['months']: ndvi = read_aligned_month('ndvi', month) temp = read_aligned_month('temp', month) rad = read_aligned_month('rad', month) sw = read_aligned_month('soil', month) if cfg['water_stress'] else None ndvi = mask_valid(ndvi, nodata=-9999, valid_range=[-0.1, 1.0]) temp = mask_valid(temp, nodata=-9999, valid_range=[-50, 50]) fpar = calc_fpar(ndvi, cfg['model']['ndvi_min'], cfg['model']['ndvi_max']) lue = calc_lue(cfg['model']['lue_max'], temp, t_opt, sw) apar = fpar * 0.5 * rad npp_m = apar * lue # gC/m2/month if annual_npp is None: annual_npp = npp_m.copy() else: annual_npp = np.nansum(np.stack([annual_npp, npp_m]), axis=0) write_geotiff(cfg['output']['out_dir'], f"npp_{cfg['year']}_{month:02d}.tif", npp_m) write_geotiff(cfg['output']['out_dir'], f"npp_{cfg['year']}_annual.tif", annual_npp)主流程的执行顺序是：先计算全年的 t_opt 基准场，再按月逐次读取数据、计算 FPAR、LUE、NPP，最后把 12 个月结果累加写出。两个容易看漏的细节：annual_npp 的累加用的是 nansum，这样某个月份某些像元因为无效值产生 NaN 时，不会被污染到全年结果里；write_geotiff 的 profile 需要复用输入数据的 transform、crs，否则输出栅格的地理参考会丢失。内存方面，按月读取再累加，避免一次性把 12 个月的数组全部堆进内存，这是脚本能跑大面积区域的原因。
5. 计算避坑与常见问题排查：现象、原因和后悔药
这部分是实际跑这个脚本最容易翻车的点，我把它们按"现象→原因→解决"的格式整理出来。每一条都是真实出现过的，不是穷举猜测。
5.1 输出全为 0 或整片常数：无效值没清干净
现象：跑完以后 NPP 栅格几乎全是 0，或者全是一个相同的低值，但输入 NDVI 明明有空间变化。
原因：最常见的是 nodata 值 -9999 没有进入 mask_valid 流程，直接参与 FPAR 计算。比如 -9999 被当作有效 NDVI 参与极值统计后，会把 ndvi_max 拉大到几百，正常像元的拉伸区间全部被压缩到接近 0。另一种情况是气温数据没做单位换算，原始数据是 0.1℃ 缩放的整数，直接当摄氏度用后，全年温度都低于 -10℃，温度胁迫系数 Ts 被归零，LUE 成为 0。
解决：读入每个变量后先用 mask_valid 处理，再打印一个像元的中间量（ndvi、fpar、ts、ws、npp）核对量级。脚本里预留了 debug 模式，会输出首行首列像元的过程值。养成先看单像元中间量、再看全图均值的习惯，能避免在最后阶段才发现问题。
5.2 影像错位或边缘黑边：投影和分辨率不一致
现象：输出的 NPP 栅格和原始 NDVI 影像在地理上对不上，比如水体区域在高值区，山区边界出现锯齿和黑边。
原因：不同来源数据的投影或像元网格没有对齐。NDVI 产品可能是 UTM，气温和辐射是经纬度网格，直接重新采样后数组运算，坐标系不同的数组按行列对齐，必然错位。黑边则是重采样后覆盖范围小于目标网格，边缘被 nodata 填充，参与了统计导致结果异常。
解决：在预处理阶段强制把每个月的输入统一重采样到 NDVI 的网格上。先对比几个控制点坐标是否吻合，再用 gdalinfo 检查投影、像元大小、原点坐标，三个参数完全一致才允许进入计算。这一步没有自动化捷径，一次性做干净，后面整年的计算都不用担心空间错位。
5.3 内存爆掉：一次读入全年数据
现象：跑小范围示例数据很快，换成全省或全国范围后，运行到第 5-6 个月直接 OutOfMemory 崩溃。
原因：代码在读取时一次性把 12 个月的 NDVI、气温、辐射全部读入内存列表，再开始循环计算，50 万平方公里、1km 分辨率就是几十万行乘几十万列的浮点数组，每层几百 MB，12 个月累积就是 10GB 以上。
解决：主流程改成按月读取、逐月计算、边算边写。年度累计只保留一个 annual_npp 数组，不保留每个月的原始数据。如果单月数据仍然过大，再用 rasterio.Window 分块处理，每读一块就计算一块。脚本的主流程已经支持按月循环，内存占用基本等于单月四层栅格加一个年度累计层。
5.4 结果量级不合理：单位换算是头号凶手
现象：年 NPP 算出来上万 g C/m²，远超植被碳固定的物理上限；或者反过来，结果小得只有个位数。
原因：太阳辐射单位没换算。如果原始辐射是 J/m²，而公式里用的是 MJ/m²，结果会大 1e6 倍；如果气温原始数据是 0.1℃ 缩放的短整型没除 10，温度胁迫系数 Ts 会整体偏低。NDVI 也一样，0-200 缩放的 DN 值没除以 200，FPAR 会被大幅低估。
解决：在 config.yaml 里增加 unit_params 小节，显式声明每个变量的原始单位和缩放系数，脚本读取后统一转成内部标准单位。例如 rad_scale=0.001 表示原始单位是 J/m²、除以 1000 转成 MJ/m²，temp_scale=0.1 表示原始单位为 0.1℃。不确定时先对单像元坐标做数值核对，查一下该地区该月的辐射和气温实测值是否在同一量级。
6. 落地技巧：用统计均值和年际对比给计算结果做体检
NPP 计算跑出结果只是第一步，结果可信不可信需要有验证手段。这里分享两个我每次跑完都会做的检查，以及对应的小技巧。
6.1 统计均值和已知参考区间做量级核对
把输出的年度 NPP 栅格做基本的 min/mean/max 统计，再对照区域生态系统的常见量级区间。阔叶林年 NPP 一般在 800-1500 g C/m² 附近，草地多在 200-600，农田正常波动在 300-700。如果统计均值明显超出对应生态系统的区间，优先回头查单位和 LUE_max，而不是怀疑模型本身。
可以快速跑一个统计验证：
import rasterio import numpy as np with rasterio.open('output/npp/npp_2019_annual.tif') as src: annual = src.read(1) valid = annual[(annual > 0) & (annual < 3000)] print('mean:', valid.mean(), 'median:', np.median(valid))这一步的意义在于快速暴露量级错误。吃过一次亏以后，我把这步做成固定流程，任何新区域数据到手，先跑通一年的统计再谈进一步分析。统计结果偏差超过 30% 时，不要急着调参数，先检查数据预处理和单位换算。
6.2 年际对比：提取像元改变方向，而不只是数值差
NPP 年际变化分析里比较实用的一个技巧是计算两期差值后，叠加变化方向掩膜，把显著增加、显著减少和相对稳定的像元分开统计。简单实现方式：
delta = npp_2020 - npp_2019 increase = (delta > threshold).astype(np.int16) decrease = (delta < -threshold).astype(np.int16) stable = 1 - increase - decreasethreshold 建议取本次计算标准差的一半，也可以取区域均值的 5%-10%，具体看研究目的。这样处理后的三个类别栅格，后续做变化检测和空间制图都比直接用连续差值直观。
从那以后我每次跑 NPP 都强制走一遍这个流程：先小范围试算、核对单像元中间量、统计均值对照文献区间，最后才放全量。这套脚本最让我省心的地方就是中间过程都能输出检查，不把模型当黑匣子。希望帮到你。
本文还有配套的精品资源，点击获取
