简介面向遥感、地理信息科学研究与国土空间分析人员这份2020年中国1km分辨率地表温度LST数据集源自美国NASA定期发布的MODIS MOD11A2 8天合成产品经提取子数据集、拼接、投影栅格、单位换算与裁剪再对全年栅格做平均处理可免去自行下载分块数据并拼接的流程。数据覆盖中国全境时间分辨率为年空间分辨率为1km投影采用Albers等积投影与WGS84椭球可直接用于生态监测、气候评估与热环境制图。压缩包共11个文件约66MB核心为开氏与摄氏两个tif栅格另有xml元数据、tfw坐标文件、txt说明和ovr金字塔文件分别对应元数据、地理配准、文档资料与大范围快速显示功能便于在ArcGIS、QGIS等软件中读取与缩放。同时提供kelvin、celsius两种温标并附原始MOD11A2引用信息与处理步骤说明便于在论文或报告中准确标注数据来源。目前已有3545人学习下载适合遥感、地理信息研究人员及规划分析人员直接使用全国尺度LST栅格数据。1. 2020年中国1km地表温度数据集一张能打开但要用对的栅格一个做城市热岛研究的同事从共享平台下载了“MODIS 2020年中国1km地表温度LST空间分布数据集”ArcGIS里打开灰蒙蒙一片直接算全国年平均温度得到一堆-9999。这个数据集听起来“已经是成品”实际底层是MODIS的MOD11A2/MYD11A2 8天合成产品白天和夜间分开、单位为开尔文且带0.02缩放因子只有把单位换算、QC质量控制和投影处理对才能用于区域统计与制图。它适合做区域气候、城市热岛、生态与农业遥感的中尺度分析不适合做严格日尺度的单点对比。下面从数据底层逻辑讲起把下载到能用的链路拆开。2. 数据底子MOD11A2与MYD11A2的8天合成、1km分辨率怎么理解2.1 MODIS LST的“1km”是星下点标称分辨率不等于每个像元都是正南正北的1km×1kmMODIS热红外波段第31、32波段的星下点空间分辨率约为1km官方LST产品基于分裂窗算法反演。到了产品网格里标称“1km”但实际像元大小与扫描角有关卫星在边缘扫描时瞬时视场会被拉大产品分发前按正弦投影做了重采样所以严格说不是每一点刚好1km。这一点在做面积统计时影响不大但在与高分辨率土地利用数据叠加时要注意尺度匹配不要按像元数直接折算面积。另外要特别注意原始MODIS陆地产品包括MOD11A2/MYD11A2使用正弦投影Sinusoidal。你在共享平台下载到处理过的中国区栅格通常是别人已经重投影成WGS84经纬度或Albers等积投影如果拿到的是原始HDF则要先gdalinfo确认子数据集和投影。我的经验是拿到任何一份LST数据先把投影和缩放因子写进处理笔记否则半年后回来看图单位是摄氏度还是开尔文都会搞混。2.2 8天合成把云遮掉的时间“匀”到一起不是连续8天的完整平均MOD11A2对应Terra星上午约10:30过境MYD11A2对应Aqua星下午约13:30过境。热红外无法穿透云层单日产品在中国南方经常出现大片无值区。为了得到较完整的空间覆盖MODIS把一年按8天一段分成约46个合成期每个合成期从所有晴空观测中生成白天的LST值和夜间的LST值。所以它不是“连续8天的平均气温”而是“8天内所有合格晴空观测的综合”。晴天多时这个值更接近晴空地表温度阴雨天多时它代表的是少数晴空窗口的采样。这会带来一个隐含问题合成期内没有云覆盖的地区才参与计算云覆盖多的地方参与天数少。后续做月、季、年统计时如果不做有效观测天数加权统计结果会偏向晴空多的时段。这正是“下载的MODIS地表温度数据能不能直接用”的根源能用于相对分布和晴空地表温度研究不能当成完全连续的时间序列或气象站气温直接用。2.3 为什么做2020年中国区域年尺度分析这版数据集比单日产品更顺手单日LST在中国东部经常缺值要做全国拼接需要处理大量空洞2020年版作为一个已经裁好边界的年序列省掉了从NASA批量下载数百景HDF再拼接的重复劳动。同时MOD11A2/MYD11A2的算法比较成熟同类产品里引用量最大做文献对比时有共同基准。适合用它做的事全国或省级尺度的年均、季节均LST制图城市热岛强度估算地表温度与NDVI/EVI的散点分析长时间序列的趋势初筛。不适合的事研究具体某一天的热浪过程或要求逐日连续输出的生态模型输入——那些场合请换MOD11A1逐日产品或再分析气温数据。做决定的标准不是“哪个文件更小”而是“时间粒度能不能对上研究问题”。3. 把2020年LST数据“拿”到本地GEE导出和本地裁剪两条路线3.1 路线AGEE导出全国年合成图附JS代码如果只是要一张全国年均LST栅格或要在省界、流域里做统计GEE是最快的路线不用下载原始HDF直接在云端处理。核心思路是构建2020年1月1日到12月31日的ImageCollection把每个8天合成期换算成摄氏度再用mean()聚合。下面是一段可以直接改的GEE JavaScript代码// 中国边界替换成你自己的矢量asset var china ee.FeatureCollection(users/your_name/china_border); // MOD11A2 白天地表温度2020年全年 var mod11 ee.ImageCollection(MODIS/061/MOD11A2) .filterBounds(china) .filterDate(2020-01-01, 2021-01-01) .select([LST_Day_1km, QC_Day]); // 把DN换算成摄氏度并用QC最低两位保留质量合格的像元 function toCelsius(img) { var lst img.select(LST_Day_1km).multiply(0.02).subtract(273.15) .rename(LST_C); var qc img.select(QC_Day); var good qc.bitwiseAnd(3).eq(0); // 最低两位为00表示质量良好 return lst.updateMask(good); } var lstC mod11.map(toCelsius); // 全年平均每个像元只对有效值求均值有空值则忽略 var annualMean lstC.mean().clip(china); // 导出到云盘坐标系用WGS84像元大小1000米 Export.image.toDrive({ image: annualMean, description: China_LST_Day_2020_annual_mean, scale: 1000, crs: EPSG:4326, maxPixels: 1e13 });代码逻辑说明select([LST_Day_1km, QC_Day])把温度和质检信息同时取到toCelsius里先乘0.02得到开尔文再减273.15得到摄氏度随后用QC最低两位筛选像元。MODIS的DN是整数单位是0.02K如果不乘比例因子后面所有统计都会差在一个“倍率”上不做QC掩膜云污染像元会混进年均值。参数说明scale1000表示导出栅格像元大小1000m与标题里的1km一致crsEPSG:4326是经纬度坐标导出后若要算面积请再重投影到Albers等积投影maxPixels要改大因为全国范围在1km尺度下仍有上千万个有效像元。如果想同时做夜间把波段换成LST_Night_1km与QC_Night再跑一套即可。3.2 路线B本地处理原始HDF附gdalwarp批量命令如果你的后续流程必须和本地土地利用、气象观测数据一起跑或者对云端环境不放心可以走本地路线。第一步先拿一个HDF文件查看子数据集名称gdalinfo MOD11A2.A2020001.h23v04.061.2020*.hdf | grep -A2 Subdatasetgdalinfo输出的Subdataset是类似HDF4_EOS:EOS_GRID:...:MODIS_Grid_8Day_1km_LST:LST_Day_1km的字符串。这个网格名不要拍脑袋写不同版本MODIS产品有差异以输出为准。然后对中国范围做投影转换和裁剪gdalwarp -overwrite \ -t_srs EPSG:4326 \ -te 73 18 135 54 \ -tr 0.01 0.01 \ -r bilinear \ HDF4_EOS:EOS_GRID:MOD11A2.A2020001.h23v04.061.2020*.hdf:MODIS_Grid_8Day_1km_LST:LST_Day_1km \ LST_day_2020001.tif这里的-te指定了中国边界的经纬度范围西到东、南到北-tr 0.01度约为1km-r bilinear是双线性重采样。MODIS原始LST的无效值通常为0处理成GeoTIFF后建议顺手设一个NoData否则后续会看见大量像元为0°C把平均温度拉低。批量处理一年46个周期要写循环常见做法是在shell里对文件名做展开mkdir -p out for hdf in MOD11A2.A2020*.hdf; do base$(basename $hdf .hdf) gdalwarp -overwrite -t_srs EPSG:4326 -te 73 18 135 54 \ -tr 0.01 0.01 -r bilinear \ HDF4_EOS:EOS_GRID:${hdf}:MODIS_Grid_8Day_1km_LST:LST_Day_1km \ out/${base}_lst.tif done参数说明这个循环假设你已经把2020年Terra星所有瓦片下载好放在当前目录如果做夜间把LST_Day_1km换成LST_Night_1km。“h23v04”这类编号是正弦投影的瓦片索引中国区域通常跨多个瓦片下载前要提前算好覆盖范围否则会缺一块。3.3 两条路线怎么选先看后续分析在哪个环境里完成路线适用情况最怕的坑GEE全国统计、时间序列探索、快速出图边界asset权限、maxPixels不够、导出后仍用4326投影算面积本地与自有栅格/样本点做pipeline、批量训练样本HDF子数据集名搞错、瓦片覆盖不全、NoData设置遗漏我的习惯探索阶段用GEE发现数据质量问题和统计口径后再用本地路线把需要的波段固化成GeoTIFF归档。不要让原始HDF直接散落在硬盘里文件名里加上“比例因子已乘/未乘”的标注这是避免返工的关键。4. 从DN到摄氏度缩放因子、QC位掩码与重投影参数4.1 缩放因子0.02为什么处理代码里都写着乘0.02MODIS LST科学数据集的存储类型是16位整数单位是开尔文官方元数据给的比例因子是0.02偏移量是0。也就是说文件里的数值乘0.02才是开尔文温度换算到摄氏度再减273.15。常见错误是只乘0.02忘记减273.15或者把它当成摄氏度直接用。下面是一段本地tif的Python换算import rasterio import numpy as np with rasterio.open(LST_Day_1km.tif) as src: dn src.read(1) nodata src.nodata profile src.profile # DN - 开尔文 - 摄氏度无效值保留为NaN kelvin dn.astype(np.float32) * 0.02 celsius kelvin - 273.15 if nodata is not None: celsius[dn nodata] np.nan # 很多MODIS产品的无效值直接是0也需要排除 celsius[dn 0] np.nan profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(LST_Day_1km_degC.tif, w, **profile) as dst: dst.write(celsius, 1)逻辑说明先转float32再乘0.02是因为16位整数乘以小数会丢精度减273.15放在乘完之后顺序不能反。无效值判断要在换算后统一转NaN这样后续统计用的numpy mean或rasterstats都会自动跳过NaN。参数说明src.nodata是读取数据自带的无值标记有的文件是-9999有的是0。MODIS原始HDF转出来的GeoTIFFNoData经常没写进文件所以代码里额外加了celsius[dn 0] np.nan兜底。如果你拿到的文件把NoData设为-9999而DN是无符号整数dn nodata会因类型不匹配不生效务必先看元数据。4.2 QC位掩码不是“QC为0就能用”这么简单QC_Day是逐像元的质量标识每个字节的不同bit表示不同信息。最低两位bit0-1代表LST生产质量00表示数据质量好01表示其他质量问题10表示云污染11表示未生产。最低两位为00的像元才是定量统计的候选。其余bit还包含云覆盖标识和平均发射率误差做严格要求时应按MOD11用户手册逐位解析。在Python里提取“最低两位为00”的掩膜import rasterio import numpy as np with rasterio.open(LST_with_QC.tif) as src: qc src.read(2) # 假设QC_Day是第二个波段 quality_mask (qc 3) 0 print(质量良好像元占比, quality_mask.mean())参数说明qc 3是位与操作取最低两位 0表示值是00。在GEE里对应qc.bitwiseAnd(3).eq(0)。实际工作中常有人直接把QC等于0当成好值这不完全对QC为0只表示质量合理其他bit可能标记云边界等最好结合业务要求决定比如做城市热岛时可以把QC为1的边界像元也放进来做严格反演验证则只保留00。4.3 区域均值不能直接算在经纬度网格上中国横跨约60度经度WGS84经纬度网格上的像元面积随纬度变化同样一个像元在黑龙江和海南对应的真实面积差不少。若直接用EPSG:4326的栅格做省级zonal statistics均值会被低纬地区面积权重轻微影响通常数值影响不大但做严格面积加权或总量估算时不可接受。常见做法是重投影到Albers等积投影。中国区域常用双标准纬线Albersgdalwarp -overwrite \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 \ -tr 1000 1000 -r bilinear \ LST_Day_1km_degC.tif LST_Day_1km_albers.tif参数解释lat_1、lat_2是两条标准纬线lon_0105°E大致处于中国中间使全中国变形最小-tr 1000 1000把像元重采样成等积投影下真实的1000米后续统计面积时每个像元面积一致。双线性重采样适合连续变量若用最邻近地表温度在像元边缘会出现台阶状突变。区域统计用rasterstatsfrom rasterstats import zonal_stats stats zonal_stats(provinces.shp, LST_Day_1km_albers.tif, stats[mean, count, min, max])这里的mean是投影后按像元等权计算count表示参与计算的有效像元数。做年均温度时同时输出count它反映数据覆盖如果某个省count很少就要怀疑云污染严重不能把它和别的省直接比较。5. 避坑下载的地表温度数据不能直接用的5个典型现场5.1 现场一把8天合成当成了“8天的平均温度”现象我用2020年46期LST做折线图发现时间序列呈明显阶梯状相邻两期温度差能到5℃。 原因8天合成值是合成期内晴空观测的合成结果不是8天的连续平均每一期的观测样本数量和时刻完全不同合成周期边界处自然出现台阶。 解决先把时间戳对齐成合成周期的中间日起始日4天再看趋势如果研究要求逐日连续改用MOD11A1逐日产品。做月均值时建议用每期有效像元数做加权平均而不是简单平均。5.2 现场二图层打开黑漆漆一片或统计出来全是-9999现象把tif拖进ArcGIS没有显示拉伸后要么全黑要么全红统计均值出现-9999。 原因通常是没有处理NoData也没有应用0.02缩放因子。很多分发版栅格把无效值设为-9999或0软件把NoData当普通数值参与显示拉伸动态范围被拉坏。 解决在图层属性里把NoData设置为文件元数据中的无效值勾选“忽略NoData”要做计算时用前面第4章代码把无效值置为NaN。显示问题不算数据问题但它最容易吓退新手。5.3 现场三和气象站气温对比RMSE大得离谱现象用LST和全国气象站2m气温对比白天RMSE到8℃以上怀疑数据分辨率或投影做错了。 原因LST是地表辐射温度2m气温是空气温度物理量本就不同晴空白天地表温度通常远高于气温夜间可能低于气温。MODIS过境是某个时刻站点观测是整点气温时间也对不齐。 解决不要做逐日点对点对比。可以比较“白天晴空LST”与同类地表温度产品或做时间窗口的统计对比在论文里务必写清楚比较的是物理意义不同的温度并标明过境时刻。5.4 现场四文件名里的日期不是观测日期现象我要取2020年1月1日的温度直接翻阅MOD11A2.A2020001开头的文件以为它就是1月1日的值。 原因MODIS 8天合成产品的日期标记是合成周期起始日A2020001对应周期是1月1日到1月8日。 解决做时间序列前把时间戳改成起始日期加合成周期长度的一半例如第1期改为1月4日。这会显著减少“相位错位”对趋势估计很重要。5.5 现场五算年均温时空洞多的地方结果很“虚”现象省级年均温度排序里青藏高原某些县比四川盆地还高明显不合理。 原因热红外LST在云多的季节大量缺测年均值默认忽略NoData像元实际统计的只是“晴空期的样本平均”某个区域冬天几乎没有合格像元时年均值其实偏向春夏晴空期。 解决统计均值的同时统计有效像元数有效像元数少于一定阈值比如全年46期里不足20期的区域只做参考。可以把白天和夜间分开统计或用每期有效天数做加权避免直接用单一年均值。这些现场合在一起说明MODIS LST数据是“下载后要做很多处理才能用”的遥感产品。能不能直接用理解8天合成、经过QC筛选、单位换算和投影转换后它能用于区域和季节尺度的地表温度分析直接拉伸显示、直接和气温对比、直接算年均那基本都要翻车。6. 用站点观测给LST做交叉验证一个可复用的最小脚本6.1 提取栅格像元并与观测对比拿到处理好的2020年一景或一期LST不要急着跑全流程。我的习惯是先选10—20个站点做交叉验证。下面这段Python用rasterstats从站点坐标提取像元值和观测序列一起计算RMSEfrom rasterstats import point_query import numpy as np # 站点经纬度示例换成你自己的 lons [116.4, 121.4, 108.9] lats [39.9, 31.2, 34.3] # 同一时段的观测地表温度℃从气象站或野外仪器读取 obs [23.5, 28.1, 26.8] # 从LST栅格提取像元值双线性插值避免像元边缘抖动 lst_vals point_query( zip(lons, lats), LST_Day_1km_degC.tif, interpolatebilinear ) lst_vals np.array([v if v is not None else np.nan for v in lst_vals]) valid ~np.isnan(lst_vals) rmse float(np.sqrt(np.mean((lst_vals[valid] - obs[valid]) ** 2))) print(RMSE(℃):, rmse)参数说明interpolatebilinear表示像元值采用双线性插值站点落在多个像元边界时结果更平滑如果站点对应像元被云污染point_query会返回None代码里统一转成NaN避免参与RMSE计算。对比时站点观测和MODIS过境时刻尽量对齐Terra约上午10:30Aqua约下午13:30站点数据最好取最近整点。6.2 我平时怎么判断一份LST数据能不能进入后续建模先抽一个8天周期做三件事打印QC掩膜的有效像元比例画一张QC良好像元的空间分布图将LST与同区域的气温做一次散点。如果有效比例低于30%或散点完全没有相关性我不会继续批量处理而是先排查单位、投影、版本。这个“先抽一景验证再全量跑批”的习惯帮我避开了很多无效投入LST的坑多在数据制备阶段不在算法阶段早一点把QC和缩放因子处理对后面才是真正的分析工作。希望帮到你。本文还有配套的精品资源点击获取