简介本资源是面向生态建模研究者与环境科学学习者的CASACarbon 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, band1): 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_min0.01, ndvi_max0.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写GeoTIFFwith 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, colorgray, linewidth0.5) plt.scatter(points[:, 0], points[:, 1], cpoints[:, 2], cmapterrain, s5) plt.colorbar(labelElevation (m)) plt.axis(equal) plt.show()这里有个细节scipy.spatial.Delaunay默认是在矩形包络内做全自动三角剖分但它的基架会包含“超三角形”产生的假三角形有些三角形的顶点可能远在数据范围之外直接画出来会有一条长线横跨空白区域。因此需要过滤掉那些边长过长的三角形。我常用的方法是设置一个最大边长阈值超过阈值就丢弃def filter_triangles(tri, points, max_edge_len500): 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.griddatafrom 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), methodlinear) # 提取等高线比如每2米一条 contour plt.contour(grid_x, grid_y, grid_z, levelsnp.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(nodatanp.nan, dtypefloat32)。很多新手在这里卡住其实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, alpha0.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, alpha10) filtered_simplices [] for simplex in tri.simplices: centroid points[simplex, :2].mean(axis0) 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%以内。本文还有配套的精品资源点击获取