Landsat年度合成影像数据集(1985-2023)入门指南：从数据获取到分析实战

1次阅读

共计 1911 个字符，预计需要花费 5 分钟才能阅读完成。

Landsat 年度无云合成影像就像是地球的年度体检报告。1985 年至今的连续观测数据，在以下领域特别有用：

农业监测：通过 NDVI 指数追踪作物长势，辅助产量预估
森林变化：检测非法砍伐或森林再生（比如巴西雨林变化）
城市扩张：量化城镇化进程（参考深圳 1985vs2023 的对比）
冰川退缩：阿拉斯加冰川每年后退几十米的清晰证据

优点：
直接调用 LANDSAT/LT05/C02/T1_L2 等数据集
内置年度合成函数 median() 快速生成无云影像
无需下载，云端计算省流量

# GEE 获取 2000 年合成影像示例代码
import ee
ee.Initialize()

collection = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2') \
    .filterDate('2000-01-01', '2000-12-31') \
    .median() \
    .select('SR_B.')  # 选择地表反射率波段

优点：
原始数据最全，包含 L1TP 级精密校正数据
可批量下载整景影像
缺点：
需要手动拼接和云过滤
下载速度较慢（建议用 wget 批量下载）

关键两步：

DN 值转 TOA 反射率（大气顶层）

# rasterio 处理示例
with rasterio.open('LT05_2000.tif') as src:
    dn = src.read(3)  # 读取红波段
    radiance = dn * 0.0003342 + 0.1  # Landsat5 的增益 / 偏置参数
    toa = radiance / np.cos(sun_zenith_rad)  # 太阳高度角校正

转地表反射率(SR) – 使用 LEDAPS 或 LaSRC 算法

推荐使用 QA 波段：

# 提取云掩膜（以 Landsat8 为例）qa = src.read(9)  # QA 波段
cloud_mask = (qa & 0b10000000) > 0  # 高位云标志
clean_data = np.where(cloud_mask, np.nan, toa)  # 云区置空

# 计算 2000 年 NDVI
with rasterio.open('LT05_2000_B3.tif') as red, \
     rasterio.open('LT05_2000_B4.tif') as nir:

    ndvi = (nir.read(1) - red.read(1)) / (nir.read(1) + red.read(1))

    # 保存结果
    profile = red.profile
    profile.update(dtype=rasterio.float32)
    with rasterio.open('ndvi_2000.tif', 'w', **profile) as dst:
        dst.write(ndvi.astype(np.float32), 1)

# 计算植被减少区域
ndvi_diff = ndvi_2020 - ndvi_2000
loss_area = (ndvi_diff < -0.2).sum() * 30 * 30  # 30m 分辨率转面积㎡

# 512x512 分块读取
block_shape = (512, 512)
for ji, window in rasterio.windows.WindowMethods.block_windows(src.transform, *block_shape):
    chunk = src.read(window=window)
    # 处理分块数据...

import dask.array as da

data = da.from_zarr('landsat_stack.zarr', chunks=(256,256))
ndvi = (data[3] - data[2]) / (data[3] + data[2])  # 延迟计算
ndvi.compute(num_workers=4)  # 启动 4 核并行