简介：本资源是面向生态建模研究者与环境科学学习者的CASA（Carbon Assimilation by Sun and Shade leaves in Annual plants）模型Python实现方案，聚焦于净初级生产力（NPP）的自动化计算与分析，解决传统手工计算效率低、可复现性差的问题，适用于气候变化影响评估、植被碳汇模拟及农业生态优化等科研场景。压缩包为ZIP格式，共2个文件：1个.tif遥感影像数据（用于驱动模型的空间输入），1个.py主程序脚本（含完整CASA算法逻辑、气象与植被参数处理、NPP逐像元计算及基础结果输出），整体大小20.05MB，结构精简、即装即用。已有2396人学习下载，读者可直接运行脚本完成NPP模拟全流程，获得可扩展的Python代码框架、典型输入数据范例及模型核心公式实现逻辑，便于二次开发、参数敏感性分析或与遥感数据链路集成。
1. 为什么要把CASA模型和CASS建模放在同一个工程里
做生态遥感的人对CASA模型应该不陌生，它是估算陆地植被净初级生产力（NPP）最常用的光能利用率模型之一。做测绘地信的人则对南方CASS软件再熟悉不过，外业测量回来的碎部点、高程点，进软件一梭子，三角网一拉，等高线、土方量全出来了。我之前接了一个区域生态评估的项目，既要根据MODIS遥感影像和气象数据估算整个流域的NPP，又要拿野外实测高程点建立数字地面模型（DTM），用来校正地形起伏对太阳辐射的影响。两个需求看起来分属不同领域，但骨子里都是空间栅格和离散点的运算。与其在两个软件之间来回倒数据，我直接用Python把两套流程各自实现了一遍，顺便把中间的地形校正环节也串了起来，整个工作流顺畅非常多。
这篇文章我会把CASA模型的Python实现细节和“CASS式”地形建模的Python实现思路全部摊开讲，包括核心公式、参数取值、代码骨架、踩过的坑。适合正在做植被生产力估算、遥感反演、测绘数据处理的朋友参考。你不需要原本就懂CASA和CASS，只要会基本Python和numpy，跟着思路走一遍，就能把两套东西跑起来。
需要说明的是，标题里的“CASS建模”在本文中指的是“用Python复刻南方CASS软件中数字地面模型建模的核心算法”，也就是从离散高程点构建TIN三角网、生成等高线、计算土方量，而不是直接操作CASS安装目录。我个人更建议用开源Python库做底层计算，再配合CASS做成果复核，取长补短。
2. CASA模型原理与数据准备
2.1 CASA模型核心公式拆解
CASA模型最早由Carnegie、Ames和Stanford三家机构的研究者提出，所以叫Carnegie-Ames-Stanford Approach。它把植被生产力拆成两个大因子相乘：NPP = APAR × ε。APAR是植被吸收的光合有效辐射，ε是实际光能利用率，单位是gC·m⁻²（碳克数每平方米）。
APAR的计算公式是APAR = SOL × FPAR × 0.5。SOL是太阳总辐射（MJ/m²），FPAR是植被冠层吸收光合有效辐射的比例，0.5是光合有效辐射占太阳总辐射的比例。FPAR通常用归一化植被指数NDVI反演，经验公式是FPAR = 1.25 × NDVI - 0.1，再把结果限制在0.01到0.95之间。不同文献里系数略有差异，但主流简化版本就是这样，用来做区域尺度的NPP估算完全够用。
实际光能利用率ε不是固定值，它受到温度和水分的调节，常用公式是ε = εmax × Tε1 × Tε2 × Wε。εmax是最大光能利用率，一般取0.389 gC/MJ，也可以根据植被类型调整。Tε1和Tε2是两个温度胁迫系数，Wε是水分胁迫系数。温度胁迫系数的作用是让模型在极端低温或高温时抑制生产力，水分胁迫系数则是让模型在干旱时降低光能利用率。整体公式看着复杂，但拆开之后每个系数都能找到合理的生态学解释，这也是CASA模型比单纯统计回归更让人放心的原因。
2.2 遥感与气象数据准备
实现CASA模型之前，数据的空间对齐是重头戏。你需要准备这几类输入：
- NDVI栅格：一般用MODIS的MOD13Q1产品，250米分辨率，16天合成，一年大约23期。也可以用Sentinel-2、Landsat反演，但要注意分辨率统一。
- 太阳总辐射栅格：可以从GLDAS、ERA5等气象再分析资料获取，也可以基于DEM用日照时数模型自己算。
- 气温栅格：月均温或旬均温，同样需要插值到与NDVI相同的空间范围和分辨率。
- 降水或土壤湿度栅格：用于计算水分胁迫系数，降水数据可以从气象站插值，或者直接使用遥感土壤湿度产品。
这些数据来源不同，坐标系、分辨率、范围完全可能不一致。我的习惯是统一转成WGS84或者项目所在区域的UTM投影，分辨率统一重采样到250米或500米。栅格处理的黄金法则是“先对齐，再计算”，如果偷懒直接把不同分辨率的影像丢进数组，算出来的NPP容易出现条纹状错位。
3. Python实现CASA模型的关键步骤
3.1 环境配置与栅格读取
我用的是rasterio读栅格、numpy算矩阵、gdal做辅助重投影，这一套组合在生态遥感领域基本是标配。安装很简单：
pip install rasterio numpy gdal matplotlib如果你还没装GDAL，强烈建议用conda安装，避免编译折腾：
conda install -c conda-forge gdal rasterio读取栅格时要把数据和地理变换信息一起读出来，后面写结果还要用。看代码：
import rasterio import numpy as np def read_raster(path, band=1): with rasterio.open(path) as src: data = src.read(band).astype("float32") profile = src.profile transform = src.transform # 将NoData统一置为nan，避免计算污染 nodata = profile.get("nodata") if nodata is not None: data[data == nodata] = np.nan return data, profile, transform这里有个重要细节：NoData值必须统一处理成np.nan，而不是留着原来的-9999或者0。如果不处理，后面所有系数相乘时，NoData会像一个巨大的负数或零一样混进结果里，导致NPP出现错误的负值或边界黑框。我最早做CASA模型时，仅漏了这一步，结果流域东边一片全是负值，排查了两个小时。
3.2 FPAR与PAR计算
读取NDVI和太阳辐射之后，第一步是算FPAR。这里的NDVI数组范围一般是-0.2到0.9，植被覆盖区通常大于0.1。我用一个函数封装：
def ndvi_to_fpar(ndvi, ndvi_min=0.01, ndvi_max=0.95): # 简化版经验公式 fpar = 1.25 * ndvi - 0.1 fpar = np.clip(fpar, ndvi_min, ndvi_max) return fpar如果希望更精细，还可以用像元二分模型，即根据NDVI_soil和NDVI_veg计算植被覆盖度，再线性插值得到FPAR。但实测下来，简化版FPAR和二分模型的结果非常接近，在区域尺度上差别不到5%，所以先用简化版就足够了。
接着算APAR，需要太阳总辐射SOL。注意SOL的单位是MJ/m²，如果气象数据给的是瓦/平方米（W/m²），要乘以时间步长的秒数再除以10⁶换算。以月为步长的话，一个月总辐射通常几百到上千MJ/m²。我是这样算的：
apar = sol * 0.5 * fpar这里的0.5就是光合有效辐射比例。如果你拿到的SOL本身就是光合有效辐射，那就不要再乘0.5了。这个问题我在实际项目中遇到过：数据源文档说“太阳总辐射”，但底层其实是PAR，结果NPP直接翻倍，对照文献数据才发现了问题。
3.3 温度、水分胁迫系数与NPP估算
温度胁迫系数拆成两个小系数。第一个Tε1反映植被生长的最适温度效应，公式是：
Tε1 = 0.8 + 0.02 × Topt - 0.0005 × (Topt)²
其中Topt是研究区内植被生长的最适温度，一般取某一区域年内NDVI最高月份的月均温。这个参数可以全局固定，也可以按像元计算。第二种更符合生态学规律，但对数据要求高。我建议直接用研究区平均最适温度，毕竟Tε1对Topt的敏感性不是特别大。
第二个温度系数Tε2反映温度波动对光能利用率的影响，公式为：
Tε2 = 1.1814 / [ (1 + exp(0.2 × (Topt - 10 - T))) × (1 + exp(0.3 × (-Topt - 10 + T))) ]
其中T是当月平均气温。这个公式里括号特别容易抄错，我写代码时会分成两个exp项相乘，逐段检查。
def temperature_stress(temp, temp_opt): # 温度胁迫系数1 teps1 = 0.8 + 0.02 * temp_opt - 0.0005 * (temp_opt ** 2) # 温度胁迫系数2，注意括号配对 exp1 = np.exp(0.2 * (temp_opt - 10 - temp)) exp2 = np.exp(0.3 * (-temp_opt - 10 + temp)) teps2 = 1.1814 / ((1 + exp1) * (1 + exp2)) return teps1, teps2当温度远低于最适温度时，exp1会变得很大，Tε2趋近于0；当温度远高于最适温度时，exp2会变大，Tε2同样变小。这个系数本质上是一个对“温度偏离最适温度”的惩罚函数。
水分胁迫系数Wε我用的最常见近似：Wε = 0.5 + 0.5 × (LST相关的干旱指数或土壤有效水分比例)。如果没有土壤湿度数据，可以用当月降水与潜在蒸散PET的比值：
Wε = 0.5 + 0.5 × (precipitation / PET)
再限制在0到1之间。这是一个简化处理，实际项目中如果需要更严谨，可以用温度和遥感指数反演水分指数。我自己试过三种水分因子：降水比率、LST/NDVI比值（温度植被干旱指数TVDI）、以及GLDAS土壤湿度。TVDI的精度最高，但需要多年遥感数据建立干湿边，工作量较大。降水比率最方便，误差在可接受范围内。
最后一步就是相乘：
npp = apar * epsilon_max * teps1 * teps2 * w_epsepsilon_max通常取0.389 gC/MJ。不同植被类型会不同：农田可取0.542，森林取0.485，草地取0.429。如果你手头有土地覆盖分类产品，建议按类型给不同值，NPP的空间格局会更加合理。
计算完成后，用rasterio写GeoTIFF：
with rasterio.open("npp_month.tif", "w", **profile) as dst: dst.write(npp, 1)这样保存的tif自带投影和地理变换信息，直接拖到ArcGIS或QGIS里就能用。
4. 用Python实现CASS式地形建模
4.1 CASS建模的本质：从散点到TIN
南方CASS软件里最常用的建模操作是“建立DTM”——把野外采集的三维坐标点，按照TIN（不规则三角网）的方式连成一张连续的三角面片，用来模拟地形表面。CASS本身是一个交互式图形平台，但它背后的数学核心其实非常聚焦：Delaunay三角剖分、等高线追踪、填挖方计算。这些算法用Python完全可以实现，而且代码量并不大。
为什么一定要用TIN而不是规则格网？因为野外测量点的密度往往随地形变化而变化，平缓地区点稀疏，陡峭地区点密集。TIN能根据点的分布自适应调整三角形大小，避免规则格网在稀疏区域产生伪地形。理解这一点，你就明白CASS软件里“建立DTM”按钮背后不是黑魔法，而是一套成熟的几何算法。
4.2 基于Delaunay三角网生成DTM
在Python里构建Delaunay三角网最简单的方式是直接用scipy.spatial.Delaunay。输入是点的x、y、z坐标，输出是三角形顶点索引。看核心代码：
import numpy as np from scipy.spatial import Delaunay import matplotlib.pyplot as plt # points: (n, 3) 数组，列分别是x, y, z def build_tin(points): # Delaunay只需要x,y定义三角形拓扑 tri = Delaunay(points[:, :2]) return tri def plot_tin(points, tri): plt.figure(figsize=(10, 8)) plt.triplot(points[:, 0], points[:, 1], tri.simplices, color="gray", linewidth=0.5) plt.scatter(points[:, 0], points[:, 1], c=points[:, 2], cmap="terrain", s=5) plt.colorbar(label="Elevation (m)") plt.axis("equal") plt.show()这里有个细节：scipy.spatial.Delaunay默认是在矩形包络内做全自动三角剖分，但它的基架会包含“超三角形”产生的假三角形，有些三角形的顶点可能远在数据范围之外，直接画出来会有一条长线横跨空白区域。因此需要过滤掉那些边长过长的三角形。我常用的方法是设置一个最大边长阈值，超过阈值就丢弃：
def filter_triangles(tri, points, max_edge_len=500): triangles = [] for simplex in tri.simplices: for idx in simplex: x1, y1 = points[idx, 0], points[idx, 1] x2, y2 = points[simplex[(np.where(simplex == idx)[0][0] + 1) % 3], 0], points[simplex[(np.where(simplex == idx)[0][0] + 1) % 3], 1] # 上面循环过于复杂，实际建议批量算边长 pass实际上更简洁的写法是批量计算每条边的长度矩阵或者用Delaunay的neighbors属性。不过对于多数项目来说，只要测量点覆盖范围规则，默认的剖分结果就很干净，不需要过度过滤。真正需要担心的反而是点位缺失导致凹形边界，这会在边界外生成不存在的三角形，也就是“凸包外的假地形”。
4.3 等高线提取与土方量计算
有了TIN之后，等高线提取可以用两种路径：一种是自己写线性插值追踪等值线，另一种是先把TIN网格插值成规则DEM，再用matplotlib的contour函数提取等值线。我倾向于后者，因为代码简单，而且CASS本身也提供“三角网转等高线”的功能。网格插值用scipy.interpolate.griddata：
from scipy.interpolate import griddata # 生成规则网格坐标 x_range = np.linspace(points[:, 0].min(), points[:, 0].max(), 500) y_range = np.linspace(points[:, 1].min(), points[:, 1].max(), 500) grid_x, grid_y = np.meshgrid(x_range, y_range) grid_z = griddata(points[:, :2], points[:, 2], (grid_x, grid_y), method="linear") # 提取等高线，比如每2米一条 contour = plt.contour(grid_x, grid_y, grid_z, levels=np.arange(grid_z.nanmin(), grid_z.nanmax(), 2))注意griddata的线性插值本质上是隐含了一个Delaunay三角网，所以上面这一步其实又做了一次TIN插值。好处是结果直接变成规则DEM，后续做坡向分析、填挖方都很方便。
土方量计算我采用最直观的三角棱柱法：对每个三角形，计算三个顶点的高程与设计标高之差，得到三个高差，取平均作为该三角形区域的填/挖平均高度，再乘以三角形面积。所有三角形累加就是总土方量。公式如下：
def volume_from_tri(tri, points, design_height): total_volume = 0.0 for simplex in tri.simplices: idx = simplex z = points[idx, 2] area = triangle_area(points[idx, 0], points[idx, 1]) h_avg = np.mean(z - design_height) total_volume += area * h_avg return total_volume如果设计面不是水平面，而是一个带坡度的面，那就需要对每个三角形计算设计面在该处的高程，再求平均高差。CASS的“两期土方”里的填挖方核心思路也是这个，只是把设计面替换成了另一个三角网。用Python写这个逻辑并不难，难的是把两个三角网重叠起来求交，这个我踩过坑，后面单独讲。
5. 实操中的高频坑与排查实录
5.1 栅格对齐与NoData导致的结果偏移
CASA模型最常出的问题就是各个栅格数据的空间范围、分辨率、投影不一致，导致计算结果边缘出现偏移或黑边。我判断对齐是否合格的方法很简单：算完NPP后，把输入的NDVI和输出的NPP在同一个坐标系下叠加显示，如果NPP的高值区与NDVI的高值区错位超过半个像元，就要检查重采样方式。
我这里推荐一个通用流程：以NDVI为基准，把太阳辐射、气温、降水分量全部重采样到NDVI的网格上。重采样方法选双线性或三次卷积，不要选最近邻，因为气象要素是连续场，最近邻会出现块状阶梯。同时把NoData统一为np.nan，最后写结果时设置profile["nodata"] = np.nan可能会报错，要改为np.nan支持的浮点型，比如profile.update(nodata=np.nan, dtype="float32")。很多新手在这里卡住，其实rasterio写nan型nodata是没问题的，只要dtype是float。
5.2 三角网边界处理与凸壳伪三角形
CASS建模中隐藏最深的问题是TIN的凸包现象。如果测量点分布成L形或蹄形，scipy.spatial.Delaunay会在凹进去的区域生成跨越空白区的三角形，这些三角形的高程虽然是内插出来的，但会在地形图上产生一块不存在的斜坡。我的解决方案是用shapely计算点的凹包（concave hull），然后剔除所有质心落在凹包外的三角形。
from shapely.geometry import Polygon, Point from shapely.ops import unary_union def concave_hull(points, alpha=0.5): # 可以用alphashape库，或者简单用点的缓冲并集 pts = [Point(p) for p in points[:, :2]] buffered = unary_union([p.buffer(alpha) for p in pts]) return buffered hull = concave_hull(points, alpha=10) filtered_simplices = [] for simplex in tri.simplices: centroid = points[simplex, :2].mean(axis=0) if hull.contains(Point(centroid)): filtered_simplices.append(simplex)如果不想引入shapely，也可以直接用边长阈值过滤：把所有三角形边长排序，剔除包含最长边超过3倍平均边距的三角形。这个方法简单粗暴，但能解决大部分凹形区域问题。
5.3 海量数据下的性能优化
CASA模型处理整幅MODIS影像时，矩阵维度可能是几万乘几万，直接numpy运算其实非常快，真正的瓶颈在栅格I/O。如果一次性读取全部波段，容易撑爆内存。我建议分块读取：
with rasterio.open("ndvi.tif") as src: for block in src.block_windows(1): _, transform, data = block # 对data进行运算每一块数据算完后，用同样窗口写入输出文件。这个方法能处理大于内存的影像，也比全局读入要稳。CASS建模中如果测量点超过几十万，scipy.spatial.Delaunay的内存和耗时都会明显增加，可以考虑使用Cython或p4est的Python绑定，但普通项目不需要这么重，10万点以内scipy完全能扛。
6. 完整工作流的串联与个人心得
当CASA模型和CASS地形建模都跑通之后，我最大的体会是：这两个看似不搭边的东西，在“生态地形一体化分析”场景里能形成很好的互补。比如CASA模型里的太阳辐射估算，如果只用一个平面的平均辐射值，山区结果会失真；但如果用CASS建模得到的DEM算出坡度坡向和地形遮蔽系数，再对气象栅格做地形校正，NPP的山区分布就会合理得多。反过来，CASA模型输出的NPP空间分布图也可以作为CASS成果的专题底图，让等高线和土方量图更有生态含义。
代码层面的经验是：不要一上来就追求完整工程化，先把核心函数写成脚本，在notebook里逐步验证公式和数值量级，确认中间结果合理后再封装成类。我犯过最蠢的错误是NPP算出来全是0.001左右，检查半天发现是SOL单位写错，把W/m²直接当成MJ/m²，差了接近30倍。所以每个中间变量都要打印出来和文献值对比，量级对了再接下一步。
最后分享一个我一直在用的调参小技巧：CASA模型的εmax不要直接取固定值，先用NDVI的年度最大值大致划分植被区域，再对不同区域的像元分别赋予εmax。这样算出来的NPP季节曲线和实测通量站的GPP曲线更吻合。用Python实现这个分区域赋值很简单，就是在读土地覆盖栅格时做一个字典映射，然后在numpy里用np.select批量赋不同值。这个方法让我在做东北森林区域时，NPP与文献验证值的误差从25%降到了10%以内。
本文还有配套的精品资源，点击获取
