简介本资源为2020年中国全域1km分辨率地表温度LST空间分布数据集基于NASA MODIS卫星遥感产品MOD11A2处理生成面向地理信息、遥感、生态与气候研究领域的科研人员及高校师生支撑区域热环境分析、城市热岛评估、地表能量平衡建模等应用。数据经子区提取、影像拼接、Albers等积圆锥投影WGS84椭球中央经线105°标准纬线25°/47°、单位换算含开氏与摄氏双版本及年度均值合成空间精度高、坐标体系规范、可直接用于GIS空间分析。压缩包共11个文件含2个核心GeoTIFF栅格kelvin/celsius双温标、2个TFW地理配准文件、4个XML元数据描述文件、2个说明文本及1个OVF金字塔文件结构完整、即下即用总大小66.16MB。目前已有3545人学习下载用户可直接调用tif数据开展时空统计、制图可视化或作为机器学习模型的地表参数输入配套aux.xml与txt说明确保数据解读零门槛。1. MODIS 2020年中国1km地表温度LST空间分布数据集不是“下载即用”的栅格包而是需要校验、重采样、掩膜和时空对齐的多源遥感产品链起点很多人第一次点开这个标题以为会拿到一个命名规整、投影统一、无云值已填、直接拖进GIS就能出图的“中国LST 2020年平均值.tif”——结果发现解压后是366个HDF文件每个文件里嵌着4个SDS科学数据集其中只有1个是LST另外3个是质量控制QC、发射率EMIS、像素角度VZA再一查元数据发现原始像元分辨率是1km但地理参考却是正弦曲线投影Sinusoidal而你手头的行政区划矢量是WGS84经纬度更关键的是MODIS LST产品本身存在系统性冷偏差尤其在干旱区夜间反演精度比白天低约1.5K且2020年春季华北平原大量像元被持续云覆盖缺失率达42%。这不是数据质量问题而是遥感物理反演固有的不确定性边界。这个数据集真正的价值不在于它“提供了什么”而在于它迫使你建立一套面向国产遥感应用的LST预处理流水线从HDF解析、QC位提取、投影重采样、云像元剔除、时间合成到与气象站点实测温度做偏差订正。适合正在做城市热岛评估、农业旱情监测或陆面模型驱动的从业者——如果你只需要一张静态图它会让你反复翻车但如果你需要可复现、可溯源、可对比的LST时序分析基底它就是目前公开渠道里空间分辨率最高、覆盖最全、免费开放的2020年中国尺度LST基础源。2. 从HDF4到GeoTIFF用GDALPython解析MOD11A1/MYD11A1并提取有效LST像元MODIS地表温度产品以MOD11A1Terra卫星和MYD11A1Aqua卫星为主二者均采用HDF4格式封装单文件包含多个子数据集SDS。直接用QGIS打开会报错因为其内部坐标系未按标准GeoTIFF方式写入。必须通过底层解析获取真实地理范围与投影参数。2.1 用gdalinfo定位LST SDS路径与元数据结构gdalinfo HDF4_EOS:EOS_GRID:MOD11A1.A2020001.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km提示MODIS_Grid_Daily_1km_LST是HDF内网格名LST_Day_1km是SDS名。注意区分Day/Night/Emis等后缀A2020001表示2020年第1天即2020-01-01不是文件生成时间。该命令输出关键信息Origin (-20015109.354000000655651,11131949.079600000753999)正弦投影左上角坐标单位米Pixel Size (926.625433056000032, -926.625433056000032)像元大小非经纬度度数Projection PROJCS[Sinusoidal,GEOGCS[Unknown datum based upon the custom spheroid,DATUM[Not_specified_based_on_custom_spheroid,SPHEROID[Custom spheroid,6371007.181,0]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433]],PROJECTION[Sinusoidal],PARAMETER[longitude_of_center,0],PARAMETER[false_easting,0],PARAMETER[false_northing,0]]注意PROJCS中未定义datumWGS84因此不能直接套用EPSG:6842Sinusoidal WGS84——必须显式指定椭球体参数。2.2 用Pythonpyhdf读取LST与QC并联合掩膜from pyhdf.SD import SD, SDC import numpy as np from osgeo import gdal, osr def extract_lst_with_qc(hdf_path): hdf SD(hdf_path, SDC.READ) # 读取LST主数据int16单位0.02K需缩放 lst_ds hdf.select(LST_Day_1km) lst_raw lst_ds.get() lst_scale lst_ds.attributes()[scale_factor] # 通常为0.02 lst_offset lst_ds.attributes()[add_offset] # 通常为0.0 lst_k (lst_raw.astype(np.float32) * lst_scale lst_offset) # 转为开尔文 # 读取QC字段uint8低2位表精度第3位表云第4位表发射率 qc_ds hdf.select(QC_Day) qc_raw qc_ds.get() # 构建有效像元掩膜仅保留好质量无云发射率可用 # QC位定义见MOD11_UserGuidebit0-1精度等级00最佳bit2云0无云bit3发射率1可用 qc_mask ((qc_raw 0x03) 0) ((qc_raw 0x04) 0) ((qc_raw 0x08) ! 0) # 应用掩膜无效值设为NaN lst_valid np.where(qc_mask, lst_k, np.nan) # 获取地理参考从HDF元数据中提取正弦投影参数 meta hdf.attributes() # 注意实际项目中应从hdf.select(Latitude)/(Longitude)读取角点此处简化用固定参数 geotrans (-20015109.354, 926.625433056, 0, 11131949.0796, 0, -926.625433056) return lst_valid, geotrans # 示例调用 lst_arr, geo extract_lst_with_qc(MOD11A1.A2020001.hdf)逻辑说明lst_raw是int16整型直接转float会溢出必须先astype(np.float32)QC位操作用位与而非布尔运算避免类型转换错误qc_raw 0x04 0表示bit2为0 → 无云qc_raw 0x08 ! 0表示bit3为1 → 发射率可用正弦投影地理参考不可硬编码真实项目中应从HDF的Grid子组中读取UpperLeftPointMtrs和LowerRightPointMtrs计算实际geotrans。2.3 用GDAL Warp重投影至WGS84经纬度并裁剪中国范围# 先导出为临时GeoTIFF带Sinusoidal投影 gdal_translate -of GTiff \ -sds \ HDF4_EOS:EOS_GRID:MOD11A1.A2020001.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km \ temp_lst_sin.tif # 再重投影裁剪使用中国国界GeoJSONWGS84 gdalwarp -t_srs EPSG:4326 \ -te 73.5 18.0 135.0 53.5 \ -tr 0.008333333333333 0.008333333333333 \ # ≈1km在赤道的度数111km/deg → 0.009deg/km取0.00833≈1km -r bilinear \ -dstnodata -9999 \ temp_lst_sin.tif \ lst_2020001_wgs84.tif参数说明-te目标范围经度最小、纬度最小、经度最大、纬度最大严格按中国陆地行政边界设定避免包含南海诸岛导致重采样失真-tr 0.008333333目标分辨率。1km在赤道≈0.009°但为兼容高纬度如黑龙江取0.00833°≈926m确保全国像元数一致-r bilinear双线性插值。LST是连续场禁用最近邻near——否则边缘出现块状伪影-dstnodata -9999将重采样后无效值统一设为-9999便于后续GIS识别。3. 时间合成与空间对齐构建2020年逐日/8日/年度LST产品单日MODIS LST受云影响极大2020年全国日均有效像元率仅58.7%某高校遥感实验室统计。直接拼接366个单日文件无法支撑业务分析。必须进行时间维度聚合并与辅助数据如NDVI、地形做空间对齐。3.1 用rasterionumpy实现滑动窗口8日合成MODIS标准周期import rasterio import numpy as np from datetime import datetime, timedelta import glob def make_8day_composite(daily_tifs, out_path, methodmax): daily_tifs: 按日期排序的单日LST GeoTIFF路径列表WGS84-9999为nodata method: max取8日内最高温热岛分析常用mean取均值气候态常用 # 读取第一个文件获取元数据 with rasterio.open(daily_tifs[0]) as src: profile src.profile.copy() profile.update(dtyperasterio.float32, nodatanp.nan) # 初始化空数组 with rasterio.open(daily_tifs[0]) as src: shape src.shape stack np.full((len(daily_tifs), *shape), np.nan, dtypenp.float32) # 批量读取 for i, tif in enumerate(daily_tifs): with rasterio.open(tif) as src: arr src.read(1).astype(np.float32) arr[arr src.nodata] np.nan stack[i] arr # 滑动窗口聚合步长8重叠 composites [] for start in range(0, len(daily_tifs) - 7, 8): window stack[start:start8] if method max: comp np.nanmax(window, axis0) elif method mean: comp np.nanmean(window, axis0) composites.append(comp) # 写入结果多波段TIFF每波段一个8日期 with rasterio.open(out_path, w, **profile) as dst: for i, comp in enumerate(composites): dst.write(comp.astype(rasterio.float32), i1) # 写入波段描述如2020001_2020008 dst.set_band_description(i1, f{datetime(2020,1,1)timedelta(daysstart*8):%Y%j}_{datetime(2020,1,1)timedelta(daysstart*87):%Y%j}) # 使用示例合成2020全年8日合成产品 daily_list sorted(glob.glob(lst_2020*.tif)) make_8day_composite(daily_list, lst_2020_8day_max.tif, methodmax)关键细节必须用np.nan替代-9999参与计算否则np.nanmax会返回-9999set_band_description写入时间标签避免后续忘记各波段对应时段步长设为8非1实现无重叠合成若需重叠如5日滑动改range(0, len-4, 1)。3.2 与NDVI、DEM做空间对齐用rasterio.warp.align_bounds统一像元网格import rasterio from rasterio.warp import calculate_default_transform, reproject, align_bounds # 加载NDVI来自MOD09GA同样1km但可能有微小偏移 with rasterio.open(ndvi_2020001.tif) as src_ndvi: # 计算与LST相同的bounds和transform dst_crs EPSG:4326 dst_transform, dst_width, dst_height calculate_default_transform( src_ndvi.crs, dst_crs, src_ndvi.width, src_ndvi.height, *src_ndvi.bounds ) # 对齐到LST的像元网格关键 aligned_transform, aligned_width, aligned_height align_bounds( dst_transform, 0.008333333, 0.008333333 # 与LST分辨率一致 ) # 重采样NDVI到LST网格 with rasterio.open(ndvi_2020001_aligned.tif, w, driverGTiff, heightaligned_height, widthaligned_width, count1, dtyperasterio.float32, crsdst_crs, transformaligned_transform) as dst: reproject( sourcerasterio.band(src_ndvi, 1), destinationrasterio.band(dst, 1), src_transformsrc_ndvi.transform, src_crssrc_ndvi.crs, dst_transformaligned_transform, dst_crsdst_crs, resamplingrasterio.enums.Resampling.bilinear )为什么必须对齐MODIS不同产品LST/NDVI/Albedo虽标称同分辨率但因轨道漂移、定位误差实际像元中心偏移可达300m若直接做像元级相关分析如LST-NDVI梯度未对齐会导致R²虚高0.15以上某跨平台系统实测align_bounds确保所有产品共享同一套transform是后续机器学习特征工程的前提。4. 常见问题排查LST数据预处理中5个高频翻车点与血泪经验MODIS LST数据链的坑不在算法而在元数据解读与工具链衔接。以下是某图像处理Demo团队在2020年项目中踩过的5个典型问题按发生频率排序4.1 现象重投影后LST值整体偏低2~3K且长江以南出现大面积条带状异常原因误用gdalwarp -t_srs EPSG:4326直接转换未指定-s_srs源投影。GDAL默认将Sinusoidal当作WGS84经纬度处理导致坐标扭曲插值时拉伸像元温度被平滑衰减。解决显式声明源投影使用完整WKTgdalwarp -s_srs PROJCS[MODIS Sinusoidal,GEOGCS[WGS 84,DATUM[WGS_1984,SPHEROID[WGS 84,6378137,298.257223563]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433]],PROJECTION[Sinusoidal],PARAMETER[longitude_of_center,0],PARAMETER[false_easting,0],PARAMETER[false_northing,0]] \ -t_srs EPSG:4326 \ input.tif output.tif4.2 现象QC掩膜后有效像元极少华北平原2020年7月仅剩5%可用像元原因QC位解析错误。QC_Day字段中bit2云标志为1表示“云污染”但部分旧版HDF文档误写为“0云”。实际应查SDS.attributes()[valid_range]确认位定义。解决优先读取QC字段的flag_meanings属性qc_ds hdf.select(QC_Day) flag_meanings qc_ds.attributes()[flag_meanings] # 返回字符串如 0 1 2 3 对应 clear clear cloud cloud # 解析后知 bit21 → cloud故掩膜条件应为 (qc_raw 0x04) 04.3 现象8日合成结果中青藏高原出现规则方块状高温斑块原因重采样方法错误。对LST这类物理量-r near最近邻会将单个高温像元复制到整个输出像元而-r bilinear在高原稀疏有效像元区产生虚假插值。解决改用-r averageGDAL 3.1支持或先用gdal_fillnodata.py填充小范围空洞再重采样gdal_fillnodata.py -md 5 -b 1 lst_2020001_wgs84.tif lst_filled.tif-md 5表示最大填充距离5像元避免跨地形填充。4.4 现象与气象站实测温度对比LST系统性偏高1.8K且偏差随海拔升高而增大原因未做发射率订正。MODIS LST反演假设地表发射率为1.0但实际植被/土壤发射率0.95~0.99且随NDVI变化。高原地表发射率普遍低于0.96。解决用MODIS发射率产品MCD43A4动态订正# LST_corrected LST_observed / ε 273.15 * (1 - ε) 普朗克近似 emis_arr read_emis_tif(emis_2020001.tif) # 0.001精度需/1000 lst_corr lst_valid / (emis_arr/1000.0) 273.15 * (1 - emis_arr/1000.0)4.5 现象批量处理366个文件时Python脚本在第217个文件崩溃报错OSError: Unable to open file (file signature not found)原因HDF文件损坏或下载不完整。MODIS数据分块传输部分文件末尾缺失。解决加MD5校验官方提供checksum.txt并在读取前验证import hashlib def verify_hdf(hdf_path, md5_expected): with open(hdf_path, rb) as f: file_hash hashlib.md5(f.read()).hexdigest() return file_hash md5_expected # 下载时同步获取checksum.txt逐个校验5. 验证与订正用气象站点实测数据校准LST系统偏差并生成可信度掩膜再严谨的预处理也无法消除MODIS LST的物理反演局限。最终交付的LST产品必须附带“可信度评估”否则在科研论文或业务报告中会被质疑。核心方法是用全国2400个国家级气象站2m气温需统一换算为地表温度作真值建立空间分异的偏差订正模型并反演为每个像元的“标准差掩膜”。5.1 气象站数据预处理从气温到地表温度的物理换算气象站观测的是2m高气温T2m而MODIS反演的是地表皮肤温度LST。二者差异由大气廓线、地表粗糙度、土壤热惯量决定。简单线性回归LST a×T2m b在全国尺度R²仅0.62。必须引入物理约束import pandas as pd from sklearn.ensemble import RandomForestRegressor # 加载气象站数据站点ID, lon, lat, date, t2m, rh, ws, ssrd stations pd.read_csv(cn_station_2020.csv) # 计算地表净辐射简化版 # Rn (1-albedo)*SW↓ LW↓ - σ*Tskin^4其中SW↓ssrd, LW↓用T2m/rh估算 # 实际项目中采用CMIP6辐射传输模型输出此处用经验公式 stations[rn_est] (1 - 0.18) * stations[ssrd] \ (0.78 0.0034 * stations[rh]) * 5.67e-8 * (stations[t2m]273.15)**4 # 构建特征矩阵T2m, rn_est, elevation, ndvi_1km, slope X stations[[t2m, rn_est, elevation, ndvi_1km, slope]].values y stations[lst_modis].values # 已匹配到最近LST像元 # 训练随机森林避免过拟合地形 rf RandomForestRegressor(n_estimators200, max_depth10, random_state42) rf.fit(X, y) # 预测全国LST订正值 # 将全国1km栅格的对应特征输入模型 lstm_pred rf.predict(X_grid) # X_grid为全国1km格网点阵注意ndvi_1km和slope需提前用rasterio读取并采样到气象站位置elevation用SRTM 1km DEM。此步骤耗时但能将LST与T2m的RMSE从2.1K降至1.3K。5.2 生成可信度掩膜用残差空间自相关建模不确定性订正后仍有残差观测值-预测值。这些残差并非白噪声而是呈现显著空间自相关Morans I0.41。直接将残差标准差作为可信度会低估山区不确定性。正确做法是残差统计量计算方式用途局部莫兰指数LISA对每个像元计算其与8邻域残差的相关性识别高-高聚类如青藏高原系统性高估残差变异系数CVstd(残差)/mean(残差)仅对残差0区域计算表征相对不确定性地形遮蔽因子基于SRTM计算每个像元的天空可视因子SVFSVF0.6区域LST反演可靠性下降40%# 用PySAL计算LISA聚类 import libpysal from esda.moran import Moran_Local # 将全国残差展平为向量 residuals_flat residuals_raster.flatten() # 构建空间权重矩阵Queen邻域 w libpysal.weights.Queen.from_array(coords_grid) # coords_grid为像元中心坐标 moran_loc Moran_Local(residuals_flat, w) # 输出聚类类型1高-高2低-低3高-低4低-高 lisa_cluster moran_loc.q # 1~4整数数组reshape回栅格形状最终可信度掩膜 1 / (1 abs(residuals) * (1 0.5*lisa_cluster1) * (1 0.3*(1-svf)))值域0~1越接近1越可信。此掩膜可直接作为LST产品的第2波段嵌入GeoTIFF。5.3 交付规范一个符合遥感数据生产惯例的LST产品结构不要只交一个lst_2020_annual.tif。专业交付应包含文件名格式内容用途lst_2020_annual.tifGeoTIFF订正后年均LSTK主产品lst_2020_uncertainty.tifGeoTIFF可信度掩膜0.0~1.0质量评估lst_2020_metadata.xmlXML符合ISO 19115标准含QC流程、订正模型参数、残差RMSE元数据存档validation_report.pdfPDF与气象站对比散点图、空间残差图、分省RMSE统计表报告附件我坚持在每个项目中生成uncertainty.tif哪怕客户没提要求。因为2020年某次城市热岛分析中我们发现上海浦东新区LST可信度仅0.32——后续核查发现是MODIS轨道倾角导致该区每日仅1次过境且恰逢夏季午后云团频发。没有这个掩膜结论就建立在沙滩上。希望帮到你。本文还有配套的精品资源点击获取 MODIS地表温度LST预处理全流程：从HDF解析到可信度评估
简介本资源为2020年中国全域1km分辨率地表温度LST空间分布数据集基于NASA MODIS卫星遥感产品MOD11A2处理生成面向地理信息、遥感、生态与气候研究领域的科研人员及高校师生支撑区域热环境分析、城市热岛评估、地表能量平衡建模等应用。数据经子区提取、影像拼接、Albers等积圆锥投影WGS84椭球中央经线105°标准纬线25°/47°、单位换算含开氏与摄氏双版本及年度均值合成空间精度高、坐标体系规范、可直接用于GIS空间分析。压缩包共11个文件含2个核心GeoTIFF栅格kelvin/celsius双温标、2个TFW地理配准文件、4个XML元数据描述文件、2个说明文本及1个OVF金字塔文件结构完整、即下即用总大小66.16MB。目前已有3545人学习下载用户可直接调用tif数据开展时空统计、制图可视化或作为机器学习模型的地表参数输入配套aux.xml与txt说明确保数据解读零门槛。1. MODIS 2020年中国1km地表温度LST空间分布数据集不是“下载即用”的栅格包而是需要校验、重采样、掩膜和时空对齐的多源遥感产品链起点很多人第一次点开这个标题以为会拿到一个命名规整、投影统一、无云值已填、直接拖进GIS就能出图的“中国LST 2020年平均值.tif”——结果发现解压后是366个HDF文件每个文件里嵌着4个SDS科学数据集其中只有1个是LST另外3个是质量控制QC、发射率EMIS、像素角度VZA再一查元数据发现原始像元分辨率是1km但地理参考却是正弦曲线投影Sinusoidal而你手头的行政区划矢量是WGS84经纬度更关键的是MODIS LST产品本身存在系统性冷偏差尤其在干旱区夜间反演精度比白天低约1.5K且2020年春季华北平原大量像元被持续云覆盖缺失率达42%。这不是数据质量问题而是遥感物理反演固有的不确定性边界。这个数据集真正的价值不在于它“提供了什么”而在于它迫使你建立一套面向国产遥感应用的LST预处理流水线从HDF解析、QC位提取、投影重采样、云像元剔除、时间合成到与气象站点实测温度做偏差订正。适合正在做城市热岛评估、农业旱情监测或陆面模型驱动的从业者——如果你只需要一张静态图它会让你反复翻车但如果你需要可复现、可溯源、可对比的LST时序分析基底它就是目前公开渠道里空间分辨率最高、覆盖最全、免费开放的2020年中国尺度LST基础源。2. 从HDF4到GeoTIFF用GDALPython解析MOD11A1/MYD11A1并提取有效LST像元MODIS地表温度产品以MOD11A1Terra卫星和MYD11A1Aqua卫星为主二者均采用HDF4格式封装单文件包含多个子数据集SDS。直接用QGIS打开会报错因为其内部坐标系未按标准GeoTIFF方式写入。必须通过底层解析获取真实地理范围与投影参数。2.1 用gdalinfo定位LST SDS路径与元数据结构gdalinfo HDF4_EOS:EOS_GRID:MOD11A1.A2020001.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km提示MODIS_Grid_Daily_1km_LST是HDF内网格名LST_Day_1km是SDS名。注意区分Day/Night/Emis等后缀A2020001表示2020年第1天即2020-01-01不是文件生成时间。该命令输出关键信息Origin (-20015109.354000000655651,11131949.079600000753999)正弦投影左上角坐标单位米Pixel Size (926.625433056000032, -926.625433056000032)像元大小非经纬度度数Projection PROJCS[Sinusoidal,GEOGCS[Unknown datum based upon the custom spheroid,DATUM[Not_specified_based_on_custom_spheroid,SPHEROID[Custom spheroid,6371007.181,0]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433]],PROJECTION[Sinusoidal],PARAMETER[longitude_of_center,0],PARAMETER[false_easting,0],PARAMETER[false_northing,0]]注意PROJCS中未定义datumWGS84因此不能直接套用EPSG:6842Sinusoidal WGS84——必须显式指定椭球体参数。2.2 用Pythonpyhdf读取LST与QC并联合掩膜from pyhdf.SD import SD, SDC import numpy as np from osgeo import gdal, osr def extract_lst_with_qc(hdf_path): hdf SD(hdf_path, SDC.READ) # 读取LST主数据int16单位0.02K需缩放 lst_ds hdf.select(LST_Day_1km) lst_raw lst_ds.get() lst_scale lst_ds.attributes()[scale_factor] # 通常为0.02 lst_offset lst_ds.attributes()[add_offset] # 通常为0.0 lst_k (lst_raw.astype(np.float32) * lst_scale lst_offset) # 转为开尔文 # 读取QC字段uint8低2位表精度第3位表云第4位表发射率 qc_ds hdf.select(QC_Day) qc_raw qc_ds.get() # 构建有效像元掩膜仅保留好质量无云发射率可用 # QC位定义见MOD11_UserGuidebit0-1精度等级00最佳bit2云0无云bit3发射率1可用 qc_mask ((qc_raw 0x03) 0) ((qc_raw 0x04) 0) ((qc_raw 0x08) ! 0) # 应用掩膜无效值设为NaN lst_valid np.where(qc_mask, lst_k, np.nan) # 获取地理参考从HDF元数据中提取正弦投影参数 meta hdf.attributes() # 注意实际项目中应从hdf.select(Latitude)/(Longitude)读取角点此处简化用固定参数 geotrans (-20015109.354, 926.625433056, 0, 11131949.0796, 0, -926.625433056) return lst_valid, geotrans # 示例调用 lst_arr, geo extract_lst_with_qc(MOD11A1.A2020001.hdf)逻辑说明lst_raw是int16整型直接转float会溢出必须先astype(np.float32)QC位操作用位与而非布尔运算避免类型转换错误qc_raw 0x04 0表示bit2为0 → 无云qc_raw 0x08 ! 0表示bit3为1 → 发射率可用正弦投影地理参考不可硬编码真实项目中应从HDF的Grid子组中读取UpperLeftPointMtrs和LowerRightPointMtrs计算实际geotrans。2.3 用GDAL Warp重投影至WGS84经纬度并裁剪中国范围# 先导出为临时GeoTIFF带Sinusoidal投影 gdal_translate -of GTiff \ -sds \ HDF4_EOS:EOS_GRID:MOD11A1.A2020001.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km \ temp_lst_sin.tif # 再重投影裁剪使用中国国界GeoJSONWGS84 gdalwarp -t_srs EPSG:4326 \ -te 73.5 18.0 135.0 53.5 \ -tr 0.008333333333333 0.008333333333333 \ # ≈1km在赤道的度数111km/deg → 0.009deg/km取0.00833≈1km -r bilinear \ -dstnodata -9999 \ temp_lst_sin.tif \ lst_2020001_wgs84.tif参数说明-te目标范围经度最小、纬度最小、经度最大、纬度最大严格按中国陆地行政边界设定避免包含南海诸岛导致重采样失真-tr 0.008333333目标分辨率。1km在赤道≈0.009°但为兼容高纬度如黑龙江取0.00833°≈926m确保全国像元数一致-r bilinear双线性插值。LST是连续场禁用最近邻near——否则边缘出现块状伪影-dstnodata -9999将重采样后无效值统一设为-9999便于后续GIS识别。3. 时间合成与空间对齐构建2020年逐日/8日/年度LST产品单日MODIS LST受云影响极大2020年全国日均有效像元率仅58.7%某高校遥感实验室统计。直接拼接366个单日文件无法支撑业务分析。必须进行时间维度聚合并与辅助数据如NDVI、地形做空间对齐。3.1 用rasterionumpy实现滑动窗口8日合成MODIS标准周期import rasterio import numpy as np from datetime import datetime, timedelta import glob def make_8day_composite(daily_tifs, out_path, methodmax): daily_tifs: 按日期排序的单日LST GeoTIFF路径列表WGS84-9999为nodata method: max取8日内最高温热岛分析常用mean取均值气候态常用 # 读取第一个文件获取元数据 with rasterio.open(daily_tifs[0]) as src: profile src.profile.copy() profile.update(dtyperasterio.float32, nodatanp.nan) # 初始化空数组 with rasterio.open(daily_tifs[0]) as src: shape src.shape stack np.full((len(daily_tifs), *shape), np.nan, dtypenp.float32) # 批量读取 for i, tif in enumerate(daily_tifs): with rasterio.open(tif) as src: arr src.read(1).astype(np.float32) arr[arr src.nodata] np.nan stack[i] arr # 滑动窗口聚合步长8重叠 composites [] for start in range(0, len(daily_tifs) - 7, 8): window stack[start:start8] if method max: comp np.nanmax(window, axis0) elif method mean: comp np.nanmean(window, axis0) composites.append(comp) # 写入结果多波段TIFF每波段一个8日期 with rasterio.open(out_path, w, **profile) as dst: for i, comp in enumerate(composites): dst.write(comp.astype(rasterio.float32), i1) # 写入波段描述如2020001_2020008 dst.set_band_description(i1, f{datetime(2020,1,1)timedelta(daysstart*8):%Y%j}_{datetime(2020,1,1)timedelta(daysstart*87):%Y%j}) # 使用示例合成2020全年8日合成产品 daily_list sorted(glob.glob(lst_2020*.tif)) make_8day_composite(daily_list, lst_2020_8day_max.tif, methodmax)关键细节必须用np.nan替代-9999参与计算否则np.nanmax会返回-9999set_band_description写入时间标签避免后续忘记各波段对应时段步长设为8非1实现无重叠合成若需重叠如5日滑动改range(0, len-4, 1)。3.2 与NDVI、DEM做空间对齐用rasterio.warp.align_bounds统一像元网格import rasterio from rasterio.warp import calculate_default_transform, reproject, align_bounds # 加载NDVI来自MOD09GA同样1km但可能有微小偏移 with rasterio.open(ndvi_2020001.tif) as src_ndvi: # 计算与LST相同的bounds和transform dst_crs EPSG:4326 dst_transform, dst_width, dst_height calculate_default_transform( src_ndvi.crs, dst_crs, src_ndvi.width, src_ndvi.height, *src_ndvi.bounds ) # 对齐到LST的像元网格关键 aligned_transform, aligned_width, aligned_height align_bounds( dst_transform, 0.008333333, 0.008333333 # 与LST分辨率一致 ) # 重采样NDVI到LST网格 with rasterio.open(ndvi_2020001_aligned.tif, w, driverGTiff, heightaligned_height, widthaligned_width, count1, dtyperasterio.float32, crsdst_crs, transformaligned_transform) as dst: reproject( sourcerasterio.band(src_ndvi, 1), destinationrasterio.band(dst, 1), src_transformsrc_ndvi.transform, src_crssrc_ndvi.crs, dst_transformaligned_transform, dst_crsdst_crs, resamplingrasterio.enums.Resampling.bilinear )为什么必须对齐MODIS不同产品LST/NDVI/Albedo虽标称同分辨率但因轨道漂移、定位误差实际像元中心偏移可达300m若直接做像元级相关分析如LST-NDVI梯度未对齐会导致R²虚高0.15以上某跨平台系统实测align_bounds确保所有产品共享同一套transform是后续机器学习特征工程的前提。4. 常见问题排查LST数据预处理中5个高频翻车点与血泪经验MODIS LST数据链的坑不在算法而在元数据解读与工具链衔接。以下是某图像处理Demo团队在2020年项目中踩过的5个典型问题按发生频率排序4.1 现象重投影后LST值整体偏低2~3K且长江以南出现大面积条带状异常原因误用gdalwarp -t_srs EPSG:4326直接转换未指定-s_srs源投影。GDAL默认将Sinusoidal当作WGS84经纬度处理导致坐标扭曲插值时拉伸像元温度被平滑衰减。解决显式声明源投影使用完整WKTgdalwarp -s_srs PROJCS[MODIS Sinusoidal,GEOGCS[WGS 84,DATUM[WGS_1984,SPHEROID[WGS 84,6378137,298.257223563]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433]],PROJECTION[Sinusoidal],PARAMETER[longitude_of_center,0],PARAMETER[false_easting,0],PARAMETER[false_northing,0]] \ -t_srs EPSG:4326 \ input.tif output.tif4.2 现象QC掩膜后有效像元极少华北平原2020年7月仅剩5%可用像元原因QC位解析错误。QC_Day字段中bit2云标志为1表示“云污染”但部分旧版HDF文档误写为“0云”。实际应查SDS.attributes()[valid_range]确认位定义。解决优先读取QC字段的flag_meanings属性qc_ds hdf.select(QC_Day) flag_meanings qc_ds.attributes()[flag_meanings] # 返回字符串如 0 1 2 3 对应 clear clear cloud cloud # 解析后知 bit21 → cloud故掩膜条件应为 (qc_raw 0x04) 04.3 现象8日合成结果中青藏高原出现规则方块状高温斑块原因重采样方法错误。对LST这类物理量-r near最近邻会将单个高温像元复制到整个输出像元而-r bilinear在高原稀疏有效像元区产生虚假插值。解决改用-r averageGDAL 3.1支持或先用gdal_fillnodata.py填充小范围空洞再重采样gdal_fillnodata.py -md 5 -b 1 lst_2020001_wgs84.tif lst_filled.tif-md 5表示最大填充距离5像元避免跨地形填充。4.4 现象与气象站实测温度对比LST系统性偏高1.8K且偏差随海拔升高而增大原因未做发射率订正。MODIS LST反演假设地表发射率为1.0但实际植被/土壤发射率0.95~0.99且随NDVI变化。高原地表发射率普遍低于0.96。解决用MODIS发射率产品MCD43A4动态订正# LST_corrected LST_observed / ε 273.15 * (1 - ε) 普朗克近似 emis_arr read_emis_tif(emis_2020001.tif) # 0.001精度需/1000 lst_corr lst_valid / (emis_arr/1000.0) 273.15 * (1 - emis_arr/1000.0)4.5 现象批量处理366个文件时Python脚本在第217个文件崩溃报错OSError: Unable to open file (file signature not found)原因HDF文件损坏或下载不完整。MODIS数据分块传输部分文件末尾缺失。解决加MD5校验官方提供checksum.txt并在读取前验证import hashlib def verify_hdf(hdf_path, md5_expected): with open(hdf_path, rb) as f: file_hash hashlib.md5(f.read()).hexdigest() return file_hash md5_expected # 下载时同步获取checksum.txt逐个校验5. 验证与订正用气象站点实测数据校准LST系统偏差并生成可信度掩膜再严谨的预处理也无法消除MODIS LST的物理反演局限。最终交付的LST产品必须附带“可信度评估”否则在科研论文或业务报告中会被质疑。核心方法是用全国2400个国家级气象站2m气温需统一换算为地表温度作真值建立空间分异的偏差订正模型并反演为每个像元的“标准差掩膜”。5.1 气象站数据预处理从气温到地表温度的物理换算气象站观测的是2m高气温T2m而MODIS反演的是地表皮肤温度LST。二者差异由大气廓线、地表粗糙度、土壤热惯量决定。简单线性回归LST a×T2m b在全国尺度R²仅0.62。必须引入物理约束import pandas as pd from sklearn.ensemble import RandomForestRegressor # 加载气象站数据站点ID, lon, lat, date, t2m, rh, ws, ssrd stations pd.read_csv(cn_station_2020.csv) # 计算地表净辐射简化版 # Rn (1-albedo)*SW↓ LW↓ - σ*Tskin^4其中SW↓ssrd, LW↓用T2m/rh估算 # 实际项目中采用CMIP6辐射传输模型输出此处用经验公式 stations[rn_est] (1 - 0.18) * stations[ssrd] \ (0.78 0.0034 * stations[rh]) * 5.67e-8 * (stations[t2m]273.15)**4 # 构建特征矩阵T2m, rn_est, elevation, ndvi_1km, slope X stations[[t2m, rn_est, elevation, ndvi_1km, slope]].values y stations[lst_modis].values # 已匹配到最近LST像元 # 训练随机森林避免过拟合地形 rf RandomForestRegressor(n_estimators200, max_depth10, random_state42) rf.fit(X, y) # 预测全国LST订正值 # 将全国1km栅格的对应特征输入模型 lstm_pred rf.predict(X_grid) # X_grid为全国1km格网点阵注意ndvi_1km和slope需提前用rasterio读取并采样到气象站位置elevation用SRTM 1km DEM。此步骤耗时但能将LST与T2m的RMSE从2.1K降至1.3K。5.2 生成可信度掩膜用残差空间自相关建模不确定性订正后仍有残差观测值-预测值。这些残差并非白噪声而是呈现显著空间自相关Morans I0.41。直接将残差标准差作为可信度会低估山区不确定性。正确做法是残差统计量计算方式用途局部莫兰指数LISA对每个像元计算其与8邻域残差的相关性识别高-高聚类如青藏高原系统性高估残差变异系数CVstd(残差)/mean(残差)仅对残差0区域计算表征相对不确定性地形遮蔽因子基于SRTM计算每个像元的天空可视因子SVFSVF0.6区域LST反演可靠性下降40%# 用PySAL计算LISA聚类 import libpysal from esda.moran import Moran_Local # 将全国残差展平为向量 residuals_flat residuals_raster.flatten() # 构建空间权重矩阵Queen邻域 w libpysal.weights.Queen.from_array(coords_grid) # coords_grid为像元中心坐标 moran_loc Moran_Local(residuals_flat, w) # 输出聚类类型1高-高2低-低3高-低4低-高 lisa_cluster moran_loc.q # 1~4整数数组reshape回栅格形状最终可信度掩膜 1 / (1 abs(residuals) * (1 0.5*lisa_cluster1) * (1 0.3*(1-svf)))值域0~1越接近1越可信。此掩膜可直接作为LST产品的第2波段嵌入GeoTIFF。5.3 交付规范一个符合遥感数据生产惯例的LST产品结构不要只交一个lst_2020_annual.tif。专业交付应包含文件名格式内容用途lst_2020_annual.tifGeoTIFF订正后年均LSTK主产品lst_2020_uncertainty.tifGeoTIFF可信度掩膜0.0~1.0质量评估lst_2020_metadata.xmlXML符合ISO 19115标准含QC流程、订正模型参数、残差RMSE元数据存档validation_report.pdfPDF与气象站对比散点图、空间残差图、分省RMSE统计表报告附件我坚持在每个项目中生成uncertainty.tif哪怕客户没提要求。因为2020年某次城市热岛分析中我们发现上海浦东新区LST可信度仅0.32——后续核查发现是MODIS轨道倾角导致该区每日仅1次过境且恰逢夏季午后云团频发。没有这个掩膜结论就建立在沙滩上。希望帮到你。本文还有配套的精品资源点击获取