简介本资源是面向生态建模与遥感研究者的CASACarbon Assimilation by Sun and Shade leaves in Annual plants模型Python实现方案聚焦于净初级生产力NPP的自动化计算与分析适用于气候变化影响评估、植被碳汇模拟及农业生态优化等科研场景。压缩包共2个文件1个.tif遥感栅格数据用于驱动模型1个.py核心脚本实现完整CASA算法流程总大小20.05MB结构精简、开箱即用。已有2396人学习下载体现其在环境科学与Python交叉领域的实用热度。用户可直接运行脚本完成NPP时空模拟代码内嵌气象因子光照、温度、CO₂、植被参数叶面积指数与阴阳叶光合区分逻辑并支持基于.tif输入的批量计算与结果可视化基础框架为后续模型调参、区域验证或教学演示提供可靠起点。 CASA模型的Python实现这个话题我拖了快一年才整理成文。之前帮课题组做过几版NPP估算的脚本也踩了不少数据处理和参数标定的坑这次把完整的思路、公式拆解、代码结构和调试经验一次讲清楚。有朋友可能在搜索时看到“cass建模”这个词先说明一下这里的CASA是Carnegie-Ames-Stanford Approach的缩写是生态遥感领域用来估算植被净初级生产力NPP的光能利用率模型跟测绘里那套CASS成图系统完全是两回事。这篇文章适合谁看第一类是刚入门遥感生态建模的研究生想把CASA模型从公式变成能跑的代码第二类是有一定Python基础、但不太清楚生态模型怎么和遥感数据打交道的朋友第三类是已经跑通过模型、想优化参数标定和验证流程的老手。内容按“原理拆解→数据准备→代码实现→参数标定→验证评估→问题排查”这个路径来组织跟着走基本能把月尺度的区域NPP估算流程搭起来。1. CASA模型的核心逻辑从光能利用率的视角看植被生产力1.1 公式拆解NPP APAR × ε 这个等式的含义CASA模型之所以长盛不衰核心在于它用一个非常简洁的线性关系描述了植被光合作用积累生物量的过程NPP APAR × ε其中APAR是植被吸收的光合有效辐射单位是MJ/m²ε是光能利用率单位是g C/MJ。两个量一乘就得到单位面积上植物通过光合作用固定的碳也就是净初级生产力NPP单位是g C/m²。这里有个关键点需要理解APAR描述的是“植物能用来做光合作用的能量有多少”ε描述的是“植物把这些能量转化成碳的效率有多高”。用生活化的比喻来说APAR相当于你手里的食材总量ε就是你的厨艺水平NPP就是最终端上桌的菜量。食材再多厨艺不行做出来的菜也有限厨艺再好没食材也不行。分两段来说APAR由太阳总辐射、植被覆盖度和光合有效辐射比例共同决定计算公式是 APAR SOL × 0.5 × FPAR。其中0.5是光合有效辐射400-700nm占太阳总辐射的比例FPAR是植被吸收的光合有效辐射比例反映植被冠层对光的拦截能力。ε在CASA模型里不是常数它受温度和水分两个胁迫因子调节ε ε_max × T_ε × W_ε。ε_max代表理想条件下植被的最大光能利用率T_ε是温度胁迫系数W_ε是水分胁迫系数。温度太低、太高或者水分不足都会让实际光能利用率打折扣。所以说到底CASA模型的Python实现核心就是写清楚这三个系数的计算过程再把它们串起来乘到NPP上。代码难度不大真正花时间的是把输入数据准备好、把参数标定对。1.2 为什么选择Python而不是ENVI/IDL传统方案早些年做CASA模型很多人是用ENVI的IDL语言或者ArcGIS的Model Builder来跑的因为生态遥感的老前辈们留下了大量IDL脚本。但我个人强烈建议用Python来做原因有三点。第一数据处理的生态更完整。CASA模型需要NDVI、气温、降水、太阳辐射四类数据这些数据通常以GeoTIFF或NetCDF格式存储Python里的gdal、rasterio、xarray可以无缝读取、裁剪、重采样、堆叠整个流程可以写成一个自动化pipeline。IDL虽然也能处理但生态圈子太小遇到问题网上的解决方案也少。第二参数标定和敏感性分析更方便。CASA模型里ε_max、FPAR最大最小值的确定都需要大量统计计算Python的numpy和pandas处理这些正合适。后面做验证时需要把模型输出的NPP和实测数据做回归分析scipy、sklearn这些库直接套上去就行。第三结果可视化和发布更顺手。matplotlib、cartopy画出来的图无论是论文插图还是项目汇报都拿得出手。ENVI出图虽然也够用但风格偏老气定制性差。当然Python也不是万能的。如果你只有几个小流域的栅格数据用ENVI手动点几下可能更快。但只要是做多年份、大数据量、区域尺度的NPP估算Python的批处理和工程化优势就非常明显了。1.3 本实现方案的适用范围我这里讲的这套实现目标场景是基于MODIS NDVI和气象栅格数据估算月尺度或年尺度的区域NPP。空间分辨率可以从250m到1km时间跨度可以是几年到几十年区域范围从一个小流域到全国尺度都可以跑。这套代码不适用于单株植物或田块尺度的生产力估算那种场景应该用过程模型如DSSAT、APSIM或者涡度相关观测来衡量。另外如果只有站点数据没有空间连续的气象场建议先把气象数据做空间插值如ANUSPLIN、克里金再跑模型否则空间连续性会出问题。2. 建模前必须搞定的数据准备输入数据是模型的命根子2.1 四类核心输入数据清单CASA模型需要四类数据我把它们列成了一张清单方便你对照准备数据类别常用数据源时间分辨率空间分辨率主要用途NDVI数据MODIS MOD13Q1/MOD13A1、SPOT/VGT、Landsat16天/月250m/500m/1km计算FPAR气温数据ERA5-Land、CN05.1格点数据、GLDAS月/日0.1°等计算温度胁迫系数降水数据ERA5-Land、CN05.1、TRMM/GPM月/日0.1°等与蒸散一起计算水分胁迫太阳辐射数据ERA5-Land向下短波辐射、GLASS辐射产品月/日0.1°/5km计算APAR这里要提醒一点不同来源数据的投影、空间范围、分辨率往往不一致比如MODIS是正弦曲线投影ERA5-Land是等经纬度网格。拿到数据后第一步就是把所有数据统一到同一个空间参考下否则后面计算时栅格对不齐出来的结果可能是错乱的。2.2 时间尺度的选择为什么推荐月尺度CASA模型理论上可以按天、按旬、按月运行但我强烈建议第一次实现用月尺度原因有三层。其一MODIS NDVI虽然有16天合成的产品但直接拿来算年尺度或日尺度的FPAR会有很多噪音月最大值合成能有效降低云污染和大气气溶胶的影响。其二气温和降水在月尺度上相对稳定T_ε和W_ε的计算不容易出现剧烈波动导致的极端值。其三月尺度的计算量适中迭代调试速度快等到流程跑通了再加密到旬尺度也不难。实际操作中通常把一年12个月的NDVI最大合成值作为模型输入对应12个月的气温、降水、太阳辐射数据。当年NPP就是12个月NPP之和。2.3 单位换算与坐标系统一最容易翻车的环节单位换算这一步是新手最容易翻车的地方我在这里栽过跟头必须展开说。首先看辐射数据。ERA5-Land的向下短波辐射单位是J/m²累计值而CASA模型需要的是MJ/m²。如果原始数据是每小时的累计辐射要先把一天内所有时次加起来得到日总量再除以1e6转成MJ/m²最后乘以当月天数得到月总量。其次看NDVI取值范围。MODIS NDVI产品理论范围是-1到1但实际数据里水体、云、雪等像元的NDVI可能为负或接近0必须事先做掩膜处理。FPAR计算时NDVI的有效范围通常在0.1到0.8之间低于0.1的值按裸地处理高于0.8的要检查是否存在异常。再看坐标系。不同投影的栅格数据如果直接用原始分辨率相除或相乘结果会有明显的接边痕迹。我的建议是统一重采样到同一个参考网格比如把MOD13Q1的250m NDVI重采样到1km与气象数据保持一致。重采样方法上NDVI这类连续变量用双线性插值或立方卷积土地利用类型这类离散变量用最近邻。最后是气温单位。ERA5-Land的气温单位是开尔文而CASA模型里温度胁迫系数的计算用的是摄氏度两者相差273.15度。这个低级错误真的有人犯过算出来的T_ε全是0因为温度都在-270℃以下。3. 核心模块的Python实现逐段写出每一行代码3.1 数据读取与预处理用rasterio把栅格读成数组正式写核心计算前先把输入数据读取的代码框架搭好。我用rasterio来读写GeoTIFF用numpy做数组运算这套组合非常稳定。import numpy as np import rasterio def read_band(path, band1): with rasterio.open(path) as src: data src.read(band).astype(np.float32) profile src.profile transform src.transform # 将无效值设为NaN注意不同产品的无效值可能不同 data[data -9999] np.nan return data, profile, transform这里有个细节值得注意。MODIS产品的无效值通常是-3000或-9999而ERA5的陆面数据在海洋区域可能是NaN。统一在读取阶段就把无效值替换成np.nan后面所有计算遇到np.nan都会自动传递不用每个函数都检查一遍有效性。3.2 FPAR的计算两种方法的平均CASA模型里FPAR的估算有两种经典方法。第一种基于NDVI的线性关系公式是FPAR_ndvi (NDVI - NDVI_min) / (NDVI_max - NDVI_min) × (FPAR_max - FPAR_min) FPAR_min其中FPAR_min取0.001FPAR_max取0.95NDVI_min和NDVI_max对应植被类型NDVI的5%和95%分位数。第二种基于植被指数SRSimple Ratio公式是SR (1 NDVI) / (1 - NDVI) FPAR_sr (SR - SR_min) / (SR_max - SR_min) × (FPAR_max - FPAR_min) FPAR_min实际研究中通常取这两种结果的算术平均值作为最终的FPAR目的是减少单一方法的系统性偏差。代码实现就按这个逻辑来写def calc_fpar(ndvi, ndvi_min, ndvi_max): fpar_min 0.001 fpar_max 0.950 # 方法一基于NDVI的线性关系 ndvi_clip np.clip(ndvi, ndvi_min, ndvi_max) fpar_ndvi (ndvi_clip - ndvi_min) / (ndvi_max - ndvi_min) fpar_ndvi fpar_ndvi * (fpar_max - fpar_min) fpar_min # 方法二基于SR指数 sr (1 ndvi_clip) / (1 - ndvi_clip 1e-8) sr_min (1 ndvi_min) / (1 - ndvi_min 1e-8) sr_max (1 ndvi_max) / (1 - ndvi_max 1e-8) fpar_sr (sr - sr_min) / (sr_max - sr_min 1e-8) fpar_sr np.clip(fpar_sr, 0, 1) * (fpar_max - fpar_min) fpar_min # 取两者的平均值 fpar (fpar_ndvi fpar_sr) / 2.0 fpar np.clip(fpar, 0, 1) return fpar代码里加1e-8是防止NDVI等于1时除数为0的极端情况。这个方法属于典型的数值稳定性处理实测中确实遇到过MODIS个别像元NDVI接近1的情况。3.3 温度胁迫系数Tε的计算公式温度对光能利用率的影响分成两部分。第一部分反映植物在最适温度下的光合效率上限公式是T_ε1 0.8 0.02 × T_opt - 0.0005 × T_opt²第二部分反映实际月均温和最适温度的偏离程度公式是T_ε2 1 / [1 exp(0.2 × (T_opt - 10 - T_month))]最终的T_ε T_ε1 × T_ε2。注意T_opt是研究区植被生长的最适温度通常取全年NDVI最大值对应月份的月均温一种简化做法是用生长季内NDVI≥0.5时对应月份的多年平均气温。def calc_temperature_stress(t_month, t_opt): # 第一部分最适温度下的效率上限 t_eps1 0.8 0.02 * t_opt - 0.0005 * t_opt**2 t_eps1 np.clip(t_eps1, 0, 1) # 第二部分温度偏离惩罚 t_eps2 1.0 / (1.0 np.exp(0.2 * (t_opt - 10 - t_month))) t_eps t_eps1 * t_eps2 t_eps np.clip(t_eps, 0, 1) return t_eps这里有个隐含假设T_opt在模拟期内是固定的。但如果研究区跨度大地形差异明显不同像元的T_opt应该不同。更严谨的做法是用12个月NDVI找每个月最大NDVI对应的温度逐像元求T_opt。3.4 水分胁迫系数Wε蒸散比的运用CASA模型中的水分胁迫系数用实际蒸散与潜在蒸散的比值来表征W_ε 0.5 0.5 × EET / PETEET是实际蒸散量PET是潜在蒸散量。当实际蒸散等于潜在蒸散时W_ε1说明水分供应充足光能利用率不受水分限制。当水分严重不足时EET接近0W_ε降到0.5光能利用率减半。实际蒸散和潜在蒸散的获取有两个途径。一是直接用再分析产品提供的蒸散变量比如ERA5有总蒸散量但需要注意单位换算和存储类型。二是用Thornthwaite公式或Penman-Monteith公式自己算潜在蒸散再用土壤水分平衡模型估算实际蒸散。第一种简单快速第二种精度更高但也更复杂。我的建议是第一次实现直接用ERA5或GLDAS的蒸散产品后续再考虑用物理公式精化。因为水分胁迫系数本身在湿润地区的影响不是特别大精化这部分投入产出比不高。def calc_water_stress(eet, pet): pet_safe np.where(pet 0, pet, np.nan) ratio eet / pet_safe ratio np.clip(ratio, 0, 1) w_eps 0.5 0.5 * ratio w_eps np.clip(w_eps, 0, 1) return w_eps3.5 ε计算与NPP汇总有了温度胁迫系数和水分胁迫系数光能利用率就很直接了。ε_max是给定植被类型的最大光能利用率后面单独用一节来讨论如何取值。这里先按Potter等1993年提出的全球平均默认值0.389 g C/MJ来写后续再替换成经过标定的值。def calc_epsilon(epsilon_max, t_eps, w_eps): return epsilon_max * t_eps * w_eps def calc_npp_annual(apar_list, fpar_list, t_eps_list, w_eps_list, epsilon_max): npp_months [] for i in range(12): apar apar_list[i] # MJ/m2/month fpar fpar_list[i] t_eps t_eps_list[i] w_eps w_eps_list[i] apar_absorbed apar * fpar eps calc_epsilon(epsilon_max, t_eps, w_eps) npp_month apar_absorbed * eps # g C/m2/month npp_months.append(npp_month) npp_annual np.sum(np.stack(npp_months, axis0), axis0) return npp_annual注意APAR的计算公式里如果输入的太阳辐射是总辐射SOL那么吸收的光合有效辐射部分应该是SOL × 0.5 × FPAR。我在代码里把0.5这层折算提前到了数据准备阶段把输入的辐射数据就当作光合有效辐射这样核心模块里直接乘FPAR就行避免每处都重复写0.5。3.6 主流程串联把模块串成流水线上面几个函数是离散的要跑通完整流程还需要一个主函数来组织数据。我习惯把一年的12个月数据放到一个列表里逐月计算后汇总。def run_casa(ndvi_monthly, tmonth_monthly, eet_monthly, pet_monthly, sol_monthly, ndvi_min, ndvi_max, t_opt, epsilon_max): npp_months [] for i in range(12): ndvi ndvi_monthly[i] tmm tmonth_monthly[i] eet eet_monthly[i] pet pet_monthly[i] sol sol_monthly[i] fpar calc_fpar(ndvi, ndvi_min, ndvi_max) t_eps calc_temperature_stress(tmm, t_opt) w_eps calc_water_stress(eet, pet) eps calc_epsilon(epsilon_max, t_eps, w_eps) apar sol * fpar # sol已经是光合有效辐射 npp_month apar * eps npp_months.append(npp_month) npp_annual np.sum(np.stack(npp_months, axis0), axis0) return npp_annual这个主流程看着简单但实际项目里每一行都要检查数据维度是否一致、异常值是否被处理。我在跑全国尺度的模拟时经常因为某个月份的降雨数据有一个像元是负值导致该像元的NPP出现极端负值最后排查了很久才发现是数据质量问题而非模型bug。4. ε_max的标定模型精度的命门4.1 默认值从哪来怎么理解它Potter等人在1993年首次提出CASA模型时给出的全球平均ε_max是0.389 g C/MJ。这个值是通过全球实测NPP数据反推出来的代表的是一个宏观均值。但你做中国东北的农田或者做亚马逊雨林直接套0.389肯定不是最优解。不同植被类型的ε_max差异很大大致规律是农作物和草地偏高常绿阔叶林居中针叶林和荒漠植被偏低。做区域研究时至少要按照IPCC土地利用分类或MODIS IGBP分类给不同植被类型赋不同的ε_max值。简单做法是参考前人在类似区域的研究文献推荐以下取值区间植被类型ε_max取值范围g C/MJ参考说明农田0.4 - 0.6作物光合效率高灌溉条件下更高草地0.35 - 0.5温带草地偏中间值常绿阔叶林0.35 - 0.45热带雨林接近0.4落叶阔叶林0.3 - 0.4生长季集中效率较高针叶林0.25 - 0.35受低温限制效率偏低荒漠灌丛0.15 - 0.25水分限制强烈4.2 用实测数据反演ε_max的完整流程如果研究区域有实测NPP数据比如样地生物量调查、涡度通量塔观测或者森林清查数据建议用实测数据反演ε_max流程分三步。第一步准备训练数据。把样地点位的坐标提取出来对应到NDVI、气温、降水、太阳辐射栅格上计算出每个样点的APAR和胁迫系数。这里要注意样点坐标和栅格像元的对应关系必须精确建议先做缓冲区分析取3×3像元的平均值降低空间配准误差。第二步构建反演方程。整理已有数据集后NPP实测值、APAR、T_ε、W_ε都是已知的未知的只有ε_max。两边相除得到每个样点的表观ε值然后统计这些表观ε值的分布。通常取中位数或生长季均值作为该区域的ε_max。# 表观ε反演示意代码 apparent_epsilon npp_observed / (apar * t_eps * w_eps) epsilon_max_calibrated np.nanmedian(apparent_epsilon)第三步验证标定效果。把反演得到的ε_max代回去重新计算NPP和留出的验证样点对比。如果相关系数提升不明显检查是不是APAR计算本身有系统偏差比如辐射数据的空间尺度不匹配。4.3 标定中常见的三个错误第一直接用稀疏的野外样点去标定大区域的ε_max样本量不够导致过拟合。建议至少要有30个以上的有效样点覆盖不同植被类型和气候条件。第二忽略了ε_max的季节变化。有些文献认为ε_max在生长季初期和末期可能偏低如果逐月反演发现明显季节趋势可以考虑按月份或按物候期设定不同ε_max。第三把ε_max标定和FPAR计算混在一起。标定ε_max时最好固定FPAR的计算方法和参数否则标定结果包含了FPAR误差模型整体的偏差来源会变得模糊。5. 验证与精度评估别让NPP结果看起来合理就交差5.1 与实测数据对比的实操要点模型估算完NPP第一件事不是画图而是和实测数据对比。实测NPP数据源包括森林普查生物量增量、通量塔的GPP估算NEPNEE分解、以及文献里发表的样地NPP测定值。对比时最常犯的错误是空间尺度不匹配。野外样方通常只有几十米见方而遥感像元是1km甚至更粗植被在空间上又不是均匀分布的两者的NPP差异可能很大。我的经验是用3×3像元的平均值去匹配样点同时把样点所在的坡向、坡度和土地利用类型信息纳入筛选剔除那些落在水体、建设用地或边界混杂像元上的样点。时间尺度上如果实测NPP是年总量模型输出也是年NPP可以直接比较。如果实测数据只有生长季的增量那必须把模型输出的非生长季NPP视为0或忽略不能把12个月都加起来再和生长季实测比。5.2 与遥感NPP产品的交叉验证没有实测数据时退而求其次可以做产品交叉验证。主流的全球NPP产品有MODIS MOD17A3500m年尺度和GLASS NPP产品。把CASA模型结果重采样到和产品相同的网格上然后做逐像元的散点图和相关分析。交叉验证要注意一个问题MOD17A3所用的光能利用率模型和CASA模型本身有相似性两种产品之间的高相关性并不能完全证明CASA结果就一定准只能说明两套模型在空间格局上的一致性较好。但反过来说如果CASA结果和MOD17A3存在明显的空间错位或量级差异那大概率是输入数据或者参数设置出了问题值得检查。5.3 三个核心评估指标无论和什么数据对比推荐三个统计指标相关系数R衡量空间格局的相似性R大于0.6算基本可用大于0.8比较理想。均方根误差RMSE反映估算值和观测值的绝对偏差对异常值比较敏感。偏差Bias反映系统性高估或低估正偏差说明模型输出偏高负偏差说明偏低。from scipy import stats def evaluate_model(obs, sim): # 去掉无效点 mask np.isfinite(obs) np.isfinite(sim) obs obs[mask] sim sim[mask] r, p stats.pearsonr(obs, sim) rmse np.sqrt(np.mean((obs - sim) ** 2)) bias np.mean(sim - obs) return {R: r, P值: p, RMSE: rmse, Bias: bias}实际研究中RMSE的量级取决于NPP本身的大小。比如森林年NPP在1000g C/m²左右RMSE 200g C/m²可以接受但若研究区主要是荒漠年NPP只有100g C/m²RMSE 200g C/m²就意味着模型完全不可用。6. 我踩过的坑常见问题与排查技巧实录6.1 数据维度不匹配导致的奇怪结果我第一版代码跑全国模拟时输出的NPP分布图上一大半区域是NaN只有零星几个像素有值。排查了半天发现是一个月份的降水数据与其他月份分辨率不一致在数组乘除时引发了广播错误。这种问题在Python里不像其他语言那样直接报错numpy会自动做广播导致结果错得莫名其妙。排查建议在进入核心计算前写一个断言函数检查所有输入数组的shape是否一致以及分辨率、投影、范围是否匹配。这一步虽然无聊但能省下大量排查时间。6.2 NDVI噪音对FPAR的冲击MODIS NDVI即便做了最大合成在某些高纬度地区或热带雨林区域仍然会残留云污染和水汽噪音。这些噪音会让NDVI瞬时骤降进而导致FPAR和NPP出现明显的凹坑。处理办法有几种最简单的是用Savitzky-Golay滤波对NDVI时间序列做平滑效果比较稳定。也可以用HANTS谐波分析重构NDVI。实测下来S-G滤波配合QA质量波段过滤能解决大部分问题。6.3 极端NPP值的处理运行完模型后一定要做一步异常值检查。NPP理论上应该大于等于0极端退化或火灾干扰时接近0如果出现大面积负值大概率是输入辐射或气温数据有误。如果出现特别大的值比如超过1500g C/m²/年在大多数区域要检查是否在计算中把单位搞错了比如把月辐射当成了年辐射。另外检查输出栅格的极值分布时我习惯画一下直方图和分位数图看到明显的拖尾分布就说明有问题而不是急着把结果扔给合作方。6.4 Python环境与依赖问题处理栅格数据离不开gdal、rasterio、numpy、scipy这几个库。新版Python折腾rasterio的依赖有时候特别痛苦特别是Windows上编译报错很常见。我的建议是直接用Anaconda创建环境用conda安装geopandas、rasterio、xarray这些库能让环境问题少很多。conda create -n casa python3.10 conda activate casa conda install -c conda-forge gdal rasterio xarray netcdf4 scipy matplotlib7. 后续可扩展的三个方向当前这套实现还是经典CASA模型的简化版实际项目中可以根据需求往三个方向扩展。第一耦合土壤水分平衡模块。经典CASA的水分胁迫系数用的是蒸散比如果研究区干旱频发建议耦合一个土壤水分平衡模型让水分胁迫直接由土壤含水量驱动精度会有明显提升。第二精细化植被分类参数。把单一ε_max改成按不同植被类型、不同物候期的参数查找表结合土地利用数据和MODIS物候产品让模型参数更具空间异质性。第三提高时空分辨率。如果计算资源充足可以把时间尺度从月降到旬空间分辨率从1km提高到250m甚至30m。但要做好心理准备计算量和数据管理复杂度都会成倍增加输出结果也要配套做瓦片化或分区存储。我个人在实际操作中的体会是CASA模型的Python实现本身不是什么难题真正的门槛在数据预处理和参数标定。多花一些时间把输入数据弄干净、把参数验证做扎实比单纯调代码更值得投入。最后再分享一个小技巧每次跑完模型把关键中间结果FPAR、T_ε、W_ε都存一份GeoTIFF方便后续排查问题时溯源这个习惯帮我省了无数返工的时间。本文还有配套的精品资源点击获取