Open Geospatial Engine教程：Landsat 8遥感计算EVI与LAI全流程

大家好，我是讯享网，很高兴认识大家。这里提供最前沿的Ai技术和互联网信息。

以下是将代码转化为的中文博客文章，包含详细的步骤解释和代码展示。

叶面积指数（LAI）是描述植被冠层结构的重要参数，广泛用于生态、农业和气候变化研究。利用卫星遥感影像，我们可以通过植被指数反演 LAI。本文将在开放地球引擎（OGE）中，以 Landsat 8 影像为例，演示如何：

对原始影像进行辐射定标和太阳高度角校正；
计算增强型植被指数（EVI）；
利用经验模型将 EVI 转换为 LAI；
可视化最终结果。

首先，导入 oge 包，初始化引擎并创建一个 Service 对象。

import oge

oge.initialize() service = oge.Service.initialize()

我们使用 Landsat 8 L1T 级地表反射率产品（产品 ID：LC08_L1T），示例影像覆盖武汉地区，拍摄日期为 2015 年 10 月 2 日。加载后转换为浮点型，便于后续运算。

lc08 = service.getCoverage(coverageID=“LC275LGN00”, productID=“LC08_L1T”) lc08_float = service.getProcess(“Coverage.toFloat”).execute(lc08)

EVI 需要红（B4）、近红外（B5）和蓝（B2）三个波段。原始 Landsat 8 L1T 产品中的数值为无量纲的 DN 值，我们需要通过增益（Gain）和偏移（Offset）将其转换为表观反射率，再乘以太阳高度角校正系数（1/sinθ）。

根据元数据，本例使用的系数为：

增益系数：2.0000E-05（所有波段相同）
偏移系数：-0.（所有波段相同）
太阳高度角校正系数：1.

# 提取波段 red = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B4”]) nir = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B5”]) blue = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B2”])

定义系数

Gain = 2.0000E-05 Offset = -0. sun_corr = 1.

校正蓝光波段

B2 = service.getProcess(“Coverage.multiplyNum”).execute(blue, Gain) B2 = service.getProcess(“Coverage.addNum”).execute(B2, Offset) B2 = service.getProcess(“Coverage.multiplyNum”).execute(B2, sun_corr)

校正红光波段

B4 = service.getProcess(“Coverage.multiplyNum”).execute(red, Gain) B4 = service.getProcess(“Coverage.addNum”).execute(B4, Offset) B4 = service.getProcess(“Coverage.multiplyNum”).execute(B4, sun_corr)

校正近红外波段

B5 = service.getProcess(“Coverage.multiplyNum”).execute(nir, Gain) B5 = service.getProcess(“Coverage.addNum”).execute(B5, Offset) B5 = service.getProcess(“Coverage.multiplyNum”).execute(B5, sun_corr)

EVI 相比 NDVI 更能减少大气和土壤背景的影响，尤其适合高植被覆盖区。其计算公式为：

ext{EVI} = G imes frac{ ext{NIR} - ext{Red}}{ ext{NIR} + C_1 imes ext{Red} - C_2 imes ext{Blue} + L}

]

其中经验常数取值为：

( G = 2.5 )（增益因子）
( C_1 = 6 )，( C_2 = 7.5 )（大气修正系数）
( L = 1 )（土壤亮度修正值）

G = 2.5 L = 1 C1 = 6 C2 = 7.5

分子: NIR - Red

numerator = service.getProcess(“Coverage.subtract”).execute(B5, B4)

分母: NIR + C1Red - C2Blue + L

red_scaled = service.getProcess(“Coverage.multiplyNum”).execute(B4, C1) blue_scaled = service.getProcess(“Coverage.multiplyNum”).execute(B2, C2) denominator = service.getProcess(“Coverage.subtract”).execute(red_scaled, blue_scaled) denominator = service.getProcess(“Coverage.add”).execute(denominator, B5) denominator = service.getProcess(“Coverage.addNum”).execute(denominator, L)

EVI = G * (分子/分母)

evi_ratio = service.getProcess(“Coverage.divide”).execute(numerator, denominator) evi = service.getProcess(“Coverage.multiplyNum”).execute(evi_ratio, G)

我们采用一个简单的线性经验模型，将 EVI 转换为 LAI：

ext{LAI} = 3.618 imes ext{EVI} - 0.118

]

该系数来源于对多类植被的统计回归。

LAI = service.getProcess(“Coverage.multiplyNum”).execute(evi, 3.618) LAI = service.getProcess(“Coverage.subtractNum”).execute(LAI, 0.118)

设置 LAI 的显示范围为 0～10，使用从蓝色到深红色的渐变色带，以便直观识别不同植被密度。

vis_params = {

"min": 0, "max": 10, "palette": [ "#0000FF", # 深蓝色 -> 低 LAI "#00FFFF", # 浅蓝色 "#00FF00", # 绿色 "#FFFF00", # 黄色 "#FF0000", # 红色 "#" # 深红色 -> 高 LAI ]

}

LAI.styles(vis_params).getMap(“LAI”)

最后将地图视图定位到武汉（经度 114.28°，纬度 30.57°），缩放级别设为 7。

oge.mapclient.centerMap(114.28, 30.57, 7)

import oge

================== 1. 初始化 ==================

oge.initialize() service = oge.Service.initialize()

================== 2. 加载 Landsat 8 数据 ==================

lc08 = service.getCoverage(coverageID=“LC275LGN00”, productID=“LC08_L1T”) lc08_float = service.getProcess(“Coverage.toFloat”).execute(lc08)

================== 3. 提取相关波段 ==================

选择红波段 (B4)

red = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B4”])

选择近红外波段 (B5)

nir = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B5”])

选择蓝波段 (B2)

blue = service.getProcess(“Coverage.selectBands”).execute(lc08_float, [“B2”])

Gain2 = 2.0000E-05 # 蓝光波段系数 Gain4 = 2.0000E-05 # 红光波段系数 Gain5 = 2.0000E-05 # 近红外波段系数 Offset2 = -0. # 蓝光波段偏移 Offset4 = -0. # 红光波段偏移 Offset5 = -0. # 近红外波段偏移 sun_corr = 1. #太阳高度角矫正系数 B2 = service.getProcess(“Coverage.multiplyNum”).execute(blue, Gain2) B2 = service.getProcess(“Coverage.addNum”).execute(B2, Offset2) B2 = service.getProcess(“Coverage.multiplyNum”).execute(B2, sun_corr) B4 = service.getProcess(“Coverage.multiplyNum”).execute(red, Gain4) B4 = service.getProcess(“Coverage.addNum”).execute(B4, Offset4) B4 = service.getProcess(“Coverage.multiplyNum”).execute(B4, sun_corr) B5 = service.getProcess(“Coverage.multiplyNum”).execute(nir, Gain5) B5 = service.getProcess(“Coverage.addNum”).execute(B5, Offset5) B5 = service.getProcess(“Coverage.multiplyNum”).execute(B5, sun_corr)

================== 4. 计算 EVI ==================

增益因子

G = 2.5

土壤亮度修正值

L = 1

常数 C1 和 C2

C1 = 6 C2 = 7.5

计算 EVI: EVI = G * ((NIR - Red) / (NIR + C1 * Red - C2 * Blue + L))

numerator = service.getProcess(“Coverage.subtract”).execute(B5, B4) red_scaled = service.getProcess(“Coverage.multiplyNum”).execute(B4, C1) blue_scaled = service.getProcess(“Coverage.multiplyNum”).execute(B2, C2) denominator = service.getProcess(“Coverage.subtract”).execute(red_scaled, blue_scaled) denominator = service.getProcess(“Coverage.add”).execute(denominator, B5) denominator = service.getProcess(“Coverage.addNum”).execute(denominator, L) evi_ratio = service.getProcess(“Coverage.divide”).execute(numerator, denominator) evi = service.getProcess(“Coverage.multiplyNum”).execute(evi_ratio, G)

LAI = service.getProcess(“Coverage.multiplyNum”).execute(evi, 3.618) LAI = service.getProcess(“Coverage.subtractNum”).execute(LAI, 0.118)

================== 5. 可视化 EVI ==================

设置可视化参数

vis_params = {

"min": 0, "max": 10, "palette": [ "#0000FF", # 深蓝色 "#00FFFF", # 浅蓝色 "#00FF00", # 绿色 "#FFFF00", # 黄色 "#FF0000", # 红色 "#" # 深红色 ]

}

显示 EVI 图层