### 使用 Python 实现 MODIS MYD11A2 数据单床算法热岛效应分析代码 MODIS MYD11A2 数据提供了地表温度（LST）信息，通常以 HDF 格式存储。为了实现基于该数据的单床算法以分析热岛效应，需要完成以下任务：读取 HDF 文件、提取 LST 数据、进行单位转换、处理云掩膜，并计算城市与周边区域的温度差异。 以下是实现这一目标的核心代码示例： ```python import os import numpy as np from pyhdf.SD import SD, SDC import rasterio from rasterio.transform import from_origin # 读取 MODIS MYD11A2 HDF 文件 def read_modis_hdf(file_path): hdf = SD(file_path, SDC.READ) datasets = hdf.datasets() lst_day_key = [key for key in datasets.keys() if 'LST_Day_1km' in key][0] lst_day = hdf.select(lst_day_key).get().astype(np.float32) # 提取日间地表温度 return lst_day # 单位转换：从科学计数法到摄氏度 def convert_to_celsius(lst_data): lst_celsius = np.where(lst_data > 0, lst_data * 0.02 - 273.15, np.nan) # 转换为摄氏度[^1] return lst_celsius # 处理云掩膜 def apply_quality_control(qc_data, lst_data): qc_flags = (qc_data >> 0) & 0b11 # 提取质量控制标志 valid_mask = (qc_flags == 0) | (qc_flags == 1) # 质量较好的像素 lst_data = np.where(valid_mask, lst_data, np.nan) # 应用掩膜 return lst_data # 计算热岛效应 def calculate_urban_heat_island(lst_data, urban_mask): urban_temps = lst_data[urban_mask] # 城市区域温度 rural_temps = lst_data[~urban_mask] # 非城市区域温度 uhi_intensity = np.nanmean(urban_temps) - np.nanmean(rural_temps) # 热岛强度 return uhi_intensity # 示例主函数 def main(modis_file, urban_mask_file, output_tif): # 读取 MODIS LST 数据 lst_data = read_modis_hdf(modis_file) lst_celsius = convert_to_celsius(lst_data) # 读取 QC 数据 qc_key = [key for key in SD(modis_file, SDC.READ).datasets().keys() if 'QC_Day' in key][0] qc_data = SD(modis_file, SDC.READ).select(qc_key).get() lst_celsius = apply_quality_control(qc_data, lst_celsius) # 加载城市掩膜 with rasterio.open(urban_mask_file) as src: urban_mask = src.read(1).astype(bool) # 计算热岛效应 uhi_intensity = calculate_urban_heat_island(lst_celsius, urban_mask) print(f"Urban Heat Island Intensity: {uhi_intensity:.2f}°C") # 输出结果为 GeoTIFF transform = from_origin(-180, 90, 0.01, 0.01) # 假设投影参数 new_dataset = rasterio.open( output_tif, 'w', driver='GTiff', height=lst_celsius.shape[0], width=lst_celsius.shape[1], count=1, dtype=rasterio.float32, crs='+proj=latlong', transform=transform, ) new_dataset.write(lst_celsius, 1) new_dataset.close() # 调用主函数 modis_file = "MYD11A2.A2023001.h23v05.006.2023002145047.hdf" urban_mask_file = "urban_mask.tif" output_tif = "lst_celsius.tif" main(modis_file, urban_mask_file, output_tif) ```