ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

MODIS地表温度LST数据处理:QC解码、Albers投影与UQ加权合成

MODIS地表温度LST数据处理:QC解码、Albers投影与UQ加权合成 简介本资源为2020年中国全域1km分辨率地表温度LST空间分布数据集面向遥感、地理信息、生态与气候研究领域的科研人员及GIS初学者用于支撑区域热环境分析、城市热岛评估、地表能量平衡建模等空间分析任务。数据基于NASA MODIS MOD11A2产品8天合成、1km分辨率经子区提取、影像拼接、Albers等积圆锥投影中央经线105°标准纬线25°/47°、单位换算含开氏与摄氏双版本TIFF及年度均值合成处理而成具备明确坐标系与规范元数据。压缩包共11个文件含2个核心GeoTIFF栅格kelvin/celsius双温标、2个TFW地理配准文件、4个XML元数据与辅助文件、2个说明文本含数据来源、处理逻辑与温度单位说明及1个OVF金字塔文件总大小66.16MB结构完整、即下即用。目前已有3545人学习下载用户可直接加载至ArcGIS/QGIS开展空间统计、时序对比或与NDVI、土地利用数据叠加分析无需额外预处理即可投入科研与教学实践。1. MODIS 2020年中国1km地表温度LST空间分布数据集不是“拿来就能画图”的栅格包而是需要校验坐标系、重采样策略和像元有效性标记的实测级遥感产品你手头拿到的这个“MODIS 2020年中国1km LST”数据集表面看是一组GeoTIFF或HDF5文件但实际它是一套经过严格质量控制、带多层掩膜与不确定性评估的科学级地表温度产品——不是气象站插值结果也不是模型模拟输出而是Terra/Aqua双星搭载的MODIS传感器在2020年全年对中国陆域含近海岛屿逐日反演的真实热红外观测。它能支撑城市热岛强度量化、农业旱情早期识别、冻土退化趋势分析等任务但前提是必须理解其空间参考框架Albers等面积投影 vs WGS84经纬度、时间合成逻辑日值/8天/月均、以及最关键的QC位标志解码规则。适合已具备GDAL/Rasterio基础、能读取HDF5属性、会处理NoData掩膜的遥感数据处理者不适合仅会用ArcGIS“拖进去就出图”的新手。很多用户下载后第一反应是“为什么打开全是-9999”——这不是损坏而是原始QC字段未解译也有人直接叠加到WGS84底图上发现偏移3–5公里这是Albers投影参数未对齐的典型表现。这份数据的价值不在“有”而在“怎么用对”。2. 数据结构解析与核心字段解码从HDF5元数据到LST主变量、QC掩膜、不确定性估值的三层映射MODIS LST产品以MYD11A1/MOD11A1 V6.1为例采用HDF-EOS5格式封装单个HDF文件内含多个SDSScientific Data Set和全局属性。中国区域1km数据集并非简单裁剪而是基于MODIS全球L2G产品经重投影、重采样、云剔除与地理配准后生成的衍生集。理解其内部结构是后续所有操作的前提。2.1 HDF5文件层级与关键SDS定位使用h5py可快速探查文件结构。以下为典型中国区域HDF5文件的顶层组织import h5py f h5py.File(MYD11A1.A2020001.h10v04.061.2020002032732.hdf, r) print(list(f.keys())) # 输出: [HDFEOS, Metadata] print(list(f[HDFEOS][GRIDS][MODIS_Grid_LST].keys())) # SDS列表关键SDS包括LST_Day_1km/LST_Night_1km地表温度主变量单位为0.02K需乘以0.02转换为开尔文K再减273.15得摄氏度℃。注意该值为亮温Brightness Temperature反演结果非辐射定标后的原始DN值。QC_Day/QC_Night8位整型质量控制字段每个bit代表不同物理条件如云、水汽、发射率误差、角度等需按NASA官方文档《MOD11_UserGuide_V6.1》Table 7解码。例如bit0–1表示“算法完成状态”bit2–3表示“云检测结果”bit4–5表示“发射率误差等级”。Emis_1km白天/夜间发射率估算值波段31/32单位为0.001需除以1000。Clear_day_cov/Clear_night_cov晴空像元覆盖比例0–100用于判断时间序列中有效观测密度。提示不要跳过f.attrs[CoreMetadata.0]和f[HDFEOS][ADDITIONAL][FILE_ATTRIBUTES]中的全局属性。其中Projection字段明确声明为Albers Equal Area Conic中央经线105°E标准纬线25°N/47°N椭球体WGS84——这是后续所有坐标系转换的基准不可用默认WGS84地理坐标系强行打开。2.2 QC位字段逐bit解码实战识别“可用但需谨慎”的像元QC字段是LST数据可信度的核心判据。直接丢弃所有非0值会损失大量有效数据全盘接受则引入系统性偏差。正确做法是构建分层掩膜import numpy as np def decode_qc(qc_array: np.ndarray) - dict: 输入uint8 QC数组如QC_Day 输出字典含各bit含义布尔掩膜 参考NASA MOD11 User Guide V6.1 Table 7 (bit0LSB) qc_uint8 qc_array.astype(np.uint8) masks {} # bit0-1: Algorithm completion flag (00good, 01other, 10not processed, 11not processed) masks[algo_ok] ((qc_uint8 0b11) 0) # bit2-3: Cloud mask (00clear, 01cloudy, 10uncertain, 11cloudy) cloud_bits (qc_uint8 2) 0b11 masks[cloud_clear] (cloud_bits 0) masks[cloud_uncertain] (cloud_bits 2) # bit4-5: Emissivity error (00low, 01medium, 10high, 11not retrieved) emis_err (qc_uint8 4) 0b11 masks[emis_err_low] (emis_err 0) # bit6-7: LST error (00low, 01medium, 10high, 11not retrieved) lst_err (qc_uint8 6) 0b11 masks[lst_err_low] (lst_err 0) return masks # 应用示例 qc_day f[HDFEOS][GRIDS][MODIS_Grid_LST][Data Fields][QC_Day][:] masks decode_qc(qc_day) valid_mask masks[algo_ok] masks[cloud_clear] masks[emis_err_low] masks[lst_err_low] lst_day_k f[HDFEOS][GRIDS][MODIS_Grid_LST][Data Fields][LST_Day_1km][:] * 0.02 lst_day_c lst_day_k - 273.15 lst_day_c[~valid_mask] np.nan # 仅保留高置信度像元这段代码的关键在于它没有简单用QC_Day 0做一刀切而是分维度评估。例如cloud_uncertain区域在干旱区可能仍具参考价值而emis_err_high在植被覆盖区影响较小但在裸土区会导致±3℃以上偏差。这种解码逻辑直接决定了你后续统计结果的物理意义是否成立。2.3 不确定性估值Uncertainty字段的物理意义与使用边界V6.1版本新增UQ_Day_1km和UQ_Night_1km字段单位为0.02K表示LST反演的1σ标准差。这不是噪声水平而是综合了大气校正残差、发射率假设误差、角度效应及仪器定标不确定性的传播结果。其值域通常为0–5K但在高纬度冬季或强逆温条件下可达8K以上。UQ值区间K典型场景建议用途1.0晴朗夏季平原、中低纬度农田可用于亚像元尺度变化检测如灌溉前后ΔT1.0–2.5丘陵林区、春季融雪期适用于区域均值统计不建议单点分析2.5高山积雪区、冬季城市边缘、薄云边缘仅作趋势指示需结合地表类型掩膜降权使用注意UQ字段本身也受QC约束——若QC_Day中algo_ok为False则对应UQ值无效。实践中我一般先用QC生成基础掩膜再对UQ进行3×3窗口中值滤波以抑制孤立噪点最后按UQ值对LST做加权平均权重1/UQ²比简单均值更能反映真实热环境空间格局。3. 空间坐标系校准与重投影Albers投影参数验证、WGS84地理坐标系转换陷阱与国产GIS软件兼容方案中国1km LST数据集采用Albers等面积圆锥投影Albers Equal Area Conic这是为保障面积统计精度如热岛面积占比、冻融范围变化而设定的强制要求。但多数开源工具默认将其识别为WGS84地理坐标系导致空间错位、面积计算失真、与矢量行政边界套合失败——这是新手踩坑最密集的环节。3.1 Albers投影参数精确提取与GDAL验证HDF文件中Projection属性仅给出名称具体参数需从GridOrigin和ProjParams中解析。使用GDAL命令行工具可直接读取gdalinfo MYD11A1.A2020001.h10v04.061.2020002032732.hdf关键输出片段Coordinate System is: PROJCRS[unnamed, BASEGEOGCRS[WGS 84, DATUM[World Geodetic System 1984, ELLIPSOID[WGS 84,6378137,298.257223563, LENGTHUNIT[metre,1]]], PRIMEM[Greenwich,0, ANGLEUNIT[degree,0.0174532925199433]]], CONVERSION[Albers Equal Area Conic, METHOD[Albers Equal Area Conic, ID[EPSG,9822]], PARAMETER[Latitude of false origin,0, ANGLEUNIT[degree,0.0174532925199433], ID[EPSG,8821]], PARAMETER[Longitude of false origin,105, ANGLEUNIT[degree,0.0174532925199433], ID[EPSG,8822]], PARAMETER[Latitude of 1st standard parallel,25, ANGLEUNIT[degree,0.0174532925199433], ID[EPSG,8823]], PARAMETER[Latitude of 2nd standard parallel,47, ANGLEUNIT[degree,0.0174532925199433], ID[EPSG,8824]], PARAMETER[Easting at false origin,0, LENGTHUNIT[metre,1], ID[EPSG,8826]], PARAMETER[Northing at false origin,0, LENGTHUNIT[metre,1], ID[EPSG,8827]]], CS[Cartesian,2], AXIS[easting,east, ORDER[1], LENGTHUNIT[metre,1]], AXIS[northing,north, ORDER[2], LENGTHUNIT[metre,1]]]重点确认四参数false_origin_lon105,std_parallel_125,std_parallel_247,datumWGS84。任何一项错误都会导致投影偏移。常见错误是将std_parallel_1误读为central_meridian或忽略false_origin为(0,0)意味着原点在投影平面中心而非经纬度(105,0)。3.2 安全重投影至WGS84地理坐标系的三步法直接用gdalwarp -t_srs EPSG:4326会因重采样算法选择不当引入几何畸变。推荐流程先提取Albers坐标系下的真实地理范围以米为单位gdalinfo -proj4 MYD11A1.A2020001.h10v04.061.2020002032732.hdf # 输出类似projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsWGS84 unitsm no_defs用gdal_translate导出带完整坐标系定义的GeoTIFF并指定-a_srs显式写入gdal_translate \ -of GTiff \ -a_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsWGS84 unitsm no_defs \ -co COMPRESSLZW \ HDF4_EOS:EOS_GRID:MYD11A1.A2020001.h10v04.061.2020002032732.hdf:MODIS_Grid_LST:LST_Day_1km \ lst_day_albers.tif再用gdalwarp重投影强制使用-r bilinear双线性并设置-tr 0.01 0.01约1km分辨率gdalwarp \ -t_srs EPSG:4326 \ -r bilinear \ -tr 0.01 0.01 \ -te 73.5 18.0 135.5 53.5 \ # 中国四至经纬度WGS84 -co COMPRESSLZW \ lst_day_albers.tif \ lst_day_wgs84.tif关键细节-tr 0.01 0.01是经验性设置1°≈111km0.01°≈1.1km比-tr 0.008983 0.008983严格1km更鲁棒-te必须用WGS84经纬度范围不能用Albers米制范围-r bilinear比-r near更能保持温度梯度连续性避免阶梯状伪影。3.3 国产GIS平台如SuperMap、MapGIS坐标系识别失败的绕过方案部分国产GIS软件无法自动解析HDF5内嵌的Albers参数显示为“未知坐标系”。此时可手动创建.prj文件PROJCS[Albers_Equal_Area_Conic, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137.0,298.257223563]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.0174532925199433]], PROJECTION[Albers], PARAMETER[False_Easting,0.0], PARAMETER[False_Northing,0.0], PARAMETER[Central_Meridian,105.0], PARAMETER[Standard_Parallel_1,25.0], PARAMETER[Standard_Parallel_2,47.0], PARAMETER[Latitude_Of_Origin,0.0], UNIT[Meter,1.0]]将此内容保存为albers_china.prj与GeoTIFF同名放置。SuperMap 10i及以上版本可正确加载。MapGIS需在“坐标系管理器”中导入该.prj并绑定至图层。4. 时间序列构建与质量加权合成从日值到月均的缺失值填补、云污染剔除与空间一致性校验单日LST受云、大气水汽、太阳高度角影响极大直接使用日值易产生虚假突变。科学做法是构建8天或月合成产品但合成过程绝非简单np.nanmean()——必须融合QC掩膜、UQ权重、时间有效性如某日无有效观测则不参与合成。4.1 日值→8天合成NASA官方合成逻辑复现MODIS官方8天产品MYD11A2采用“滑动窗口最大值合成法Maximum Value Composite, MVC”但针对LST做了优化优先选择QC最高、UQ最小、观测时间最接近地方正午的像元。我们可在本地复现该逻辑import xarray as xr import numpy as np def composite_8day(daily_files: list, var_name: str LST_Day_1km) - xr.DataArray: 输入按时间排序的8个HDF文件路径列表如2020001–2020008 输出xarray DataArray含坐标、属性、合成后LST℃及UQ lst_list, qc_list, uq_list [], [], [] for fp in daily_files: with h5py.File(fp, r) as f: # 读取主变量K与QC、UQ lst_k f[HDFEOS][GRIDS][MODIS_Grid_LST][Data Fields][var_name][:] qc f[HDFEOS][GRIDS][MODIS_Grid_LST][Data Fields][fQC_{var_name.split(_)[1]}][:] uq_k f[HDFEOS][GRIDS][MODIS_Grid_LST][Data Fields][fUQ_{var_name.split(_)[1]}][:] # 转换为℃并应用QC掩膜仅保留algo_ok且cloud_clear masks decode_qc(qc) valid masks[algo_ok] masks[cloud_clear] lst_c (lst_k * 0.02 - 273.15) lst_c[~valid] np.nan uq_c uq_k * 0.02 # UQ单位也是0.02K lst_list.append(lst_c) uq_list.append(uq_c) # 转为xarray便于时间维度操作 da_lst xr.DataArray( np.stack(lst_list, axis0), dims[time, y, x], coords{time: [int(fp[-13:-7]) for fp in daily_files]} # 儒略日 ) da_uq xr.DataArray( np.stack(uq_list, axis0), dims[time, y, x] ) # 加权合成权重 1 / (UQ² ε)ε0.01避免除零 weights 1 / (da_uq**2 0.01) weights weights / weights.sum(dimtime) # 归一化 # 合成LST Σ(weight_i × lst_i) lst_composite (da_lst * weights).sum(dimtime) return lst_composite # 使用示例合成2020年第1个8天周期 files_8day [ MYD11A1.A2020001.h10v04.061.2020002032732.hdf, MYD11A1.A2020002.h10v04.061.2020003032732.hdf, # ... 至 2020008 ] lst_8day composite_8day(files_8day)该函数的核心优势在于它没有丢弃UQ高的像元而是降低其权重同时通过decode_qc确保云污染像元完全排除。相比np.nanmean()合成结果在云区边缘更平滑在晴空区更稳定。4.2 月均合成中的“有效观测日数”阈值设定与空间校验月均LST要求当月至少有10天有效观测NASA标准否则该月标记为“数据不足”。但中国西部高原、东北冬季常低于此阈值。此时需空间校验def monthly_composite_with_validation( daily_files: list, min_valid_days: int 10, spatial_consistency_window: int 5 ) - tuple[xr.DataArray, np.ndarray]: 返回月均LST DataArray 有效日数掩膜0不足1充足 spatial_consistency_window: 对每个像元检查其5×5邻域内有效日数方差 若方差 3则判定为“空间异常”即使日数达标也标记为无效 # 步骤1统计每个像元的有效日数QC合格且非nan valid_count np.zeros_like(lst_list[0], dtypeint) for lst_c in lst_list: valid_count ~np.isnan(lst_c) # 步骤2空间一致性检验5×5窗口内方差 from scipy.ndimage import generic_filter variance_map generic_filter( valid_count.astype(float), lambda x: np.var(x), sizespatial_consistency_window ) # 步骤3生成最终掩膜 day_mask (valid_count min_valid_days) space_mask (variance_map 3) final_mask day_mask space_mask # 步骤4仅对final_maskTrue区域计算月均 lst_stack np.stack(lst_list, axis0) monthly_mean np.nanmean(lst_stack, axis0) monthly_mean[~final_mask] np.nan return xr.DataArray(monthly_mean, dims[y,x]), final_mask # 输出final_mask可用于后续空间插值或标记“数据不可靠区”此方法避免了“某县全域因1个像元有效而全境被赋值”的粗放做法符合遥感产品真实性检验规范。5. 常见问题排查5类高频翻车现场与血泪经验总结MODIS LST数据看似标准化但实际处理中存在大量隐性陷阱。以下是我在多个项目中反复验证的5类典型问题每条均按“现象→原因→解决”结构给出可立即执行的方案。5.1 现象打开GeoTIFF后LST值全为-9999或极小负数如-273.15原因未对原始LST变量单位0.02K进行缩放转换或误将QC掩膜应用于错误变量。HDF中LST_Day_1km存储的是整型DN值直接读取即为大负数。解决确认读取的是LST_Day_1km而非LST_Day_1km_QC强制执行lst_c (lst_dn * 0.02) - 273.15在应用QC前先用np.unique(lst_dn)检查DN值域是否在[0, 65535]内V6.1有效范围若出现-32768等值说明HDF读取异常需重下文件。5.2 现象Albers投影数据与中国省级矢量边界如1:100万shp套合时东部沿海偏移10–20公里原因矢量边界shp的.prj文件中定义的Albers参数与MODIS不一致——常见错误是lat_1/lat_2设为25°/45°旧版或lon_0设为110°北京中心。解决用ogrinfo -so province.shp检查矢量坐标系若参数不匹配用ogr2ogr -t_srs projaea lat_125 lat_247 lon_0105 datumWGS84重投影矢量或更稳妥将LST转为WGS84后再与WGS84矢量叠加见3.2节。5.3 现象同一地点不同日期LST值突变超过15℃如晴天35℃→阴天20℃→次日又35℃时间序列呈锯齿状原因未使用QC字段剔除cloud_uncertain或emis_err_medium像元导致部分“疑似云”或“中等发射率误差”像元混入而这些像元在反演中易受大气水汽扰动。解决严格限定QC条件algo_ok cloud_clear (emis_err 0) (lst_err 0)对时间序列添加3天移动平均滤波scipy.signal.savgol_filter窗口长度5阶数2绘制UQ时间序列若UQ同步突变证实为反演不确定性主导。5.4 现象重投影后WGS84 GeoTIFF在QGIS中显示为“马赛克块”尤其在省界附近出现明显接缝原因gdalwarp未指定-taptarget aligned pixels参数导致不同HDF瓦片重投影后像素网格未对齐。解决在gdalwarp命令中加入-tap -tr 0.01 0.01强制所有输出文件像素左上角对齐同一经纬度网格或先用gdalbuildvrt构建虚拟镶嵌再统一重投影gdalbuildvrt mosaic.vrt *.hdf gdalwarp -t_srs EPSG:4326 -tap -tr 0.01 0.01 mosaic.vrt mosaic_wgs84.tif5.5 现象用Pythonrasterio读取HDF5时抛出KeyError: HDF4_EOS:EOS_GRID或Driver not registered原因rasterio默认安装不包含HDF4驱动或GDAL版本3.3HDF4支持不稳定。解决优先使用h5py直接读取见2.1节绕过rasterio若必须用rasterio安装conda install -c conda-forge gdal rasterioconda-forge版本含HDF4或改用xarray.open_rasterio()其对HDF5兼容性更好import xarray as xr ds xr.open_rasterio(file.hdf, enginerasterio) # 注意需指定子数据集名如ds.sel(band1) or ds.isel(band0)6. 进阶技巧构建LST异常指数LSTAI实现旱情动态监测与城市热岛强度量化LST本身是绝对温度但单一数值难以直接反映异常程度。真正有业务价值的是相对于气候态的偏离——即LST异常指数Land Surface Temperature Anomaly Index, LSTAI。它能将2020年LST数据转化为可跨年份、跨区域比较的标准化指标支撑农业旱情预警与城市规划决策。6.1 LSTAI计算公式与气候态基准构建LSTAI定义为LSTAI LST2020− LSTclim其中LST_clim是2001–2020年同期如7月的多年平均LST。关键在于LST_clim必须与2020年数据在相同QC标准、相同重采样方法、相同空间分辨率下生成否则偏差可达2–5℃。我采用以下稳健流程构建气候态统一QC标准对2001–2020年全部MYD11A1数据使用2.2节decode_qc()函数仅保留algo_ok cloud_clear emis_err_low lst_err_low像元逐月合成对每年7月1–31日用4.1节composite_8day()函数生成3个8天合成再取平均得月均空间聚合对每个1km像元计算20年7月均值的中位数而非均值规避极端年份如2006年川渝大旱影响生成气候态栅格输出为GeoTIFF坐标系与2020年数据完全一致Albers。# 示例计算2020年7月LSTAI lst_2020_july xr.open_rasterio(lst_202007_albers.tif).squeeze() lst_clim_july xr.open_rasterio(lst_clim_07_albers.tif).squeeze() # 确保二者空间范围、分辨率完全一致必要时重采样 if not (lst_2020_july.shape lst_clim_july.shape): lst_clim_july lst_clim_july.interp_like(lst_2020_july, methodnearest) lstai_july lst_2020_july - lst_clim_july lstai_july.rio.to_raster(lstai_202007.tif, compressLZW)6.2 LSTAI在农业旱情监测中的阈值分级与验证LSTAI正值表示地表温度高于常年常伴随土壤水分亏缺。根据某高校农业遥感团队在华北平原的实测验证分级阈值如下LSTAI℃干旱等级土壤湿度状态0–1典型作物响应 0.5无旱0.35正常生长0.5–1.5轻旱0.25–0.35拔节期玉米叶片轻微卷曲1.5–2.5中旱0.15–0.25小麦灌浆速率下降20%2.5重旱0.15棉花蕾铃脱落率30%验证方法选取100个农田样点同步采集0–20cm土壤湿度ML2x传感器与LSTAI计算Spearman相关系数。在华北平原7–8月LSTAI与土壤湿度呈显著负相关ρ -0.72, p0.001证明其物理合理性。6.3 城市热岛强度UHI量化LSTAI空间梯度与建成区缓冲区分析城市热岛不是“城市比农村热”而是城市核心区与周边乡村的LSTAI差值。但简单取矩形框会受地形干扰。我采用以下三步法构建建成区缓冲区用2020年全球不透水面数据如GAIA提取城市建成区生成5km、10km、15km缓冲区计算各缓冲区LSTAI均值# 假设buffer_5km为二值掩膜1缓冲区内 lstai_5km_mean lstai_july.where(buffer_5km).mean().item() lstai_rural_mean lstai_july.where(~buffer_15km).mean().item() # 15km外为乡村 uhi_intensity lstai_5km_mean - lstai_rural_mean空间梯度分析沿城市中心到乡村方向每1km提取LSTAI剖面拟合线性趋势——斜率0.15℃/km即判定为强热岛。从那以后我每次处理MODIS LST都强制走一遍QC位解码Albers参数验证UQ加权合成三步检查。哪怕项目 deadline压得只剩48小时这15分钟也绝不省——因为一次坐标系错误可能导致整个县域的热岛面积统计偏差37%而返工重跑要耗掉两天。希望帮到你。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进