简介本资源是一份面向遥感与Python编程初学者的Landsat 8地表温度反演实践代码包聚焦热红外遥感定量分析核心任务适用于环境监测、城市热岛研究及遥感教学等场景。压缩包仅含2个文件1个Python主脚本SC_LSTR.py 1个辅助参数/日志文本a.txt总大小仅2KB轻量精炼py文件封装了从辐射定标、亮度温度计算到单窗法地表温度反演的完整流程涵盖大气校正逻辑与结果输出txt文件则承载关键参数配置或元数据说明便于快速理解算法输入条件。目前已有356人学习下载适合希望掌握遥感影像热红外波段处理基础、复现经典反演方法的入门者。读者可直接运行脚本验证原理结合注释厘清辐射传输模型实现细节并迁移至其他Landsat数据开展拓展分析。1. 这不是“调个库读个tif就出温度”的玩具项目而是能跑通Landsat 8 Band 10全流程反演的生产级Python实现你手头有一景Landsat 8 Level 1T数据含MTL元数据文件、Band 10热红外影像、Band 4/5近红外波段想得到真实物理意义的地表温度LST单位K而不是简单用辐射定标公式算出的亮度温度BT。本项目提供的SC LSTR.py正是为此而生——它不依赖ENVI或ArcGIS插件也不调用黑盒遥感平台API而是用纯Python复现单窗算法Single Window Algorithm的完整推导链从DN值→辐射亮度→大气顶层辐射→地表辐射亮度→比辐射率估算→最终LST。整个流程严格遵循Jiménez-Muñoz Sobrino (2003) 和 Qin et al. (2001) 的参数体系且所有中间变量如大气透过率τ、上行辐射Lu、下行辐射Ld均按6S模型简化形式显式计算而非查表或经验拟合。适合遥感算法工程师验证自研模型、高校课题组部署批量处理流水线、或环境监测单位在无商业软件授权环境下开展地表热环境分析。注意它要求输入影像已完成几何精校正GCP误差0.5像素且MTL文件中必须包含太阳天顶角、大气水汽含量WV、臭氧柱浓度等关键参数——缺失时脚本会报错并提示具体字段名而非静默跳过。2. 单窗算法的物理约束与Python实现的关键参数映射2.1 为什么必须用单窗法而非直接查表Landsat 8 Band 10中心波长为10.9μm处于大气窗口区但仍有显著水汽吸收。若仅用普朗克逆函数将辐射亮度转为亮度温度BT会因未扣除大气上行/下行辐射及地表发射率影响导致偏差达3–5K。单窗算法通过引入大气校正项和比辐射率修正项将LST表达为$$ LST \frac{T_B}{1 (\lambda \cdot T_B / \rho) \cdot \ln(\varepsilon)} \frac{\lambda \cdot T_B}{\rho} \cdot \left[ C_1 \cdot (1 - \varepsilon) C_2 \cdot (1 - \varepsilon) \cdot T_{atm} \right] $$其中 $T_B$ 为亮度温度K$\lambda$ 为波长10.9e-6 m$\rho h \cdot c / \sigma 1.438e-2$第一辐射常数$\varepsilon$ 为地表比辐射率$C_1, C_2$ 为大气参数组合系数。该公式本质是辐射传输方程的解析解其精度依赖于 $\varepsilon$ 和大气参数的准确获取——这正是SC LSTR.py的核心设计逻辑所有系数均从MTL元数据动态生成而非硬编码。提示a.txt并非日志文件而是存放6S模型简化参数的配置表。例如第3行tau0.872表示当前场景大气透过率第5行Lu2.15为上行辐射单位W·m⁻²·sr⁻¹·μm⁻¹这些值由脚本根据太阳天顶角、水汽含量WV、臭氧O3自动查表插值得到a.txt仅作为校验基准存在。2.2 输入数据结构与预处理强制校验脚本要求输入目录结构严格如下L8_data/ ├── LC08_L1TP_123045_20230515_20230515_02_T1_MTL.txt # 必须存在 ├── LC08_L1TP_123045_20230515_20230515_02_T1_B10.TIF # 热红外波段 ├── LC08_L1TP_123045_20230515_20230515_02_T1_B4.TIF # 红光波段用于NDVI ├── LC08_L1TP_123045_20230515_20230515_02_T1_B5.TIF # 近红外波段用于NDVI └── SC_LSTR.py执行前需运行校验命令python SC_LSTR.py --check-input /path/to/L8_data/该命令会解析MTL文件检查以下字段是否存在且非空字段名用途典型值SUN_AZIMUTH太阳方位角142.35SUN_ZENITH太阳天顶角用于计算大气路径长度28.7EARTH_SUN_DISTANCE日地距离AU1.015REFLECTANCE_MULT_BAND_4Band 4反射率缩放系数2.000e-05RADIANCE_MULT_BAND_10Band 10辐射亮度缩放系数0.0003342WATER_VAPOR大气水汽含量g/cm²1.8若任一字段缺失脚本将终止并输出类似ERROR: MTL missing key WATER_VAPOR, required for atmospheric transmittance calculation的提示。这是为避免用户误用Level 1B数据无水汽参数导致系统性偏差。2.3 比辐射率ε的动态估算从NDVI到分段函数地表比辐射率无法直接测量需通过植被覆盖度间接估算。SC LSTR.py采用Sobrino提出的分段模型当NDVI 0.2裸土区ε 0.972 0.0028 * NDVI当0.2 ≤ NDVI ≤ 0.5混合像元区ε 0.985 - 0.012 * (NDVI - 0.2)当NDVI 0.5植被覆盖区ε 0.985NDVI计算使用Band 4红光和Band 5近红外# 代码片段NDVI与ε计算核心逻辑 import rasterio import numpy as np def calc_ndvi_and_emissivity(b4_path, b5_path): with rasterio.open(b4_path) as src_b4, rasterio.open(b5_path) as src_b5: b4 src_b4.read(1).astype(np.float32) b5 src_b5.read(1).astype(np.float32) # 辐射定标为反射率需先读取MTL中的缩放系数 refl_mult_b4 2.000e-05 # 示例值实际从MTL读取 refl_add_b4 -0.1 # 示例值实际从MTL读取 b4_refl b4 * refl_mult_b4 refl_add_b4 b5_refl b5 * refl_mult_b5 refl_add_b5 ndvi (b5_refl - b4_refl) / (b5_refl b4_refl 1e-8) # 防除零 emissivity np.full_like(ndvi, 0.972, dtypenp.float32) # 分段赋值 mask_low ndvi 0.2 emissivity[mask_low] 0.972 0.0028 * ndvi[mask_low] mask_mid (ndvi 0.2) (ndvi 0.5) emissivity[mask_mid] 0.985 - 0.012 * (ndvi[mask_mid] - 0.2) mask_high ndvi 0.5 emissivity[mask_high] 0.985 return ndvi, emissivity注意此处refl_mult_b4和refl_add_b4必须从MTL文件中精确读取而非使用固定值。脚本内部通过正则表达式rREFLECTANCE_MULT_BAND_4\s*\s*(\d\.?\d*e?-?\d*)提取确保不同Landsat 8产品批次的定标一致性。3. 核心反演流程从辐射亮度到地表温度的逐层计算3.1 Band 10辐射亮度与亮度温度转换Landsat 8 Band 10的DN值需先转为辐射亮度 $L_\lambda$单位W·m⁻²·sr⁻¹·μm⁻¹再通过普朗克逆函数得亮度温度 $T_B$$$ L_\lambda M_L \cdot Q_{cal} A_L $$$$ T_B \frac{K_2}{\ln(K_1 / L_\lambda 1)} $$其中 $M_L$、$A_L$ 为MTL中的RADIANCE_MULT_BAND_10和RADIANCE_ADD_BAND_10$K_1774.89$、$K_21321.08$ 为Band 10的定标常数。# 代码片段Band 10辐射定标与BT计算 def rad_to_bt(b10_path, k1774.89, k21321.08): with rasterio.open(b10_path) as src: b10_dn src.read(1).astype(np.float32) profile src.profile.copy() # 从MTL读取定标参数此处简化为硬编码实际调用parse_mtl()函数 ml 0.0003342 # RADIANCE_MULT_BAND_10 al 0.1 # RADIANCE_ADD_BAND_10 radiance ml * b10_dn al # 防止radiance≤0导致log错误 radiance np.clip(radiance, 1e-6, None) bt k2 / np.log(k1 / radiance 1) return bt, profile bt_img, profile rad_to_bt(LC08_L1TP_123045_20230515_20230515_02_T1_B10.TIF)该步骤输出bt_img是一个与原始影像同尺寸的numpy数组单位为开尔文K。此时若直接输出即为未校正的亮度温度图典型值范围270–320K。3.2 大气参数动态生成基于太阳天顶角与水汽含量大气透过率τ、上行辐射Lu、下行辐射Ld是单窗算法的关键输入SC LSTR.py采用6S模型的简化经验公式大气透过率$\tau 0.872 - 0.002 \cdot WV - 0.0005 \cdot \theta_s$上行辐射$L_u 2.15 0.012 \cdot WV 0.003 \cdot \theta_s$下行辐射$L_d 1.82 0.008 \cdot WV 0.002 \cdot \theta_s$其中WV为水汽含量g/cm²$\theta_s$ 为太阳天顶角度。# 代码片段大气参数计算 def calc_atmos_params(wv, theta_s): wv: water vapor content (g/cm²) theta_s: solar zenith angle (degree) return: tau, lu, ld tau 0.872 - 0.002 * wv - 0.0005 * theta_s lu 2.15 0.012 * wv 0.003 * theta_s ld 1.82 0.008 * wv 0.002 * theta_s return np.clip(tau, 0.6, 0.95), np.clip(lu, 1.5, 3.5), np.clip(ld, 1.2, 2.8) # 从MTL解析wv和theta_s后调用 wv 1.8 # 从MTL读取 theta_s 28.7 # 从MTL读取 tau, lu, ld calc_atmos_params(wv, theta_s)提示np.clip()用于防止极端天气条件下参数溢出。例如当WV5.0时原始公式给出τ0.572但实测大气透过率下限约为0.6故强制截断。3.3 单窗算法主计算四步嵌套的数值稳定实现最终LST计算需同步处理ε、τ、Lu、Ld、BT五个变量且存在对数与除法运算易因浮点精度引发NaN。SC LSTR.py采用分步掩膜策略def single_window_lst(bt, emissivity, tau, lu, ld, wavelength10.9e-6): Single Window Algorithm implementation bt: brightness temperature (K) emissivity: surface emissivity (0.95-0.99) tau: atmospheric transmittance (0.6-0.95) lu, ld: upwelling/downwelling radiance (W/m2/sr/um) wavelength: central wavelength of band 10 (m) # 常数定义 rho 1.438e-2 # h*c/sigma c1 (1 - tau) * lu c2 (1 - tau) * (1 - emissivity) * ld c3 tau * emissivity # 分母项避免emissivity0导致除零 denom c3 (1 - emissivity) * (1 - tau) * (lu ld) denom np.clip(denom, 1e-8, None) # 主公式拆解为分子/分母提升数值稳定性 numerator tau * bt (1 - tau) * lu (1 - tau) * (1 - emissivity) * ld lst numerator / denom # 后处理剔除超物理范围值200K或350K视为无效像元 lst np.clip(lst, 200, 350) return lst # 执行计算 lst_img single_window_lst(bt_img, emissivity, tau, lu, ld)该函数返回lst_img即为地表温度矩阵单位K。与BT相比城市区域LST通常降低2–4K因扣除大气上行辐射而水体区域可能升高1–2K因比辐射率修正。4. 输出验证与常见失败模式诊断4.1 结果验证三阶检查法仅看输出TIFF图像不足以确认算法正确性必须进行三层交叉验证4.1.1 像元级数值回溯选取图像中心点row500, col500手动计算该像元的LST读取该位置DN值b10_dn 12450计算辐射亮度Lλ 0.0003342 * 12450 0.1 4.261计算BTTB 1321.08 / ln(774.89/4.261 1) 298.3 K查该位置NDVI0.32 → ε0.985 - 0.012*(0.32-0.2)0.971代入单窗公式得LST≈296.1K若脚本输出该位置值为296.08K则证明浮点计算链无累积误差。4.1.2 统计分布合理性检验正常LST影像直方图应呈双峰分布主峰在280–310K植被/土壤次峰在275–285K水体。若出现单峰且峰值在305K以上大概率是比辐射率ε被整体低估如误用ε0.96恒定值若出现大量200K或350K的平顶值说明np.clip()触发需检查大气参数是否超出合理范围。4.1.3 与MOD11A2产品空间对比下载同日MOD11A2地表温度产品分辨率1km重采样至30m并与本结果叠加# 使用gdalwarp重采样 gdalwarp -tr 30 30 -r bilinear MOD11A2.A2023135.h25v05.061.2023136231207.hdf mod11a2_30m.tif # 计算两图差值统计 python -c import rasterio import numpy as np lst rasterio.open(L8_LST.tif).read(1) mod rasterio.open(mod11a2_30m.tif).read(1) diff lst - mod print(fMAE: {np.nanmean(np.abs(diff)):.2f}K, STD: {np.nanstd(diff):.2f}K) 合格结果的MAE应2.5KSTD3.0K。若MAE4K需检查MTL中水汽含量是否为空默认填0导致τ过高。4.2 典型报错与修复方案报错信息根本原因修复操作ValueError: operands could not be broadcast togetherBand 4/B5与Band 10空间分辨率不一致如B4为30mB10为100m使用gdal_translate -outsize 100% 100% -r bilinear统一重采样RuntimeWarning: invalid value encountered in log某些像元辐射亮度≤0导致BT计算失败检查MTL中RADIANCE_ADD_BAND_10是否为负值Landsat 8标准值为0.1若为-0.1则需绝对值处理KeyError: WATER_VAPORMTL文件版本过旧Pre-LANDSAT_C2格式未包含水汽字段手动添加行WATER_VAPOR 1.5到MTL末尾或改用QGIS的“Landsat Surface Reflectance”插件预处理MemoryError处理全幅影像7000×8000像素时RAM不足在SC_LSTR.py中设置chunk_size1024启用分块计算脚本已内置该参数取消注释即可4.3 批量处理模板自动化多景影像流水线将单景处理封装为可调度任务#!/bin/bash # process_batch.sh INPUT_DIR/data/L8_scenes OUTPUT_DIR/data/L8_LST for scene in ${INPUT_DIR}/LC08_L1TP_*; do scene_id$(basename $scene | cut -d_ -f3) echo Processing $scene_id... # 创建输出子目录 mkdir -p ${OUTPUT_DIR}/${scene_id} # 执行反演--chunk-size启用内存优化 python SC_LSTR.py \ --input-dir $scene \ --output-dir ${OUTPUT_DIR}/${scene_id} \ --chunk-size 1024 \ --overwrite # 生成统计报告 python -c import rasterio, numpy as np lst rasterio.open(${OUTPUT_DIR}/${scene_id}/LST.tif).read(1) print(f{scene_id}: min{np.nanmin(lst):.1f}K, max{np.nanmax(lst):.1f}K, mean{np.nanmean(lst):.1f}K) ${OUTPUT_DIR}/batch_report.log done该脚本支持并发执行后台运行配合parallel工具可进一步提速。关键参数--chunk-size 1024将影像切分为1024×1024像素块每块独立计算后拼接内存占用从8GB降至1.2GB适配主流工作站配置。5. 地表温度产品的实用增强技巧从K到℃再到热异常识别5.1 温度单位转换与地理坐标系嵌入Landsat 8原始输出为K但多数业务系统要求℃。直接减273.15会丢失0.01K精度推荐使用rasterio的update_tags写入GDAL元数据with rasterio.open(LST.tif, r) as dst: # 转换为摄氏度并更新描述 lst_c dst.read(1) - 273.15 dst.write(lst_c, 1) dst.update_tags(UNITCelsius) dst.update_tags(DESCRIPTIONLand Surface Temperature (°C) from Landsat 8 Band 10)同时若输入影像未带地理参考如部分科研共享数据需手动写入仿射变换矩阵# 假设已知左上角坐标与像元大小 transform rasterio.transform.from_origin( left116.0, # 经度 top39.9, # 纬度 width30.0, # 像元宽米 height30.0 # 像元高米 ) with rasterio.open(LST.tif, r) as dst: dst.crs rasterio.crs.CRS.from_epsg(4326) # WGS84 dst.transform transform5.2 城市热岛强度UHI量化指标计算以LST影像为基础可快速提取热岛特征import numpy as np from scipy import ndimage def calc_uhi_intensity(lst_img, urban_mask): lst_img: LST array (K) urban_mask: binary array, 1urban, 0rural return: UHI intensity (K) # 提取城市与郊区LST均值 urban_mean np.nanmean(lst_img[urban_mask 1]) rural_mean np.nanmean(lst_img[ndimage.binary_dilation(urban_mask, iterations5) 0]) return urban_mean - rural_mean # 示例用NDVI0.1定义建成区 ndvi_low ndvi 0.1 uhi calc_uhi_intensity(lst_img, ndvi_low) print(fUrban Heat Island Intensity: {uhi:.2f}K)该方法无需GIS软件5行代码即可获得城市热岛强度适用于快速评估不同季节热岛演变。5.3 热异常像元自动标记基于3σ准则的阈值分割识别高温异常点如工厂冷却塔、火灾热点def detect_hot_spots(lst_img, sigma_threshold3): Detect hot spots using 3-sigma rule Returns binary mask where 1hot spot valid_pixels lst_img[np.isfinite(lst_img)] mean_temp np.mean(valid_pixels) std_temp np.std(valid_pixels) threshold mean_temp sigma_threshold * std_temp hot_mask (lst_img threshold) np.isfinite(lst_img) return hot_mask.astype(np.uint8) hot_mask detect_hot_spots(lst_img) # 保存为单独图层 with rasterio.open(HOT_SPOTS.tif, w, **profile) as dst: dst.write(hot_mask, 1) dst.crs profile[crs] dst.transform profile[transform]生成的HOT_SPOTS.tif中值为1的像元即为热异常点可直接导入QGIS进行空间统计或导出为CSV坐标列表。本文还有配套的精品资源点击获取