从农田到城市：手把手教你用Landsat 8/9 C2L2 ST数据做地表温度分析（含Python代码示例）

从农田到城市手把手教你用Landsat 8/9 C2L2 ST数据做地表温度分析含Python代码示例当我们需要监测城市热岛效应或评估农田灌溉效率时地表温度数据就像一把精准的测温枪。Landsat 8/9卫星提供的Collection 2 Level 2C2L2地表温度产品ST让这个过程变得前所未有的简单——不需要复杂的大气校正开箱即用的数据产品让研究者能快速投入实际分析。但正如任何遥感数据产品都有其特性ST数据在使用时也藏着不少彩蛋从神秘的块状伪影到突如其来的数据缺失这些坑我已经替你们踩过一遍。今天我们就用Python代码一步步解开这些温度密码。1. 数据获取与预处理1.1 从USGS下载ST数据获取Landsat数据最直接的途径是USGS的EarthExplorer平台。注册账号后在搜索条件中选择数据集Landsat Collection 2 Level-2产品类型L2SP包含地表温度产品时间范围根据研究需求设定典型下载文件结构LC08_L2SP_123045_20230502_20230508_02_T1/ ├── LC08_L2SP_123045_20230502_20230508_02_T1_ANG.txt ├── LC08_L2SP_123045_20230502_20230508_02_T1_MTL.json ├── LC08_L2SP_123045_20230502_20230508_02_T1_MTL.xml ├── LC08_L2SP_123045_20230502_20230508_02_T1_QA_PIXEL.TIF ├── LC08_L2SP_123045_20230502_20230508_02_T1_SR_B1.TIF └── LC08_L2SP_123045_20230502_20230508_02_T1_ST_B10.TIF # 地表温度波段1.2 Python环境配置处理遥感数据需要几个关键库# 基础数据处理 import numpy as np import pandas as pd # 地理数据处理 import rasterio from rasterio.plot import show import geopandas as gpd # 可视化 import matplotlib.pyplot as plt import seaborn as sns # 温度转换常数 KELVIN_OFFSET 273.152. 温度数据核心处理流程2.1 读取ST波段并单位转换Landsat的ST产品默认以开尔文为单位但日常研究更常用摄氏度def read_st_band(st_path): 读取地表温度波段并转换为摄氏度 with rasterio.open(st_path) as src: st_kelvin src.read(1) st_celsius st_kelvin - KELVIN_OFFSET meta src.meta.copy() # 更新元数据 meta.update(dtyperasterio.float32) return st_celsius, meta # 示例使用 st_celsius, meta read_st_band(LC08_L2SP_123045_20230502_ST_B10.TIF)2.2 质量检查与掩膜处理QA波段是遥感数据的质检报告我们需要用它过滤低质量像素def apply_qa_mask(st_array, qa_path): 应用QA波段掩膜 with rasterio.open(qa_path) as src: qa src.read(1) # 创建掩膜示例过滤云、云阴影、雪 cloud_mask (qa 0b0000000011000000) ! 0 cirrus_mask (qa 0b0000110000000000) ! 0 snow_mask (qa 0b0000000000001100) ! 0 # 应用掩膜 masked_st np.ma.masked_where(cloud_mask | cirrus_mask | snow_mask, st_array) return masked_st # 示例使用 clean_st apply_qa_mask(st_celsius, LC08_L2SP_123045_20230502_QA_PIXEL.TIF)3. 伪影识别与异常值处理3.1 典型数据问题识别ST数据常见两类异常块状伪影表现为规则的矩形区域温度异常边缘效应农田与城市交界处的温度跳变伪影检测代码def detect_artifacts(st_array, threshold5): 检测温度异常块状区域 from scipy.ndimage import generic_filter # 计算局部标准差 def std_func(values): return np.std(values) local_std generic_filter(st_array, std_func, size5) # 标记高变异区域 artifacts_mask local_std threshold return artifacts_mask3.2 异常值修正策略针对不同问题需要采用不同处理方法问题类型检测方法修正方案适用场景块状伪影局部标准差分析中值滤波大面积均匀地表边缘效应土地利用对比分区统计城乡过渡带数据缺失NaN值检测空间插值小范围缺失边缘效应修正示例def correct_edge_effect(st_array, landuse_map): 基于土地利用类型修正边缘效应 # 假设landuse_map是相同区域的分类栅格 with rasterio.open(landuse_map) as src: landuse src.read(1) # 按土地利用类型计算区域平均温度 unique_classes np.unique(landuse) class_means {} for cls in unique_classes: class_means[cls] np.nanmean(st_array[landuse cls]) # 替换异常值 corrected st_array.copy() for cls in unique_classes: class_mask landuse cls class_std np.nanstd(st_array[class_mask]) outlier_mask (st_array (class_means[cls] 2*class_std)) | \ (st_array (class_means[cls] - 2*class_std)) corrected[outlier_mask class_mask] class_means[cls] return corrected4. 时空分析与可视化4.1 时间序列分析框架分析多个时相数据时建议构建以下处理流程数据标准化统一空间分辨率匹配坐标系对齐研究区域变化检测def calculate_temp_change(st_stack, dates): 计算温度变化趋势 from sklearn.linear_model import LinearRegression # 重塑数据为时间序列 n_time, rows, cols st_stack.shape X np.array([i for i in range(n_time)]).reshape(-1, 1) # 每个像素拟合线性趋势 trends np.zeros((rows, cols)) for i in range(rows): for j in range(cols): y st_stack[:, i, j] if np.any(np.isnan(y)): trends[i,j] np.nan else: model LinearRegression().fit(X, y) trends[i,j] model.coef_[0] # 变化斜率 return trends4.2 专业级可视化技巧热力图增强技巧def enhanced_st_plot(st_array, extent, title): 增强型地表温度可视化 fig, ax plt.subplots(figsize(10, 8)) # 使用感知均匀的颜色映射 im ax.imshow(st_array, cmapinferno, vminnp.nanpercentile(st_array, 2), vmaxnp.nanpercentile(st_array, 98), extentextent) # 添加专业元素 cbar fig.colorbar(im, axax, shrink0.6) cbar.set_label(地表温度 (°C), rotation270, labelpad20) ax.set_title(title, pad20) ax.grid(True, alpha0.3, linestyle--) # 添加比例尺和指北针 scale_bar(ax, length1000) # 假设1km比例尺 add_north_arrow(ax) return fig5. 典型应用场景实现5.1 城市热岛效应分析热岛强度计算流程划分城市核心区与郊区缓冲区计算各区域平均温度分析温度梯度变化def urban_heat_island(st_array, urban_mask, rural_radius5000): 计算城市热岛强度 from skimage.morphology import disk # 创建郊区掩膜城市边界外扩指定半径 rural_mask binary_dilation(urban_mask, disk(rural_radius/30)) ^ urban_mask # 计算温差 urban_temp np.nanmean(st_array[urban_mask]) rural_temp np.nanmean(st_array[rural_mask]) uhi_intensity urban_temp - rural_temp return uhi_intensity, urban_temp, rural_temp5.2 农田水分胁迫评估水分胁迫指数计算def water_stress_index(st_array, ndvi_array): 基于温度-植被指数关系评估水分胁迫 # 计算温度-植被指数空间 valid_mask (~np.isnan(st_array)) (~np.isnan(ndvi_array)) st_valid st_array[valid_mask] ndvi_valid ndvi_array[valid_mask] # 拟合干湿边 from scipy.stats import linregress bins np.linspace(0, 1, 20) bin_centers (bins[:-1] bins[1:]) / 2 max_temp [np.max(st_valid[(ndvi_valid bins[i]) (ndvi_valid bins[i1])]) for i in range(len(bins)-1)] # 干边拟合 dry_edge linregress(bin_centers, max_temp) # 计算水分胁迫指数 wsi st_array - (dry_edge.slope * ndvi_array dry_edge.intercept) return wsi处理灌溉农田数据时特别要注意ST产品在植被-裸土过渡区域的异常值。有次分析美国中部的农田时发现明明灌溉充足的区域却显示高温异常后来才发现是产品在混合像元处理时的算法特性导致的。这种情况下建议结合NDVI数据先对植被覆盖度进行分层再分别分析各层的温度特征。