ARTICLE · INTELLIGENCE

战地情报 · 详情页

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

MODIS 2022中国1km地表温度数据预处理全指南

MODIS 2022中国1km地表温度数据预处理全指南 简介本资源为2022年中国全域1km分辨率地表温度LST空间分布数据集基于NASA MODIS MOD11A2产品加工生成面向遥感、地理信息、生态与气候研究领域的科研人员及高校师生支撑区域热环境分析、城市热岛效应评估、地表能量平衡建模等应用。数据经子区提取、影像拼接、Albers投影转换、单位换算含开氏与摄氏双版本及年度均值合成最终提供覆盖全国的栅格TIFF文件辅以元数据XML、地理配准TFW、统计说明TXT及金字塔OVF文件共7个文件总大小64.55MB。已有743人学习下载用户可直接加载至ArcGIS/QGIS开展空间分析无需额外预处理包含完整坐标系参数Albers投影中央经线105°标准纬线25°/47°WGS84椭球并附带温度单位说明与数据来源引用规范便于学术引用与成果复现。1. MODIS 2022年中国1km地表温度LST空间分布数据集不是“下载即用”的遥感产品而是需要校验、重投影、掩膜和时空对齐的工程级输入源你拿到这个名为“MODIS 2022年中国1km地表温度LST空间分布数据集.zip”的压缩包时第一反应可能是——终于不用自己下HDF、转TIFF、拼接、裁剪、重采样了恭喜省了3天预处理时间但更要警惕它大概率是某次科研项目产出的中间成果未经标准化质检坐标系混杂常见WGS84地理坐标系 vs. Albers等积投影像元值未统一标定有的存为uint16缩放整数有的直接存K氏浮点云掩膜缺失或粗粒度且2022全年365天的数据未必全量覆盖——尤其青藏高原冬季、东北林区雪盖期、长江中下游梅雨季有效像元率常低于30%。这不是一个开箱即用的“地图图层”而是一份需你亲手校准的遥感工程原材料。适合城市热岛分析、农业旱情初筛、生态模型驱动、气候验证等场景的工程师与研究生不适合只想拖进QGIS点几下就出图的用户。真正价值不在“有数据”而在你能把它变成可复现、可溯源、可嵌入业务流程的LST时间序列变量——这要求你掌握MODIS LST产品的物理约束、MOD11A2/MYD11A2产品结构、中国行政区划矢量对齐逻辑以及1km分辨率下常见的空间错位陷阱。2. 解压与结构解析识别数据组织逻辑区分MOD11A2与MYD11A2定位有效文件与元数据拿到.zip文件后不要急于双击打开。先用命令行确认内部结构——因为很多批量处理脚本依赖路径一致性而Windows资源管理器解压可能自动创建多余父目录。unzip -l MODIS_2022_China_LST_1km.zip | head -20你大概率会看到类似结构MODIS_2022_China_LST_1km/ ├── tif/ │ ├── MOD11A2.A2022001.006.2022012075925.tif │ ├── MOD11A2.A2022001.006.2022012075925.xml │ ├── MYD11A2.A2022001.006.2022012081234.tif │ └── ... ├── shp/ │ └── china_province.shp ├── docs/ │ └── README.md └── metadata/ └── product_version.txt2.1 理解MOD11A2与MYD11A2的本质差异不是“两个版本”而是“两颗卫星、两种过境时间”MOD11A2Terra卫星上午10:30左右过境提供白天LSTDaytime LST对应地表受太阳辐射加热后的峰值温度MYD11A2Aqua卫星下午1:30左右过境提供夜间LSTNighttime LST反映地表长波辐射冷却状态。提示二者不能简单平均城市热岛研究中夜间LST更敏感农业蒸散发估算中白天LST更关键若需日较差Diurnal Temperature Range必须严格配对同一天的MOD11A2白天与MYD11A2夜间——但注意Aqua下午过境时部分区域已进入黄昏辐射定标存在边缘误差需检查QC_Day/QC_Night质量标志波段。2.2 检查TIFF文件的GDAL信息确认坐标系、数据类型、NoData值与波段语义用gdalinfo查看任意一个TIFF以MOD11A2为例gdalinfo tif/MOD11A2.A2022001.006.2022012075925.tif重点关注输出中的这几行Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433]] Origin (73.000000000000000,54.000000000000000) Pixel Size (0.008333333333333,-0.008333333333333) Band 1 Block12000x1 TypeUInt16, ColorInterpGray NoData Value0 Metadata: STATISTICS_MINIMUM7500 STATISTICS_MAXIMUM15000这里暴露三个关键事实坐标系是WGS84地理坐标系经纬度非投影坐标系→ 后续做面积统计、缓冲区分析前必须重投影像元大小0.008333° ≈ 1km在赤道附近成立高纬度实际距离收缩→ 这是MODIS LST产品标准分辨率但注意1km是等角分辨率非等距分辨率数据类型为UInt16NoData值为0且值域7500–15000→ 这是MODIS标准缩放格式真实LSTK DN × 0.02 150.0见MOD11A2官方文档Table 8。若直接当℃读取会得到荒谬的75–150℃。2.3 验证元数据XML与README的一致性避免被“2022年全量”误导打开tif/MOD11A2.A2022001.006.2022012075925.xml搜索RangeBeginningDate和RangeEndingDate确认该文件是否真对应2022年1月1日A2022001。再检查docs/README.md中是否声明数据是否经过云掩膜Cloud Mask处理若未处理QC_Day波段必须参与筛选是否剔除“前景像素”Foreground Pixel以外的像元MODIS LST产品中仅QC_Day 0二进制00000000表示“最佳质量”其余含云、边缘、水体等是否已进行“陆地掩膜”Land/Water Mask原始MODIS LST包含湖泊、水库但多数应用需纯陆地LST。血泪经验某次城市热岛分析中因未检查README直接使用了未掩膜水体的LST导致太湖周边出现虚假高温带——水面LST在夏季可达35℃但城市热岛定义中“地表”特指不透水面与植被水体应排除。3. 坐标系统一与重投影为什么Albers等积投影是分析中国LST的默认选择中国国土跨越5个时区、地形起伏剧烈从海平面到海拔5000m若坚持用WGS84地理坐标系计算面积、做缓冲区、叠加行政边界会产生系统性畸变在黑龙江北部1°经度≈55km在海南1°经度≈110km → 同一像元在不同纬度代表的实际面积差一倍圆形缓冲区在高纬度拉成椭圆影响热岛半径测算与省级矢量如shp/china_province.shp叠加时因矢量常为Albers投影地理坐标系TIFF会轻微错位尤其边境线。3.1 选择Albers Equal Area Conic中国常用而非UTM或Web MercatorAlbers投影在中国全境变形最小且保持面积比例不变——这对计算“某省高温像元占比”、“城市建成区平均LST”至关重要。ESRI定义的中国Albers参数为标准纬线125°N标准纬线247°N中央经线105°E原点纬度0°N单位米用GDAL重投影保留原始分辨率不重采样gdalwarp \ -t_srs projaea lat_125 lat_247 lon_0105 datumWGS84 unitsm no_defs \ -r near \ # 最近邻重采样避免LST值被插值污染 -tr 1000 1000 \ # 强制1km像元尺寸单位米 -te 1800000 2000000 6200000 6300000 \ # 覆盖中国全域的Albers范围米 -co COMPRESSLZW \ tif/MOD11A2.A2022001.006.2022012075925.tif \ tif_albers/MOD11A2.A2022001.006.2022012075925_albers.tif参数说明-t_srs目标投影字符串严格匹配中国Albers标准-r near必须用最近邻LST是物理量双线性或三次卷积会生成不存在的中间值如28.37℃-tr 1000 1000明确指定1km×1km像元避免GDAL自动计算导致精度漂移-te目标范围左、下、右、上单位米此范围覆盖东起135°E、西至73°E、北至54°N、南至18°N的中国陆域——比WGS84经纬度框更紧凑减少空像元。3.2 批量重投影脚本用PythonGDAL自动化处理全年数据# batch_reproject.py import os import glob from osgeo import gdal input_dir tif/ output_dir tif_albers/ os.makedirs(output_dir, exist_okTrue) # Albers China projection string albers_srs projaea lat_125 lat_247 lon_0105 datumWGS84 unitsm no_defs # China bounding box in Albers (meters) te [1800000, 2000000, 6200000, 6300000] # xmin, ymin, xmax, ymax for tif_path in glob.glob(os.path.join(input_dir, *.tif)): if MOD11A2 in tif_path or MYD11A2 in tif_path: output_path os.path.join(output_dir, os.path.basename(tif_path).replace(.tif, _albers.tif)) # Use gdal.Warp for better control than subprocess options gdal.WarpOptions( dstSRSalbers_srs, xRes1000, yRes1000, outputBoundste, resampleAlggdal.GRA_NearestNeighbour, creationOptions[COMPRESSLZW] ) gdal.Warp(output_path, tif_path, optionsoptions) print(fReprojected {os.path.basename(tif_path)} → {os.path.basename(output_path)})运行后所有TIFF将转换为Albers投影像元中心坐标单位为米可直接与shp/china_province.shp确保其也为Albers叠加分析。4. LST值反演与质量控制把DN值转为物理温度并用QC波段剔除无效像元MODIS LST产品中每个TIFF实际包含多个波段GDAL默认只读Band 1但关键质量信息藏在附加波段里。以MOD11A2为例一个HDF文件解包后通常含Band 1LST_Day_1km白天LSTDN值Band 2QC_Day白天质量控制8-bitBand 3Emis_3131波段发射率Band 4Emis_3232波段发射率但你手上的TIFF很可能只保留了Band 1LST——这是最大隐患。若无QC波段你无法区分“真实高温”与“云污染高温”。4.1 从DN到物理LST必须执行的缩放公式官方文档MOD11A2 ATBD v6.1明确LST(K) Scale_Factor × DN Add_Offset其中Scale_Factor 0.02,Add_Offset 150.0→LST(℃) (DN × 0.02 150.0) − 273.15用Rasterio批量转换同时写入NoDataimport rasterio from rasterio.enums import Resampling import numpy as np def dn_to_celsius(tif_path, output_path): with rasterio.open(tif_path) as src: profile src.profile.copy() data src.read(1).astype(np.float32) # Read Band 1 # Apply scaling: DN - K - ℃ lst_k data * 0.02 150.0 lst_c lst_k - 273.15 # Set NoData (original DN0 becomes -273.15℃, but we want explicit NoData) lst_c[data 0] np.nan profile.update( dtyperasterio.float32, nodatanp.nan, compresslzw ) with rasterio.open(output_path, w, **profile) as dst: dst.write(lst_c, 1) # Example usage dn_to_celsius(tif_albers/MOD11A2.A2022001.006.2022012075925_albers.tif, lst_celsius/MOD11A2.A2022001.006.2022012075925_celsius.tif)4.2 QC波段缺失时的补救方案用外部云掩膜与地形约束若ZIP中无QC波段常见于简化版数据集必须引入外部约束NASA MCD12Q1土地覆被产品提取“水体”class 11和“永久冰雪”class 15像元设为NoDataSRTM DEM海拔5000m区域如喜马拉雅主峰LST信噪比极低建议屏蔽ERA5逐日云量若需高精度可下载2022年ERA5 hourly cloud cover对80%云量日期的LST整体设为NoData。简易Python实现基于土地覆被# 使用rasterio叠加土地覆被掩膜 with rasterio.open(lst_celsius/MOD11A2.A2022001.006.2022012075925_celsius.tif) as lst_src: lst lst_src.read(1) profile lst_src.profile with rasterio.open(mcd12q1_2022_china_albers.tif) as lc_src: lc lc_src.read(1) # 设水体11和冰雪15为nan lst[(lc 11) | (lc 15)] np.nan # 写回 profile.update(nodatanp.nan) with rasterio.open(lst_masked/MOD11A2.A2022001.006.2022012075925_masked.tif, w, **profile) as dst: dst.write(lst, 1)玄学提醒不要相信“云少数据好”。MODIS LST在薄卷云下仍会严重高估温度云顶辐射被误判为地表辐射因此即使QC波段显示“部分云”也建议剔除。保守做法仅保留QC_Day 0的像元二进制00000000其他全设为NoData。5. 常见问题与避坑指南那些让LST分析结果集体翻车的隐蔽陷阱5.1 现象同一省份不同日期的LST均值波动剧烈±10℃远超气候学合理范围原因未统一处理“晨昏效应”与“季节性太阳高度角”。TerraMOD上午过境AquaMYD下午过境但2022年1月1日与7月1日的太阳高度角差异巨大导致相同地物LST可差15℃以上。若混用MOD11A2与MYD11A2计算“日均LST”本质是把苹果和橙子平均。解决明确分析目标——若研究“日变化”必须配对同一天的MOD白天与MYD夜间若研究“季节变化”只用MOD11A2白天或只用MYD11A2夜间不可混用若需“日均”用气象站实测日均温校准而非简单平均。5.2 现象叠加省级矢量后LST统计值在省界处突变且河北与山东交界处出现条带状异常原因矢量边界与栅格像元未对齐。.shp文件若为WGS84而TIFF为Albers即使坐标系声明正确GDAL在rasterstats.zonal_stats中仍可能因重采样插值引入边界伪影。解决确保矢量与栅格同投影、同分辨率用ogr2ogr将china_province.shp转为Albers并用gdal_rasterize生成1km省级ID栅格统计时用exact_extract非rasterstats它基于像元中心点判断归属避免插值或改用“大像元法”将LST重采样为10km再统计牺牲精度换稳定性。5.3 现象青藏高原西部LST常年35℃与实测气象站数据15℃矛盾原因高原地区大气稀薄MODIS 31/32波段反演LST的算法Split-Window在低水汽条件下失效且发射率假设ε0.97不适用裸岩/冻土。官方文档明确标注海拔4000m区域LST不确定性3K。解决对海拔4000m像元添加警告标签或直接设为NoData引入ASTER GEDv3发射率数据替代MODIS默认发射率重运行LST反演需专业遥感软件在论文方法部分必须声明“本研究LST数据在青藏高原高海拔区存在系统性高估结论不适用于该区域”。5.4 现象导出的GeoTIFF在QGIS中显示为全黑或颜色条范围异常min-273, max150原因未设置正确的nodata值或QGIS默认拉伸为“MinMax”而LST有效值集中在20–45℃-273℃的NoData拉伸后掩盖全部细节。解决在GDAL写入时显式设置nodatanp.nan非0或-9999QGIS中右键图层→Properties→Symbology→Band Rendering→Min/Max→手动设Min10, Max50或用gdal_translate -a_nodata nan input.tif output.tif强制声明。5.5 现象用rasterio.mask按省界裁剪后文件体积暴增3倍原因裁剪生成的GeoTIFF默认不压缩且NoData区域被填充为0而非跳过。解决裁剪时传入compresslzw和predictor2针对浮点型或用gdalwarp -crop_to_cutline -co COMPRESSLZW -co PREDICTOR2替代Python裁剪。6. 进阶技巧构建可复现的LST时间序列分析流水线用DAG调度每日质检单日LST处理只是起点。真正工程价值在于构建2022全年LST时间序列并支持回溯、对比、报警。我一般用AirflowPython构建轻量DAG核心是三个原子任务6.1 任务1每日LST质检Quality Gate对每个日期的LST TIFF自动运行三项检查有效像元率 10%排除全云日像元值范围 ∈ [−20℃, 60℃]剔除缩放错误与ERA5日均温空间相关系数 0.6用scipy.stats.pearsonr计算滑动窗口10天。def quality_gate(date_str): lst_path flst_masked/MOD11A2.A{date_str}.006.*_masked.tif era5_path fera5_daily/{date_str}_t2m.tif with rasterio.open(glob.glob(lst_path)[0]) as lst, \ rasterio.open(era5_path) as era5: lst_data lst.read(1).flatten() era5_data era5.read(1).flatten() # Remove NaN mask ~np.isnan(lst_data) ~np.isnan(era5_data) r, _ pearsonr(lst_data[mask], era5_data[mask]) if r 0.6: raise ValueError(fLow correlation {r:.3f} on {date_str})6.2 任务2省级LST统计表CSV输出供BI工具接入用exact_extract替代zonal_stats确保精度exactextract -r lst_masked/MOD11A2.A2022001.006.*_masked.tif \ -p shp_albers/china_province_albers.shp \ -f csv \ -o stats_2022001.csv \ --include-geometryfalse \ --summarymean,std,min,max,count输出CSV含列province_name,mean_lst,std_lst,min_lst,max_lst,count_valid_pixels6.3 任务3热岛强度指数UHI Index计算定义UHI LST_city − LST_rural其中LST_rural取城市建成区外5km缓冲区内的平均LST。用rasterio.features.shapes提取建成区栅格再用rasterio.mask裁剪缓冲区# Step 1: Rasterize city boundary (from OSM or RESIDENTIAL landuse) # Step 2: Buffer by 5000m in Albers # Step 3: Mask LST with buffer, exclude city pixels # Step 4: Mean of remaining pixels LST_rural我的习惯绝不保存中间TIFF所有统计结果直出CSV/Parquet每天凌晨2点触发DAG邮件推送当日质检报告原始LST TIFF按日期分文件夹存档但永远不删除——因为2022年某日LST异常半年后才发现是MODIS传感器临时校准偏移需回溯原始数据重处理。这套流程跑满一年后你会拥有一份真正可信的中国1km LST时间序列——它不是“数据集”而是你亲手锻造的气候观测基础设施。希望帮到你。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

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