保姆级教程：用GEE（Google Earth Engine）复现CASA模型计算NPP，附完整数据集清单与Python代码
从零实现CASA模型：GEE平台NPP计算全流程解析与避坑指南
第一次在GEE上跑通CASA模型的那个深夜，我盯着屏幕上终于正确的NPP输出结果，突然意识到——那些看似简单的公式背后，藏着无数个可能出错的细节。本文将带你完整走通这个技术闭环，从数据集选择到代码调试，每个环节都经过实战验证。
1. 环境配置与数据准备
1.1 GEE Python API初始化
确保已安装最新版earthengine-api并完成认证：
!pip install earthengine-api --upgrade import ee ee.Authenticate() # 按提示完成浏览器认证 ee.Initialize()常见报错处理：
EEException: Invalid credentials→ 重新运行ee.Authenticate()TimeoutError→ 检查网络代理设置
1.2 核心数据集清单
当前可用的最新数据组合（2024年验证）：
| 参数 | 数据集 | 版本 | 时间范围 | 替代方案 |
|---|---|---|---|---|
| FPAR | MCD15A3H | v6.1 | 2002-2024 | GLASS FPAR |
| PAR | MCD18C2 | v6.1 | 2002-2024 | CERES SYN1deg |
| 温度 | ERA5-Land | 月度 | 1950-今 | TerraClimate |
| NDVI | MOD09A1 | v6.1 | 2000-今 | Sentinel-2 MSI |
重要提示：MODIS v6.1数据已停止更新，但GEE尚未接入v6.2。若需要最新数据，建议通过NASA Earthdata直接下载后导入GEE Assets。
2. 关键参数计算实战
2.1 NDVI计算优化方案
传统NDVI计算在GEE中的性能瓶颈问题：
def calculate_ndvi(image): # 优选经过大气校正的波段 nir = image.select('sur_refl_b02') red = image.select('sur_refl_b01') ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI') return image.addBands(ndvi) # 使用8天合成数据减少计算量 modis_ndvi = ee.ImageCollection('MODIS/061/MOD09A1') \ .filterDate('2020-01-01', '2020-12-31') \ .map(calculate_ndvi)质量检查技巧：
- 添加云掩膜：
image.updateMask(image.select('State_1km').bitwiseAnd(8).eq(0)) - 季节滤波：
.filter(ee.Filter.calendarRange(6,9,'month'))只保留生长季数据
2.2 PAR数据时空聚合
处理3小时分辨率PAR数据的正确姿势：
def daily_par(image): # 取当地时间12:00的数据（GMT 1200对应各地正午） return image.select('GMT_1200_PAR') par_collection = ee.ImageCollection('MODIS/061/MCD18C2') \ .filterDate('2020-01-01', '2020-12-31') \ .map(daily_par) # 月度聚合时避免简单平均 monthly_par = par_collection.reduce(ee.Reducer.sum()).divide(24)3. CASA模型完整实现
3.1 光能利用率(LUE)计算
整合温度胁迫因子和水分胁迫因子：
def calculate_lue(temp_img, pet_img): # 温度胁迫因子（0-1范围） temp = temp_img.select('temperature_2m').subtract(273.15) # K转℃ t1 = temp.subtract(-10).divide(30) # 经验参数 t_factor = t1.min(1).max(0) # 水分胁迫因子（使用TerraClimate数据） pet = pet_img.select('pet') aet = pet_img.select('aet') w_factor = aet.divide(pet).min(1) return t_factor.multiply(w_factor).rename('LUE') # 获取气候数据 terra_climate = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE') era5_temp = ee.ImageCollection('ECMWF/ERA5_LAND/MONTHLY_AGGR')3.2 NPP最终计算
将各参数整合为CASA公式：
def casa_npp(fpar, par, lue): # 转换单位：PAR从W/m²转换为MJ/m² par_mj = par.multiply(0.0036) return fpar.multiply(par_mj).multiply(lue).multiply(0.5) \ .rename('NPP') # 0.5为碳转换系数 # 组装最终结果 annual_npp = casa_npp( fpar_mean, par_sum, lue_mean )4. 验证与结果导出
4.1 结果可视化检查
创建交互式质量检查工具：
def add_layer_with_scale(collection, band, min_val, max_val, name): vis_params = { 'min': min_val, 'max': max_val, 'palette': ['blue', 'yellow', 'green'] } Map.addLayer(collection.select(band), vis_params, name) # 在GEE地图上叠加检查 Map.centerObject(roi, 8) add_layer_with_scale(annual_npp, 'NPP', 0, 1000, 'Annual NPP')4.2 数据导出策略
针对不同用途的导出方案：
| 用途 | 格式 | 分辨率 | 推荐方法 |
|---|---|---|---|
| 区域分析 | GeoTIFF | 原始 | ee.batch.Export.image.toDrive |
| 时间序列 | CSV | 点位 | ee.batch.Export.table.toDrive |
| 网页展示 | PNG | 缩略图 | getThumbURL |
性能优化技巧：
- 分块导出大区域：
scale=500时单次导出不超过1e8像素 - 使用
clip()限制研究区范围 - 夜间提交任务避免API限额
记得在导出前添加元数据：
annual_npp = annual_npp.set({ 'units': 'gC/m²/year', 'source': 'CASA model v1.0', 'author': 'YOUR_NAME' })