【Python实例】读取/处理 Landsat LST数据

下载Landsat LST数据的详细过程可参见另一博客-【数据集】EarthExplore下载Landsat LST 数据。本博客提供详细Python代码进行LST数据的可视化以及处理。

数据准备

以 Landsat 8 的地表温度波段（ST_B10.TIF）为例，实现以下两个功能：
一、读取并可视化地表温度（ST_B10）图像+云层覆盖信息（QA）图像
二、提取影像详细信息（包括元数据、云层、空间信息等）

地表温度单位说明：
Landsat 8 ST_B10 是以 数字数值（DN） 存储的，转换为 开尔文（K） 的公式为：

ST = DN * 0.00341802 + 149.0

💡 QA 掩码如何工作：

Landsat 地表温度 QA 波段是一个 位掩码（bitmask），每个像素的值是一个整数，具体位代表不同特征：

位段	含义
0–1	填充值标记（无效像素）
2–3	云层
4–5	云影
6–7	雪/冰

基于Python读取Landsat LST数据

1_LST信息读取及可视化

绘制的地块 LST 数据如下：
此图像出现了一条 垂直白色带状区域，表示该区域的像素值为 缺失值或无效值（NaN），所以在绘图时显示为白色。

🔍 可能原因分析（按优先级排序）

编号	原因	说明
🟡 1	卫星扫描缝隙 / 数据缺失（缺口）	Landsat 卫星扫描为“摆扫式”，在某些路径行可能会出现写入失败、断带或缺失
🟠 2	云层遮挡被掩膜掉了	如果 QA 掩膜将该区域标记为云或云影，温度图就被设置为 NaN
🔵 3	卫星姿态调整或角度问题	有些场景在极端角度下成像，边缘存在畸变或数据空白
⚫ 4	数据产品问题 / 处理失败	处理过程中某波段损坏、投影失败，也可能导致区域为空

完整Python实现代码如下：

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from rasterio.plot import show
from rasterio.plot import show_hist
from rasterio.plot import plotting_extent
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.patches import Rectangle
import os# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'# === 1. 设置路径 ===
lst_path = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\LC08_L2SP_127039_20190117_20200829_02_T1_ST_B10.TIF"
qa_path = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\LC08_L2SP_127039_20190117_20200829_02_T1_ST_QA.TIF"save_dir = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\Figures"if not os.path.exists(save_dir):os.makedirs(save_dir)# === 读取 LST 波段 ===
with rasterio.open(lst_path) as src:band = src.read(1)metadata = src.metabounds = src.boundscrs = src.crstransform = src.transformnodata = src.nodataextent = plotting_extent(src)# 打印元信息
print("📄 File metadata:")
print(f" - Data type: {metadata['dtype']}")
print(f" - Resolution: {src.res}")
print(f" - CRS: {crs}")
print(f" - Bounds: {bounds}")
print(f" - NoData value: {nodata}")
print("\n🌡 Pixel value range (DN):")
print(f" - Min: {np.nanmin(band)}")
print(f" - Max: {np.nanmax(band)}")# 转换为温度（摄氏度）
lst_kelvin = band * 0.00341802 + 149.0
lst_kelvin = np.where(band == nodata, np.nan, lst_kelvin)
lst_celsius = lst_kelvin - 273.15# === 3. 可视化地表温度 + 经纬度 ===
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(lst_celsius, cmap='plasma', extent=extent)
ax.set_title("Landsat 8 Surface Temperature (°C)", fontsize=15)
ax.set_xlabel("Longitude", fontsize=14)
ax.set_ylabel("Latitude", fontsize=14)
ax.grid(True, linestyle='--', linewidth=0.5)# ✅ 设置 colorbar 与主图高度相同
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label("Temperature (°C)", fontsize=14)# 保存图像
output_path_lst = os.path.join(save_dir, "Landsat8_Surface_Temperature.png")
plt.tight_layout()
plt.savefig(output_path_lst, dpi=300)
plt.show()

输出信息如下：

📄 File metadata:- Data type: uint16- Resolution: (30.0, 30.0)- CRS: EPSG:32648- Bounds: BoundingBox(left=642885.0, bottom=3239685.0, right=869415.0, top=3470715.0)- NoData value: 0.0🌡 Pixel value range (DN):- Min: 0- Max: 47143

2_QA信息读取及可视化

这幅图是根据 TIF 文件 ST_QA.TIF 中的 “位掩码（bitmask）” 生成的，所显示的是：

区域	含义	在图中颜色
黑色区域	被标记为云层或云影的像素	黑色（值为 True）
白色区域	被标记为清晰无云的像素	白色（值为 False）
灰色区域（若有）	可能是图像压缩造成混色或边缘过渡	视觉效果

完整Python实现代码如下：

# === 4. 读取 QA 波段并解析云掩膜 ===
with rasterio.open(qa_path) as qa_src:qa_band = qa_src.read(1)qa_extent = plotting_extent(qa_src)# Landsat 8 ST_QA 位掩码含义（示例：常用位段）
# Bit 0-1: Designated Fill
# Bit 2-3: Cloud
# Bit 4-5: Cloud shadow
# Bit 6-7: Snow/Ice
# Bit 8-9+: Cirrus, etc.# 提取云层 & 云影掩膜（以下是简化示例，具体可根据官方文档调整）
cloud_mask = np.bitwise_and(qa_band, 0b00001100) > 0
shadow_mask = np.bitwise_and(qa_band, 0b00110000) > 0# 创建综合掩膜
cloud_and_shadow = np.logical_or(cloud_mask, shadow_mask)# 可视化掩膜
plt.figure(figsize=(8, 6))
plt.imshow(cloud_and_shadow, cmap='gray', extent=extent)
plt.title("Cloud and Shadow Mask (QA Band)", fontsize=15)
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)
plt.grid(True, linestyle='--', linewidth=0.5)"""
# ✅ 添加图像边框（外框）
rect = Rectangle((qa_extent[0], qa_extent[2]),  # 左下角 (minX, minY)qa_extent[1] - qa_extent[0],   # 宽度 (maxX - minX)qa_extent[3] - qa_extent[2],   # 高度 (maxY - minY)linewidth=1.5,edgecolor='black',facecolor='none'
)
ax.add_patch(rect)
"""# ✅ 方法：手动用 plot 画边框
x_min, x_max, y_min, y_max = qa_extent[0], qa_extent[1], qa_extent[2], qa_extent[3]# 画4条边线
ax.plot([x_min, x_max], [y_min, y_min], color='black', linewidth=1.5)  # 下边
ax.plot([x_min, x_max], [y_max, y_max], color='black', linewidth=1.5)  # 上边
ax.plot([x_min, x_min], [y_min, y_max], color='black', linewidth=1.5)  # 左边
ax.plot([x_max, x_max], [y_min, y_max], color='black', linewidth=1.5)  # 右边# 保存图像
output_path_qa = os.path.join(save_dir, "Landsat8_QA_Cloud_Mask.png")
plt.tight_layout()
plt.savefig(output_path_qa, dpi=300)
plt.show()# 输出云像素统计
total_pixels = cloud_and_shadow.size
cloud_pixels = np.sum(cloud_and_shadow)
print(f"\n☁️ Cloud/Shadow Pixels: {cloud_pixels} ({cloud_pixels / total_pixels * 100:.2f}%)")print(f"✅ QA mask image saved to: {output_path_qa}")
print(f"✅ LST image saved to: {output_path_lst}")

输出信息如下：

☁️ Cloud/Shadow Pixels: 55351269 (95.19%)

3_基于云层/阴影掩膜图层 (ST_QA.TIF)得到去云温度图

目标任务

1. 使用 cloud_and_shadow 掩膜，将受影响像素设为 NaN
2. 绘制“去云”后的地表温度图（单位：°C）
3. 导出处理后的图像为 .png 图像（可视化图）和 .tif（GeoTIFF）

绘制的去云温度图如下：

完整实现代码如下（Python）：

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from rasterio.plot import plotting_extent
from mpl_toolkits.axes_grid1 import make_axes_locatable
from rasterio.transform import Affine
from rasterio.crs import CRS
import os
import numpy.ma as ma# 设置字体
plt.rcParams['font.family'] = 'Times New Roman'# === 文件路径 ===
lst_path = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\LC08_L2SP_127039_20190117_20200829_02_T1_ST_B10.TIF"
qa_path = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\LC08_L2SP_127039_20190117_20200829_02_T1_ST_QA.TIF"
save_dir = r"D:\0 DataBase\13_EarthExplore_LC08_L2SP\Figures"if not os.path.exists(save_dir):os.makedirs(save_dir)# === 读取 LST 波段 ===
with rasterio.open(lst_path) as src:lst_dn = src.read(1)lst_meta = src.meta.copy()lst_extent = plotting_extent(src)lst_crs = src.crslst_transform = src.transformnodata = src.nodata# 转换为摄氏度
lst_kelvin = lst_dn * 0.00341802 + 149.0
lst_kelvin = np.where(lst_dn == nodata, np.nan, lst_kelvin)
lst_celsius = lst_kelvin - 273.15# === 读取 QA 波段并创建掩膜 ===
with rasterio.open(qa_path) as qa_src:qa_band = qa_src.read(1)cloud_mask = np.bitwise_and(qa_band, 0b00001100) > 0
shadow_mask = np.bitwise_and(qa_band, 0b00110000) > 0
cloud_and_shadow = np.logical_or(cloud_mask, shadow_mask)# === 应用掩膜：将云层和阴影区域设为 NaN ===
lst_masked = np.where(cloud_and_shadow, np.nan, lst_celsius)# 使用掩膜数组隐藏 NaN
masked_array = ma.masked_invalid(lst_masked)# === 可视化去云后的地表温度图 ===
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(masked_array, cmap='plasma', extent=lst_extent)
ax.set_title("Cloud-masked Surface Temperature (°C)", fontsize=15)
ax.set_xlabel("Longitude", fontsize=14)
ax.set_ylabel("Latitude", fontsize=14)
ax.grid(True, linestyle='--', linewidth=0.5)from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label("Temperature (°C)")plt.tight_layout()
plt.savefig(os.path.join(save_dir, "LST_Masked_by_QA_fixed.png"), dpi=300)
plt.show()png_path = os.path.join(save_dir, "LST_Masked_by_QA.png")
print(f"✅ Cloud-masked temperature PNG saved to:\n{png_path}")# === 导出 GeoTIFF（去云地表温度）===
tif_path = os.path.join(save_dir, "LST_Masked_by_QA.tif")# 更新元数据
lst_meta.update({"dtype": "float32","nodata": np.nan,"count": 1,"crs": lst_crs,"transform": lst_transform
})with rasterio.open(tif_path, "w", **lst_meta) as dst:dst.write(lst_masked.astype(np.float32), 1)print(f"✅ Cloud-masked temperature GeoTIFF saved to:\n{tif_path}")

参考