用Python实现CASA模型NPP计算：从NDVI到光能利用率的完整流程

简介这是一套基于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 算 NPPCASA 模型从公式到可跑脚本的完整闭环遥感生态方向拿到 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 到植被吸收的有效光合辐射FPARFraction of Photosynthetically Active Radiation表示植被冠层吸收光合有效辐射的比例取值范围一般在 0 到 0.95 之间。最常用的估算路径是归一化植被指数 NDVI原理是 NDVI 对植被冠层叶面积和吸收效率敏感可以直接用极值线性拉伸映射到 FPARFPAR_NDVI (NDVI - NDVI_min) / (NDVI_max - NDVI_min) × (FPAR_max - FPAR_min) FPAR_min不过在实践中单靠 NDVI 一条路径容易出现饱和问题尤其是茂密植被区 NDVI 接近 0.9 以后增长乏力FPAR 会被低估。所以很多实现会引入简单比值 SRSimple 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_min0.05、NDVI_max0.95、SR_min1.08、SR_max5.0、FPAR_min0.001、FPAR_max0.95这些是文献里最常见的经验取值。这里有一个容易犯的错跑到茂密森林区域时如果发现 FPAR 常年顶到 0.94 以上不要急着怀疑公式先检查 NDVI_max 是不是被个别异常像元比如云边或水体的高值噪声拉大了导致正常像元的拉伸区间被压缩。NDVI 极值建议用数据自身分位数去替代全局 min/max脚本在配置里预留了这个开关。用 numpy 可以很直接地实现这个双路径计算向量化处理比逐像元 for 循环快几个量级import numpy as np def calc_fpar(ndvi, ndvi_min0.05, ndvi_max0.95, sr_min1.08, sr_max5.0, fpar_min0.001, fpar_max0.95): ndvi np.clip(ndvi, 0.0, 1.0) # 路径 1NDVI 线性拉伸 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就可以计算月 NPPAPAR 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 数据文件结构与命名约定示例数据按类别放在不同子目录下命名规则统一为变量_年份_月份.tifexample_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( driverGTiff, crscrs, transformtransform, widthwidth, heightheight ) dst np.empty((height, width), dtypenp.float32) reproject( sourcedata, src_transformsrc.transform, src_crssrc.crs, destinationdst, dst_transformtransform, dst_crscrs, resamplingResampling.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_rangeNone): 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)对 NDVIvalid_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, swNone, sw_wilt0.1, sw_field0.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, dtypenp.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]), axis0) write_geotiff(cfg[output][out_dir], fnpp_{cfg[year]}_{month:02d}.tif, npp_m) write_geotiff(cfg[output][out_dir], fnpp_{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 分辨率就是几十万行乘几十万列的浮点数组每层几百 MB12 个月累积就是 10GB 以上。解决主流程改成按月读取、逐月计算、边算边写。年度累计只保留一个 annual_npp 数组不保留每个月的原始数据。如果单月数据仍然过大再用 rasterio.Window 分块处理每读一块就计算一块。脚本的主流程已经支持按月循环内存占用基本等于单月四层栅格加一个年度累计层。5.4 结果量级不合理单位换算是头号凶手现象年 NPP 算出来上万 g C/m²远超植被碳固定的物理上限或者反过来结果小得只有个位数。原因太阳辐射单位没换算。如果原始辐射是 J/m²而公式里用的是 MJ/m²结果会大 1e6 倍如果气温原始数据是 0.1℃ 缩放的短整型没除 10温度胁迫系数 Ts 会整体偏低。NDVI 也一样0-200 缩放的 DN 值没除以 200FPAR 会被大幅低估。解决在 config.yaml 里增加 unit_params 小节显式声明每个变量的原始单位和缩放系数脚本读取后统一转成内部标准单位。例如 rad_scale0.001 表示原始单位是 J/m²、除以 1000 转成 MJ/m²temp_scale0.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 都强制走一遍这个流程先小范围试算、核对单像元中间量、统计均值对照文献区间最后才放全量。这套脚本最让我省心的地方就是中间过程都能输出检查不把模型当黑匣子。希望帮到你。本文还有配套的精品资源点击获取