尧图网站建设 
简介一套基于Python实现的CASA模型NPP估算脚本面向生态遥感入门者与区域碳汇评估场景解决植被净初级生产力逐像元计算问题。核心文件casa_npp.py封装了光能利用率、APAR植被吸收光合有效辐射等关键计算逻辑支持读取GeoTIFF格式输入并附有示例影像及对应输出结果便于对照验证。代码结构清晰变量命名规范兼容主流Python环境仅需替换输入路径并配置太阳辐射、温度系数等基础参数即可复现流程同时提供create_demo_data.py用于生成演示数据降低测试门槛。资源包共19个文件其中包含14个tif遥感数据、2个Python脚本、依赖列表文本及工程配置文件压缩包大小约20.57MB目录划分明确。目前已有50人学习下载适合作为CASA建模入门参考或轻量级碳汇评估工具。 做区域碳收支评估那阵子我反复被同一件事折腾每次想算植被净初级生产力NPP都要打开不同的软件先做波段运算再套公式出结果换个时间窗口或换一套遥感影像数据就要重来一遍。后来干脆把CASA模型的计算过程整理成一套Python脚本输入示例遥感影像数据就能批量跑出NPP结果从NDVI读入到GeoTIFF输出代码结构不复杂但省下的时间相当可观。这篇帖子就把脚本的公式拆解、数据准备、核心实现和踩坑点完整写出来适合已经有遥感基础、想自己用Python复现CASA模型NPP计算的同学参考。1. CASA模型的公式拆解一个公式串起NPP的天气和植被输入1.1 NPP APAR × ε为什么光能利用率模型这么流行CASA模型属于光能利用率模型家族核心公式非常简洁NPP APAR × ε其中APAR是植被吸收的光合有效辐射ε是实际光能利用率。这个公式之所以能被遥感领域广泛接受是因为它把复杂的植被生产力过程压缩成两个可观测的变量一个与植被结构相关FPAR即植被吸收光合有效辐射的比例一个与环境胁迫相关温度、水分对光能利用率的限制。也就是说只要拿到植被指数影像和气象栅格就能在像元尺度上估算NPP不需要构建过程模型里那些土壤、冠层、光合酶等一大堆参数。CASA里APAR通常这样展开APAR SOL × FPAR × 0.5这里SOL是太阳总辐射0.5是光合有效辐射占太阳总辐射的比例。FPAR则从NDVI反演经典做法是先算SR简单比值SR (1 NDVI) / (1 - NDVI)再用SR极值做线性拉伸FPAR (SR - SR_min) / (SR_max - SR_min) × (FPAR_max - FPAR_min) FPAR_min实际计算中SR_min取1.08SR_max取4.46左右FPAR_min取0.001FPAR_max取0.95。如果NDVI偏低FPAR也可能出现负值所以要加下限保护。ε部分也有固定套路ε ε_max × T_ε1 × T_ε2 × W_ε其中ε_max是最大光能利用率常取0.389 gC/MJ不同植被类型可以调整。T_ε1和T_ε2是温度胁迫系数W_ε是水分胁迫系数。这几个限制因子的设计逻辑很直观植物在最适温度附近光合效率最高偏离太远就会下降水分不足或过少都会压低光能转化效率。1.2 输入数据需求哪些是必须的哪些可以简化根据公式脚本至少需要四类输入NDVI影像、太阳总辐射影像、月均温影像以及一个代表水分状态的指数影像。严格版的CASA还会用到实际蒸散和潜在蒸散但如果手上只有遥感数据可以用土壤湿度指数或从NDVI归一化得到的湿润指数来近似水分胁迫。下面的表格是我预先设置好的输入变量变量含义典型单位示例数据来源NDVI归一化植被指数无量纲-1到1MODIS MOD13Q1SRAD太阳总辐射MJ/m²/dayERA5、气象站插值TMP月平均气温摄氏度MODIS LST、气象站点T_OPT最适温度摄氏度由年内NDVI峰值月份气温确定WET可选水分指数无量纲0到1蒸散比、土壤湿度产品这里特别提醒温度胁迫计算里需要T_OPT它不是生长季最高温度而是NDVI达到年内最大值那一月对应的月均温。这个细节很关键很多初学者直接把夏季最高温填进去结果T_ε1被压得很低。脚本里我会根据输入的NDVI时间序列自动提取T_OPT避免这个坑。2. 数据准备阶段示例影像的目录结构与统一栅格化处理2.1 示例数据目录怎么放拿到一个区域的示例遥感影像后第一件事不是写公式而是把数据整理成固定目录结构。我习惯用一个data目录装原始影像用output目录放结果。示例数据的组织方式如下casapp/ ├── casa_npp.py ├── data/ │ ├── ndvi/ │ │ └── ndvi_2020_07.tif │ ├── srad/ │ │ └── srad_2020_07.tif │ ├── tmp/ │ │ └── tmp_2020_07.tif │ ├── topt/ │ │ └── topt_2020.tif │ └── wet/ │ └── wet_2020_07.tif └── output/这种命名方式的好处是后续脚本可以通过文件夹和日期关键字自动匹配同月的多幅影像。如果你的数据是月度合成产品建议文件名里带上年月比如ndvi_2020_07.tif这样批量处理多个月份时能省去大量手动配置。2.2 统一栅格元数据的核心逻辑遥感影像最麻烦的问题之一是两个栅格的坐标系、范围、分辨率不一致。直接读进numpy就开算哪怕只差一个像元都会导致结果错位。所以我专门写了一个检查函数在执行CASA计算前先统一所有输入栅格的元数据。import numpy as np import rasterio from rasterio.warp import reproject, Resampling def check_and_align(path_dict, reference_pathNone): path_dict: {ndvi: path, srad: path, tmp: path, wet: path, topt: path} 全部统一到参考栅格的网格上返回对齐后的字典 {var: (transform, array)} with rasterio.open(reference_path) as ref: ref_transform ref.transform ref_crs ref.crs ref_width ref.width ref_height ref.height ref_profile ref.profile aligned {} for name, path in path_dict.items(): with rasterio.open(path) as src: if src.crs ! ref_crs or src.transform ! ref_transform or src.width ! ref_width or src.height ! ref_height: dst_array np.zeros((ref_height, ref_width), dtypenp.float32) reproject( sourcerasterio.band(src, 1), destinationdst_array, src_transformsrc.transform, src_crssrc.crs, dst_transformref_transform, dst_crsref_crs, resamplingResampling.bilinear ) else: dst_array src.read(1).astype(np.float32) aligned[name] dst_array return aligned, ref_profile这段代码只用rasterio和numpy没有额外依赖。最容易被忽视的地方是重采样方法NDVI这类连续变量用双线性插值没问题但如果是土地利用类型这种类别型数据就不要用双线性得改成最近邻。本脚本涉及的输入都是连续变量所以统一用bilinear。3. Python核心代码从NDVI到NPP的四步计算逻辑3.1 先算FPAR和APAR真正进入CASA核心计算时我习惯分成四步每步对应一个函数。第一步是根据NDVI计算植被吸收的光合有效辐射比例FPAR再结合太阳辐射得到APAR。注意NDVI影像如果是从MODIS等产品读取的整数型数据很多产品的NDVI做了10000倍缩放读入后要先除以10000否则FPAR可能会跑到几十倍。def calc_fpar(ndvi): 从NDVI计算FPAR ndvi np.where((ndvi 0), 0.0, ndvi) sr (1.0 ndvi) / (1.0 - ndvi 1e-10) fpar (sr - 1.08) / (4.46 - 1.08) * (0.95 - 0.001) 0.001 fpar np.clip(fpar, 0.001, 0.95) return fpar def calc_apar(fpar, srad_mj, days1): 光合有效辐射APAR单位MJ/m²/month apar 0.5 * srad_mj * fpar * days return apar这里有个细节srad输入如果是“日总量”则先乘以本月天数再乘0.5如果srad已经是“月总量”就把days参数设为1。示例脚本我采用日总量乘以月天数的方式这样气象数据更通用。3.2 温度与水分的限制系数第二步是温度胁迫。T_ε1只和最适温度T_OPT有关T_ε2则同时受T_OPT和当月均温TMP影响。公式里有exp指数运算注意分母不能为0我给参数加了保护。def calc_temperature_stress(tmp, topt): 温度胁迫系数 t1 0.8 0.02 * topt - 0.0005 * topt ** 2 exp1 np.exp(0.2 * (topt - 10.0 - tmp)) exp2 np.exp(0.3 * (-topt - 10.0 tmp)) t2 1.1814 / ((1.0 exp1) * (1.0 exp2)) t1 np.clip(t1, 0, 1) t2 np.clip(t2, 0, 1) return t1 * t2 def calc_water_stress(wet): 水份胁迫系数wet范围0-1 w 0.5 0.5 * wet w np.clip(w, 0.1, 1.0) return w水分胁迫这里我做了简化如果只有NDVI数据没有蒸散产品可以直接用归一化植被指数做一个湿润代理比如把NDVI归一化到0到1后作为wet输入。严格版本应该用实际蒸散与潜在蒸散的比值但示例脚本更侧重流程演示所以保留了可替换接口。第三步是把所有系数乘起来得到实际光能利用率εdef calc_epsilon(tstress, wstress, epsilon_max0.389): 实际光能利用率单位gC/MJ epsilon epsilon_max * tstress * wstress return epsilon第四步就是最终NPPdef calc_npp(apar, epsilon): NPP单位gC/m²/month npp apar * epsilon return npp这四个函数就是整个脚本的核心。日常使用中我会根据研究的植被类型微调epsilon_max比如常绿阔叶林可以调到0.485农田有时会更高但这个参数需要有实测或文献支撑不建议随意改。4. 整跑脚本与输出控制把结果写成GeoTIFF4.1 主函数和命令行参数设计为了让脚本可以复用在多个月份我建议把读取、对齐、计算、写出都串到主函数里用argparse接收文件路径。主函数逻辑如下import argparse def main(): parser argparse.ArgumentParser(descriptionCASA模型NPP计算脚本) parser.add_argument(--ndvi, requiredTrue, helpNDVI影像路径) parser.add_argument(--srad, requiredTrue, help太阳总辐射影像路径单位MJ/m2/day) parser.add_argument(--tmp, requiredTrue, help月均温影像路径单位℃) parser.add_argument(--wet, requiredTrue, help水分指数影像路径0-1) parser.add_argument(--topt, requiredTrue, help最适温度影像路径单位℃) parser.add_argument(--out, requiredTrue, help输出NPP GeoTIFF路径) parser.add_argument(--days, typeint, default30, help计算天数用于月尺度累加) parser.add_argument(--eps_max, typefloat, default0.389, help最大光能利用率gC/MJ) args parser.parse_args() path_dict { ndvi: args.ndvi, srad: args.srad, tmp: args.tmp, wet: args.wet, topt: args.topt } aligned, ref_profile check_and_align(path_dict, reference_pathargs.ndvi) ndvi aligned[ndvi].astype(np.float64) srad aligned[srad].astype(np.float64) tmp aligned[tmp].astype(np.float64) wet aligned[wet].astype(np.float64) topt aligned[topt].astype(np.float64) # 处理无效值 valid (ndvi ! -9999) np.isfinite(ndvi) (srad 0) ndvi[~valid] np.nan fpar calc_fpar(ndvi) apar calc_apar(fpar, srad, daysargs.days) tstress calc_temperature_stress(tmp, topt) wstress calc_water_stress(wet) eps calc_epsilon(tstress, wstress, epsilon_maxargs.eps_max) npp calc_npp(apar, eps) out_profile ref_profile.copy() out_profile.update(dtyperasterio.float32, nodatanp.nan, count1) with rasterio.open(args.out, w, **out_profile) as dst: dst.write(np.float32(npp), 1) print(f输出完成: {args.out}) if __name__ __main__: main()命令行运行示例python casa_npp.py \ --ndvi data/ndvi/ndvi_2020_07.tif \ --srad data/srad/srad_2020_07.tif \ --tmp data/tmp/tmp_2020_07.tif \ --wet data/wet/wet_2020_07.tif \ --topt data/topt/topt_2020.tif \ --out output/npp_2020_07.tif \ --days 314.2 输出结果怎么读单位与波段关系输出的NPP单位是gC/m²/month即每月每平方米植被固定的碳克数。空间参考和像元大小与参考栅格完全一致方便后续在ArcGIS或QGIS里做区域统计。如果想换算成年度NPP把12个月结果相加即可如果想换算成kgC/m²/month除以1000即可。有一个容易混淆的地方NPP的栅值出现负值或者0到底正常吗在裸土、水体、城市不透水面等区域NDVI接近0FPAR被clip到0.001NPP会趋近于0这是正常的。但如果大面积出现负值那基本可以断定输入NDVI没有缩放正确或者NoData没有遮罩掉。5. 结果检查与常见错误单位、坐标参考系和Nodata是重灾区5.1 跑完后的合理性检查清单脚本能跑通不代表结果正确。我每算完一个区域都会做三件来验证第一打开NPP输出的直方图。正常情况下森林区域的月NPP应该在100-400 gC/m²/month草地和农田在20-150 gC/m²/month左右。如果全场都超过1000先看是不是srad单位写错了比如把W/m²当成MJ/m²/day用数值会膨胀近百倍。第二把结果和MODIS的MOD17A3年NPP产品做对比。MOD17A3是全球500米分辨率NPP产品虽然不是绝对真值但可以当参考基准。我通常把CASA结果聚合到年尺度再与MOD17A3做相关性分析相关系数在0.6-0.8之间基本说明流程是对的。第三检查边界像元和NoData位置。如果输出影像在原本是NoData的区域出现了明显的异常值多半是遍历性计算时没有排除无效像元。下面这段代码可以快速统计有效像元占比with rasterio.open(output/npp_2020_07.tif) as src: data src.read(1) valid_ratio np.count_nonzero(np.isfinite(data)) / data.size print(f有效像元占比: {valid_ratio:.2%})如果有效像元占比低于80%说明预处理阶段掩膜丢了太多有效区域需要回到对齐函数检查重采样范围。5.2 我最常踩的几个坑第一个坑是坐标系不一致。两个栅格虽然看起来都在同一个经纬度范围但一个用WGS84地理坐标系一个用Web墨卡托投影坐标系rasterio里直接比较crs会报错但如果数据本身是文本格式的.prj坐标系统定义有细微差异直接比较字符串可能不相等。解决办法是在check_and_align里统一用epsg编号判断而不是字符串比较。第二个坑是NoData被当作0参与计算。很多影像产品把背景值设为0NDVI影像的背景0会被误判成稀疏植被FPAR变成0.001NPP跟着变成接近0这样倒不会让结果爆表但会让输出影像的背景区域变成0而不是NoData。后续做面积加权统计时会严重低估总量。我在主函数里统一把所有记录为NoData的值先置为NaN再参与公式计算。第三个坑是SRAD的单位。ERA5下行短波辐射的单位是J/m²有些版本是W/m²的小时数据直接拿来做MJ/m²/day会差几个数量级。我建议在读取SRAD后立刻输出一个像元值做人工核对比如青藏高原夏季晴天某像元的日辐射量应该在20-35 MJ/m²/day这个量级如果读出来是几百万那一定是单位没换算。第四个坑是T_OPT的尺度不匹配。如果TMP是逐月数据T_OPT也应该是同一时期的月均温影像如果TMP是单月数据T_OPT就用多年平均的当月最适温度不要把TMP自身同时当作T_OPT。示例脚本里我为T_OPT单独留了一个参数位目的就是防止这种偷懒导致的系统性误差。最后分享一个我实际改进过的点原始脚本使用固定FPAR参数后来在华南阔叶林区域验证时常绿林的NPP总比MODIS产品低15%左右。我把epsilon_max从0.389调整到0.429并把水分胁迫换成基于土壤湿度产品的蒸散比版本之后相关性提高了不少。这说明CASA模型参数有很强的地方性套用全球默认值只能保证流程通畅做精细研究时必须用实测数据重新标定。后续如果你想按植被类型细分epsilon_max可以用土地利用类型栅格做分区计算在矩阵运算里用np.where按类别选择参数即可逻辑和这套脚本完全兼容。本文还有配套的精品资源点击获取 Python实现CASA模型计算植被净初级生产力（NPP）
简介一套基于Python实现的CASA模型NPP估算脚本面向生态遥感入门者与区域碳汇评估场景解决植被净初级生产力逐像元计算问题。核心文件casa_npp.py封装了光能利用率、APAR植被吸收光合有效辐射等关键计算逻辑支持读取GeoTIFF格式输入并附有示例影像及对应输出结果便于对照验证。代码结构清晰变量命名规范兼容主流Python环境仅需替换输入路径并配置太阳辐射、温度系数等基础参数即可复现流程同时提供create_demo_data.py用于生成演示数据降低测试门槛。资源包共19个文件其中包含14个tif遥感数据、2个Python脚本、依赖列表文本及工程配置文件压缩包大小约20.57MB目录划分明确。目前已有50人学习下载适合作为CASA建模入门参考或轻量级碳汇评估工具。 做区域碳收支评估那阵子我反复被同一件事折腾每次想算植被净初级生产力NPP都要打开不同的软件先做波段运算再套公式出结果换个时间窗口或换一套遥感影像数据就要重来一遍。后来干脆把CASA模型的计算过程整理成一套Python脚本输入示例遥感影像数据就能批量跑出NPP结果从NDVI读入到GeoTIFF输出代码结构不复杂但省下的时间相当可观。这篇帖子就把脚本的公式拆解、数据准备、核心实现和踩坑点完整写出来适合已经有遥感基础、想自己用Python复现CASA模型NPP计算的同学参考。1. CASA模型的公式拆解一个公式串起NPP的天气和植被输入1.1 NPP APAR × ε为什么光能利用率模型这么流行CASA模型属于光能利用率模型家族核心公式非常简洁NPP APAR × ε其中APAR是植被吸收的光合有效辐射ε是实际光能利用率。这个公式之所以能被遥感领域广泛接受是因为它把复杂的植被生产力过程压缩成两个可观测的变量一个与植被结构相关FPAR即植被吸收光合有效辐射的比例一个与环境胁迫相关温度、水分对光能利用率的限制。也就是说只要拿到植被指数影像和气象栅格就能在像元尺度上估算NPP不需要构建过程模型里那些土壤、冠层、光合酶等一大堆参数。CASA里APAR通常这样展开APAR SOL × FPAR × 0.5这里SOL是太阳总辐射0.5是光合有效辐射占太阳总辐射的比例。FPAR则从NDVI反演经典做法是先算SR简单比值SR (1 NDVI) / (1 - NDVI)再用SR极值做线性拉伸FPAR (SR - SR_min) / (SR_max - SR_min) × (FPAR_max - FPAR_min) FPAR_min实际计算中SR_min取1.08SR_max取4.46左右FPAR_min取0.001FPAR_max取0.95。如果NDVI偏低FPAR也可能出现负值所以要加下限保护。ε部分也有固定套路ε ε_max × T_ε1 × T_ε2 × W_ε其中ε_max是最大光能利用率常取0.389 gC/MJ不同植被类型可以调整。T_ε1和T_ε2是温度胁迫系数W_ε是水分胁迫系数。这几个限制因子的设计逻辑很直观植物在最适温度附近光合效率最高偏离太远就会下降水分不足或过少都会压低光能转化效率。1.2 输入数据需求哪些是必须的哪些可以简化根据公式脚本至少需要四类输入NDVI影像、太阳总辐射影像、月均温影像以及一个代表水分状态的指数影像。严格版的CASA还会用到实际蒸散和潜在蒸散但如果手上只有遥感数据可以用土壤湿度指数或从NDVI归一化得到的湿润指数来近似水分胁迫。下面的表格是我预先设置好的输入变量变量含义典型单位示例数据来源NDVI归一化植被指数无量纲-1到1MODIS MOD13Q1SRAD太阳总辐射MJ/m²/dayERA5、气象站插值TMP月平均气温摄氏度MODIS LST、气象站点T_OPT最适温度摄氏度由年内NDVI峰值月份气温确定WET可选水分指数无量纲0到1蒸散比、土壤湿度产品这里特别提醒温度胁迫计算里需要T_OPT它不是生长季最高温度而是NDVI达到年内最大值那一月对应的月均温。这个细节很关键很多初学者直接把夏季最高温填进去结果T_ε1被压得很低。脚本里我会根据输入的NDVI时间序列自动提取T_OPT避免这个坑。2. 数据准备阶段示例影像的目录结构与统一栅格化处理2.1 示例数据目录怎么放拿到一个区域的示例遥感影像后第一件事不是写公式而是把数据整理成固定目录结构。我习惯用一个data目录装原始影像用output目录放结果。示例数据的组织方式如下casapp/ ├── casa_npp.py ├── data/ │ ├── ndvi/ │ │ └── ndvi_2020_07.tif │ ├── srad/ │ │ └── srad_2020_07.tif │ ├── tmp/ │ │ └── tmp_2020_07.tif │ ├── topt/ │ │ └── topt_2020.tif │ └── wet/ │ └── wet_2020_07.tif └── output/这种命名方式的好处是后续脚本可以通过文件夹和日期关键字自动匹配同月的多幅影像。如果你的数据是月度合成产品建议文件名里带上年月比如ndvi_2020_07.tif这样批量处理多个月份时能省去大量手动配置。2.2 统一栅格元数据的核心逻辑遥感影像最麻烦的问题之一是两个栅格的坐标系、范围、分辨率不一致。直接读进numpy就开算哪怕只差一个像元都会导致结果错位。所以我专门写了一个检查函数在执行CASA计算前先统一所有输入栅格的元数据。import numpy as np import rasterio from rasterio.warp import reproject, Resampling def check_and_align(path_dict, reference_pathNone): path_dict: {ndvi: path, srad: path, tmp: path, wet: path, topt: path} 全部统一到参考栅格的网格上返回对齐后的字典 {var: (transform, array)} with rasterio.open(reference_path) as ref: ref_transform ref.transform ref_crs ref.crs ref_width ref.width ref_height ref.height ref_profile ref.profile aligned {} for name, path in path_dict.items(): with rasterio.open(path) as src: if src.crs ! ref_crs or src.transform ! ref_transform or src.width ! ref_width or src.height ! ref_height: dst_array np.zeros((ref_height, ref_width), dtypenp.float32) reproject( sourcerasterio.band(src, 1), destinationdst_array, src_transformsrc.transform, src_crssrc.crs, dst_transformref_transform, dst_crsref_crs, resamplingResampling.bilinear ) else: dst_array src.read(1).astype(np.float32) aligned[name] dst_array return aligned, ref_profile这段代码只用rasterio和numpy没有额外依赖。最容易被忽视的地方是重采样方法NDVI这类连续变量用双线性插值没问题但如果是土地利用类型这种类别型数据就不要用双线性得改成最近邻。本脚本涉及的输入都是连续变量所以统一用bilinear。3. Python核心代码从NDVI到NPP的四步计算逻辑3.1 先算FPAR和APAR真正进入CASA核心计算时我习惯分成四步每步对应一个函数。第一步是根据NDVI计算植被吸收的光合有效辐射比例FPAR再结合太阳辐射得到APAR。注意NDVI影像如果是从MODIS等产品读取的整数型数据很多产品的NDVI做了10000倍缩放读入后要先除以10000否则FPAR可能会跑到几十倍。def calc_fpar(ndvi): 从NDVI计算FPAR ndvi np.where((ndvi 0), 0.0, ndvi) sr (1.0 ndvi) / (1.0 - ndvi 1e-10) fpar (sr - 1.08) / (4.46 - 1.08) * (0.95 - 0.001) 0.001 fpar np.clip(fpar, 0.001, 0.95) return fpar def calc_apar(fpar, srad_mj, days1): 光合有效辐射APAR单位MJ/m²/month apar 0.5 * srad_mj * fpar * days return apar这里有个细节srad输入如果是“日总量”则先乘以本月天数再乘0.5如果srad已经是“月总量”就把days参数设为1。示例脚本我采用日总量乘以月天数的方式这样气象数据更通用。3.2 温度与水分的限制系数第二步是温度胁迫。T_ε1只和最适温度T_OPT有关T_ε2则同时受T_OPT和当月均温TMP影响。公式里有exp指数运算注意分母不能为0我给参数加了保护。def calc_temperature_stress(tmp, topt): 温度胁迫系数 t1 0.8 0.02 * topt - 0.0005 * topt ** 2 exp1 np.exp(0.2 * (topt - 10.0 - tmp)) exp2 np.exp(0.3 * (-topt - 10.0 tmp)) t2 1.1814 / ((1.0 exp1) * (1.0 exp2)) t1 np.clip(t1, 0, 1) t2 np.clip(t2, 0, 1) return t1 * t2 def calc_water_stress(wet): 水份胁迫系数wet范围0-1 w 0.5 0.5 * wet w np.clip(w, 0.1, 1.0) return w水分胁迫这里我做了简化如果只有NDVI数据没有蒸散产品可以直接用归一化植被指数做一个湿润代理比如把NDVI归一化到0到1后作为wet输入。严格版本应该用实际蒸散与潜在蒸散的比值但示例脚本更侧重流程演示所以保留了可替换接口。第三步是把所有系数乘起来得到实际光能利用率εdef calc_epsilon(tstress, wstress, epsilon_max0.389): 实际光能利用率单位gC/MJ epsilon epsilon_max * tstress * wstress return epsilon第四步就是最终NPPdef calc_npp(apar, epsilon): NPP单位gC/m²/month npp apar * epsilon return npp这四个函数就是整个脚本的核心。日常使用中我会根据研究的植被类型微调epsilon_max比如常绿阔叶林可以调到0.485农田有时会更高但这个参数需要有实测或文献支撑不建议随意改。4. 整跑脚本与输出控制把结果写成GeoTIFF4.1 主函数和命令行参数设计为了让脚本可以复用在多个月份我建议把读取、对齐、计算、写出都串到主函数里用argparse接收文件路径。主函数逻辑如下import argparse def main(): parser argparse.ArgumentParser(descriptionCASA模型NPP计算脚本) parser.add_argument(--ndvi, requiredTrue, helpNDVI影像路径) parser.add_argument(--srad, requiredTrue, help太阳总辐射影像路径单位MJ/m2/day) parser.add_argument(--tmp, requiredTrue, help月均温影像路径单位℃) parser.add_argument(--wet, requiredTrue, help水分指数影像路径0-1) parser.add_argument(--topt, requiredTrue, help最适温度影像路径单位℃) parser.add_argument(--out, requiredTrue, help输出NPP GeoTIFF路径) parser.add_argument(--days, typeint, default30, help计算天数用于月尺度累加) parser.add_argument(--eps_max, typefloat, default0.389, help最大光能利用率gC/MJ) args parser.parse_args() path_dict { ndvi: args.ndvi, srad: args.srad, tmp: args.tmp, wet: args.wet, topt: args.topt } aligned, ref_profile check_and_align(path_dict, reference_pathargs.ndvi) ndvi aligned[ndvi].astype(np.float64) srad aligned[srad].astype(np.float64) tmp aligned[tmp].astype(np.float64) wet aligned[wet].astype(np.float64) topt aligned[topt].astype(np.float64) # 处理无效值 valid (ndvi ! -9999) np.isfinite(ndvi) (srad 0) ndvi[~valid] np.nan fpar calc_fpar(ndvi) apar calc_apar(fpar, srad, daysargs.days) tstress calc_temperature_stress(tmp, topt) wstress calc_water_stress(wet) eps calc_epsilon(tstress, wstress, epsilon_maxargs.eps_max) npp calc_npp(apar, eps) out_profile ref_profile.copy() out_profile.update(dtyperasterio.float32, nodatanp.nan, count1) with rasterio.open(args.out, w, **out_profile) as dst: dst.write(np.float32(npp), 1) print(f输出完成: {args.out}) if __name__ __main__: main()命令行运行示例python casa_npp.py \ --ndvi data/ndvi/ndvi_2020_07.tif \ --srad data/srad/srad_2020_07.tif \ --tmp data/tmp/tmp_2020_07.tif \ --wet data/wet/wet_2020_07.tif \ --topt data/topt/topt_2020.tif \ --out output/npp_2020_07.tif \ --days 314.2 输出结果怎么读单位与波段关系输出的NPP单位是gC/m²/month即每月每平方米植被固定的碳克数。空间参考和像元大小与参考栅格完全一致方便后续在ArcGIS或QGIS里做区域统计。如果想换算成年度NPP把12个月结果相加即可如果想换算成kgC/m²/month除以1000即可。有一个容易混淆的地方NPP的栅值出现负值或者0到底正常吗在裸土、水体、城市不透水面等区域NDVI接近0FPAR被clip到0.001NPP会趋近于0这是正常的。但如果大面积出现负值那基本可以断定输入NDVI没有缩放正确或者NoData没有遮罩掉。5. 结果检查与常见错误单位、坐标参考系和Nodata是重灾区5.1 跑完后的合理性检查清单脚本能跑通不代表结果正确。我每算完一个区域都会做三件来验证第一打开NPP输出的直方图。正常情况下森林区域的月NPP应该在100-400 gC/m²/month草地和农田在20-150 gC/m²/month左右。如果全场都超过1000先看是不是srad单位写错了比如把W/m²当成MJ/m²/day用数值会膨胀近百倍。第二把结果和MODIS的MOD17A3年NPP产品做对比。MOD17A3是全球500米分辨率NPP产品虽然不是绝对真值但可以当参考基准。我通常把CASA结果聚合到年尺度再与MOD17A3做相关性分析相关系数在0.6-0.8之间基本说明流程是对的。第三检查边界像元和NoData位置。如果输出影像在原本是NoData的区域出现了明显的异常值多半是遍历性计算时没有排除无效像元。下面这段代码可以快速统计有效像元占比with rasterio.open(output/npp_2020_07.tif) as src: data src.read(1) valid_ratio np.count_nonzero(np.isfinite(data)) / data.size print(f有效像元占比: {valid_ratio:.2%})如果有效像元占比低于80%说明预处理阶段掩膜丢了太多有效区域需要回到对齐函数检查重采样范围。5.2 我最常踩的几个坑第一个坑是坐标系不一致。两个栅格虽然看起来都在同一个经纬度范围但一个用WGS84地理坐标系一个用Web墨卡托投影坐标系rasterio里直接比较crs会报错但如果数据本身是文本格式的.prj坐标系统定义有细微差异直接比较字符串可能不相等。解决办法是在check_and_align里统一用epsg编号判断而不是字符串比较。第二个坑是NoData被当作0参与计算。很多影像产品把背景值设为0NDVI影像的背景0会被误判成稀疏植被FPAR变成0.001NPP跟着变成接近0这样倒不会让结果爆表但会让输出影像的背景区域变成0而不是NoData。后续做面积加权统计时会严重低估总量。我在主函数里统一把所有记录为NoData的值先置为NaN再参与公式计算。第三个坑是SRAD的单位。ERA5下行短波辐射的单位是J/m²有些版本是W/m²的小时数据直接拿来做MJ/m²/day会差几个数量级。我建议在读取SRAD后立刻输出一个像元值做人工核对比如青藏高原夏季晴天某像元的日辐射量应该在20-35 MJ/m²/day这个量级如果读出来是几百万那一定是单位没换算。第四个坑是T_OPT的尺度不匹配。如果TMP是逐月数据T_OPT也应该是同一时期的月均温影像如果TMP是单月数据T_OPT就用多年平均的当月最适温度不要把TMP自身同时当作T_OPT。示例脚本里我为T_OPT单独留了一个参数位目的就是防止这种偷懒导致的系统性误差。最后分享一个我实际改进过的点原始脚本使用固定FPAR参数后来在华南阔叶林区域验证时常绿林的NPP总比MODIS产品低15%左右。我把epsilon_max从0.389调整到0.429并把水分胁迫换成基于土壤湿度产品的蒸散比版本之后相关性提高了不少。这说明CASA模型参数有很强的地方性套用全球默认值只能保证流程通畅做精细研究时必须用实测数据重新标定。后续如果你想按植被类型细分epsilon_max可以用土地利用类型栅格做分区计算在矩阵运算里用np.where按类别选择参数即可逻辑和这套脚本完全兼容。本文还有配套的精品资源点击获取