基于Landsat的1985-2024年地表温度LST分析，能否帮我写出关键步骤代码

<think>好的，我需要帮助用户解决如何使用Python和GDAL处理Landsat影像来分析1985-2024年的地表温度（LST）。首先，我需要回顾用户提供的引用内容，看看有没有相关的步骤或代码可以参考。 用户提到了引用[2]和[3]，其中引用[2]讨论了Landsat计算地表温度和TVDI的方法，但提到劈窗算法还不完善，适用于干旱区。引用[3]则详细介绍了如何用Python提取ECOSTRESS的地表温度数据，包括单位转换和数据处理。引用[4]提到了使用GDAL处理HDF格式数据，转换为TIFF并进行重投影，这可能对处理Landsat数据有用。 首先，处理Landsat数据的关键步骤可能包括数据下载、预处理（如辐射定标、大气校正）、地表温度反演（可能使用单窗算法或劈窗算法）、时间序列分析和结果输出。我需要确保每个步骤都有对应的代码示例，并且符合用户要求的GDAL和Python环境。 数据下载部分，用户可能需要从USGS Earth Explorer获取数据，这里可以建议使用API或手动下载。预处理部分，辐射定标和大气校正需要将DN值转换为辐射亮度，再转换为反射率或亮度温度。这里需要参考Landsat的元数据，比如ML和AL参数。 地表温度反演部分，引用[3]提到将开尔文转换为摄氏度，可能需要应用单窗算法，这需要地表比辐射率和大气参数。代码示例应包括读取热红外波段，应用辐射定标，转换为亮度温度，然后计算LST。 时间序列分析可能需要遍历多年的数据，计算年均温度或变化趋势，使用GDAL处理多期影像，并生成统计结果。 需要注意单位转换，比如引用[3]中的K转摄氏度，以及scale_factor的处理。此外，引用[4]中的HDF转TIFF和重投影可能适用于Landsat数据，特别是如果用户需要处理不同时期的影像，确保空间对齐。 在编写代码示例时，要使用GDAL读取和写入GeoTIFF文件，处理波段，进行数学运算，并确保投影一致。可能还需要处理不同的Landsat版本（如Landsat 5, 7, 8, 9），因为时间跨度从1985到2024年，不同卫星的波段编号和参数可能不同，需要分别处理。 最后，需要提醒用户注意不同算法的适用性，比如劈窗算法在干旱区的限制，以及单窗算法的假设条件。此外，质量控制波段的使用，如引用[4]中提到的，可能需要特别处理数据类型。 总结步骤后，我需要按照用户的要求，分步骤给出代码示例，并确保LaTeX格式正确，行内公式用$...$，独立公式用$$...$$，同时生成相关问题，并正确引用提供的参考资料。</think>以下是如何使用Python和GDAL处理Landsat数据实现地表温度(LST)反演的关键步骤及示例代码： --- ### 一、数据准备阶段 1. **数据下载** 从USGS Earth Explorer获取Landsat Level-2数据（已包含地表反射率产品）或Level-1数据（需自行辐射定标） 2. **数据预处理** ```python from osgeo import gdal import numpy as np # 读取热红外波段（以Landsat 8为例，B10为热红外波段） ds = gdal.Open('LC08_L1TP_123045_20220101_20220101_01_T1_B10.TIF') thermal_band = ds.GetRasterBand(1).ReadAsArray() ``` --- ### 二、地表温度反演流程 #### 1. 辐射定标 将DN值转换为辐射亮度（$L_{\lambda}$）： $$ L_{\lambda} = ML \times Q_{cal} + AL $$ ```python ML = 3.3420E-04 # 元数据中RADIANCE_MULT_BAND_10 AL = 0.10000 # 元数据中RADIANCE_ADD_BAND_10 radiance = ML * thermal_band + AL # [W/(m²·sr·μm)] ``` #### 2. 亮度温度转换 $$ T_{B} = \frac{K_{2}}{\ln(\frac{K_{1}}{L_{\lambda}} + 1)} $$ ```python K1 = 774.89 # 元数据中K1_CONSTANT_BAND_10 K2 = 1321.08 # 元数据中K2_CONSTANT_BAND_10 bt = K2 / np.log(K1 / radiance + 1) # 单位：开尔文(K) ``` #### 3. 地表温度计算（单窗算法[^2]） $$ LST = \frac{T_{B}}{1 + (\lambda \cdot T_{B}/\rho) \cdot \ln \epsilon} $$ ```python # 参数设置（需根据具体影像调整） wavelength = 10.9e-6 # Landsat 8 B10波长(μm) rho = 1.438e-2 # 常数 epsilon = 0.96 # 地表比辐射率 lst_kelvin = bt / (1 + (wavelength * bt / rho) * np.log(epsilon)) lst_celsius = lst_kelvin - 273.15 # 转换为摄氏度[^3] ``` --- ### 三、时间序列处理 ```python # 遍历多时相数据（示例框架） years = range(1985, 2025) for year in years: input_path = f'Landsat_{year}.tif' output_path = f'LST_{year}.tif' # 此处插入上述温度反演代码 # 结果保存 driver = gdal.GetDriverByName('GTiff') ds_out = driver.Create(output_path, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) ds_out.GetRasterBand(1).WriteArray(lst_celsius) ds_out.SetGeoTransform(ds.GetGeoTransform()) ds_out.SetProjection(ds.GetProjection()) ds_out = None ``` --- ### 四、结果验证与可视化 ```python # 温度统计 print(f"Max LST: {np.nanmax(lst_celsius):.2f}°C") print(f"Min LST: {np.nanmin(lst_celsius):.2f}°C") # 生成温度分级图 import matplotlib.pyplot as plt plt.imshow(lst_celsius, cmap='jet', vmin=0, vmax=50) plt.colorbar(label='Temperature (°C)') plt.show() ``` --- ### 关键注意事项 1. Landsat 5/7/8/9的热红外波段编号不同（如Landsat 5为B6） 2. 地表比辐射率$\epsilon$需通过NDVI计算[^2] 3. 建议使用Level-2数据可跳过辐射定标步骤 4. 长时间序列分析需统一空间分辨率（建议保持30m）