简介这份数据集提供2020年中国区域1km空间分辨率的地表温度LST栅格成果面向遥感、地理信息、气候与生态等方向的研究人员和学生可用于城市热岛、地表能量平衡、干旱监测及长时间序列温度变化分析等场景。数据源自NASA MODIS MOD11A2产品经提取子数据集、拼接、投影转换、单位换算与裁剪后按8天合成再平均得到年度LST采用Albers等面积投影与WGS84椭球中央经线105°、标准纬线25°与47°。压缩包共11个文件约66.16MB包含2个tif栅格开氏与摄氏两套温度值、2个tfw坐标文件、4个xml元数据、2个txt说明及1个ovr金字塔文件便于直接加载与快速渲染。目前已有3545人学习下载适合需要现成中国区年度LST底图、希望省去预处理环节并快速开展空间分析的用户。1. MODIS 2020年中国1km地表温度数据集从下载到出图的完整路径2020年夏季长江中下游持续高温很多做城市热岛、农业干旱、电力负荷分析的同行都在找同一份底图覆盖全国、1km分辨率、能直接算年均和季节均值的2020年地表温度栅格。MODIS LST 产品恰好卡在这个需求点上——它每天过境两次空间分辨率1km时间跨度从2000年至今2020年这一年的数据完整性和质量都相当好。但真正上手会发现从NASA下载的原始HDF文件到一张能放进报告里的中国区LST分布图中间隔着投影转换、无效值处理、云掩膜、拼接裁剪、单位换算五道坎。这篇笔记就按我自己跑过一遍的流程把每一步的命令、参数和翻车点讲清楚适合做遥感应用、GIS分析、气候相关研究的从业者照着复现。2. MODIS LST 产品选型与2020年数据获取2.1 MOD11A1 还是 MYD11A1先搞清楚你要哪颗星MODIS搭载在Terra和Aqua两颗卫星上地表温度产品对应两套Terra的MOD11系列和Aqua的MYD11系列。两者算法基本一致但过境时间不同——Terra地方时约10:30和22:30过境Aqua约13:30和01:30过境。做城市热岛分析午后温度更有代表性通常选Aqua的MYD11A1做干旱监测或蒸散发估算上午的Terra数据更常用。2020年中国区1km LST数据集绝大多数论文和业务系统用的是MOD11A1因为它的时间序列更长、与历史数据衔接更好。产品层级上MOD11A1是L3级日产品每天一景包含LST_Day_1km和LST_Night_1km两个主要波段以及质量控制QC波段。如果你要的是年均或月均需要自己从日产品合成如果不想做合成可以直接用MOD11A28天合成或MOD11C3月合成但空间分辨率会降到0.05度。1km分辨率这个硬指标决定了你必须从MOD11A1或MYD11A1出发。2.2 用 Earthdata 账号批量下载 2020 年全年数据下载入口是NASA的LAADS DAAC。先注册Earthdata账号然后在Find Data里搜MOD11A1时间范围选2020-01-01到2020-12-31空间范围用中国区经纬度框选大约73°E–135°E18°N–54°N。网页端一次最多下几百景全年366天、每天4个瓦片覆盖中国区手动点不现实。常见做法是用wget脚本批量拉取。# 先在LAADS DAAC生成下载链接列表保存为urls.txt # 每行格式类似 # https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MOD11A1/2020/001/MOD11A1.A2020001.h26v05.061.2021003191946.hdf # 用wget批量下载带Earthdata认证 while read url; do wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies \ --keep-session-cookies --auth-no-challengeon \ --user你的用户名 --password你的密码 \ -P ./MOD11A1_2020/ $url done urls.txt这段脚本的核心是--load-cookies和--save-cookies配合让wget在多次请求间保持Earthdata的会话状态。--auth-no-challengeon避免服务器返回401时wget反复重试。-P指定下载目录建议按年份建文件夹。参数上如果你只做中国区可以在生成urls.txt时只保留h26v05、h27v05、h28v04等覆盖中国的瓦片编号能把数据量从几十GB压到几GB。提示LAADS DAAC的下载链接有时效性生成后尽快下载。如果中断重新生成链接比续传更省事。2.3 瓦片编号与中国区覆盖范围对照MODIS采用正弦投影瓦片划分中国区主要涉及以下瓦片瓦片编号大致覆盖区域是否必选h26v05华东、华中是h27v05华南、西南是h28v04西北、青藏是h25v04华北北部是h26v04华北、东北南部是h27v04西南、藏南是h23v04新疆北部视范围h24v04新疆中部视范围实际做全国图建议把h23–h28、v04–v05全部下载避免边缘缺失。每个瓦片约10–20MB全年下来大约20–40GB取决于你保留的波段数。3. 从HDF到中国区LST栅格投影、拼接与裁剪3.1 用GDAL读取HDF子数据集并转投影MODIS HDF文件里包含多个子数据集GDAL用subdataset语法访问。先看一个文件里有什么gdalinfo MOD11A1.A2020001.h26v05.061.2021003191946.hdf输出里会列出SUBDATASET_1_NAME到SUBDATASET_12_NAME其中LST_Day_1km通常是第1个QC_Day是第3个左右。读取并转成GeoTIFF# 读取LST_Day_1km子数据集转成GeoTIFF gdal_translate -of GTiff \ HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h26v05.061.2021003191946.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km \ LST_Day_2020001_h26v05.tif # 重投影到WGS84经纬度 gdalwarp -t_srs EPSG:4326 -r bilinear \ LST_Day_2020001_h26v05.tif \ LST_Day_2020001_h26v05_wgs84.tifgdal_translate的-of GTiff指定输出格式subdataset路径必须用双引号包住否则冒号会被shell解析。gdalwarp的-t_srs EPSG:4326把正弦投影转成经纬度-r bilinear用双线性插值对连续型温度数据比最近邻更合适。注意MODIS LST的缩放因子是0.02单位开尔文转成摄氏度需要value * 0.02 - 273.15这一步在后续用Python处理。3.2 按天拼接瓦片并裁剪中国区边界每天有多个瓦片需要先拼接再裁剪。用GDAL的VRT虚拟拼接# 生成VRT拼接文件 gdalbuildvrt LST_Day_2020001.vrt \ LST_Day_2020001_h26v05_wgs84.tif \ LST_Day_2020001_h27v05_wgs84.tif \ LST_Day_2020001_h28v04_wgs84.tif \ LST_Day_2020001_h25v04_wgs84.tif # 用中国区shapefile裁剪 gdalwarp -cutline china_boundary.shp -crop_to_cutline \ -dstnodata -9999 LST_Day_2020001.vrt \ LST_Day_2020001_china.tifgdalbuildvrt不实际复制像素只建索引速度快。-cutline指定裁剪边界-crop_to_cutline让输出范围紧贴边界-dstnodata -9999设无效值。中国区边界shapefile可以从公开地理数据平台获取注意用WGS84坐标系。3.3 用Python批量处理全年数据并合成年均值逐日处理366天手动跑命令不现实。用Python脚本串起来import os import glob import numpy as np from osgeo import gdal def hdf_to_lst(hdf_path, out_path): 读取HDF中的LST_Day_1km转投影并裁剪中国区 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_Daily_1km_LST:LST_Day_1km # 转GeoTIFF temp_tif out_path.replace(.tif, _temp.tif) gdal.Translate(temp_tif, subdataset, formatGTiff) # 重投影裁剪 gdal.Warp(out_path, temp_tif, dstSRSEPSG:4326, cutlineDSNamechina_boundary.shp, cropToCutlineTrue, dstNodata-9999, resampleAlggdal.GRA_Bilinear) os.remove(temp_tif) def lst_to_celsius(tif_path): 读取LST栅格转摄氏度无效值保持-9999 ds gdal.Open(tif_path) band ds.GetRasterBand(1) arr band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() # MODIS LST缩放因子0.02单位K arr arr * 0.02 - 273.15 arr[arr -100] -9999 # 异常值过滤 return arr, ds.GetGeoTransform(), ds.GetProjection() # 批量处理 hdf_list sorted(glob.glob(./MOD11A1_2020/*.hdf)) for hdf in hdf_list: date_str os.path.basename(hdf).split(.)[1] # A2020001 out_tif f./LST_china/LST_{date_str}.tif hdf_to_lst(hdf, out_tif)gdal.Translate和gdal.Warp是GDAL Python绑定的核心函数参数与命令行一一对应。arr * 0.02 - 273.15完成开尔文到摄氏度的换算arr[arr -100] -9999过滤掉换算后明显异常的值比如原始填充值0换算后是-273.15。批量循环里按日期命名输出文件方便后续按时间索引。年均值合成用numpy的nanmean# 读取全年栅格堆叠后求年均 tif_files sorted(glob.glob(./LST_china/LST_*.tif)) stack [] for f in tif_files: arr, gt, proj lst_to_celsius(f) arr[arr -9999] np.nan stack.append(arr) stack np.array(stack) annual_mean np.nanmean(stack, axis0) annual_mean[np.isnan(annual_mean)] -9999 # 写出GeoTIFF driver gdal.GetDriverByName(GTiff) out_ds driver.Create(./LST_china_annual_2020.tif, annual_mean.shape[1], annual_mean.shape[0], 1, gdal.GDT_Float32) out_ds.SetGeoTransform(gt) out_ds.SetProjection(proj) out_ds.GetRasterBand(1).WriteArray(annual_mean) out_ds.GetRasterBand(1).SetNoDataValue(-9999) out_ds Nonenp.nanmean忽略NaN适合处理有云缺失的日数据。如果某天整景都是无效值该天不参与均值计算这是合理的。输出时把NaN还原成-9999保持GIS软件兼容。4. 质量控制与无效值处理别让云像素污染你的均值4.1 读懂QC波段哪些像素能用哪些必须扔MOD11A1的QC_Day波段是一个8位整数不同位段编码不同信息。关键位是bit 0–1LST质量和bit 2–3LST误差。常见做法是只保留QC值0高质量误差1K和1次高质量误差2K其余全部当无效值。用Python解析def mask_by_qc(lst_arr, qc_arr): 根据QC波段掩膜LST只保留QC 0和1的像素 # QC bit 0-1: 00高质量, 01次高质量, 10低质量, 11未生产 qc_quality qc_arr 0b11 mask (qc_quality 0) | (qc_quality 1) lst_masked lst_arr.copy() lst_masked[~mask] np.nan return lst_maskedqc_arr 0b11取出最低两位判断是否在0或1。这一步必须在转摄氏度之前做因为QC波段和LST波段是独立子数据集需要分别读取后对齐。4.2 云掩膜与异常值过滤的实操参数除了QC还可以用LST值本身做粗过滤。2020年中国区地表温度合理范围大约在-50°C到70°C之间超出这个范围的像素大概率是云污染或算法异常。另外水体、冰雪的LST与陆地差异大做全国均值时如果不区分下垫面结果会被水体拉低。常见做法是叠加土地覆盖数据如MCD12Q1做分层统计或者至少用NDVI阈值把水体剔掉。# 异常值过滤 lst_masked[lst_masked -50] np.nan lst_masked[lst_masked 70] np.nan # 可选用NDVI掩膜水体 # ndvi 0 通常是水体 ndvi_arr ... # 读取同期NDVI lst_masked[ndvi_arr 0] np.nan参数上-50和70是经验阈值可以根据研究区调整。青藏高原冬季LST可能低于-50如果研究区包含高原下限要放宽到-60。4.3 无效值比例统计判断某天数据能不能用处理完每天的数据后统计无效值占比。如果某天全国范围内无效值超过70%这天基本不可用建议从年均值计算中剔除。用Python快速统计valid_ratio np.sum(~np.isnan(lst_masked)) / lst_masked.size if valid_ratio 0.3: print(f{date_str} 有效像素不足30%建议剔除)2020年夏季长江流域梅雨季偏长6–7月云覆盖严重那段时间的日数据有效像素可能只有40%–50%。年均值合成时这些天的权重自然降低但如果做月均值需要额外注意。5. 避坑与排查MODIS LST 处理中最容易翻车的 5 个地方5.1 现象转投影后图像整体偏移几十公里原因MODIS HDF文件自带正弦投影参数但某些GDAL版本读取HDF4时投影信息解析不完整导致gdalwarp用了错误的源投影。解决先用gdalinfo确认源文件的投影字符串如果显示为空白或异常手动指定-s_srs projsinu R6371007.181 nadgridsnull wktext no_defs。这个正弦投影参数是MODIS的标准定义R6371007.181是球体半径。5.2 现象LST值全是0或异常大原因读取子数据集时路径写错GDAL返回了填充值而不是实际数据。MODIS HDF的填充值是0换算后是-273.15°C。解决检查subdataset路径中的MODIS_Grid_Daily_1km_LST是否与实际文件一致不同版本061 vs 006的网格名可能不同。用gdalinfo列出所有SUBDATASET确认。5.3 现象拼接后瓦片接缝处有明显条带原因不同瓦片的LST值在重叠区不一致gdalbuildvrt默认取最后一个瓦片的值导致接缝。解决用gdalwarp的-r bilinear在拼接时做插值融合或者用gdal_merge.py的-n -9999参数让无效值不覆盖有效值。更好的做法是在VRT阶段用-resolution highest统一分辨率。5.4 现象年均值图上有大片异常低温区原因冬季高纬度地区LST本身很低但QC掩膜把有效低温也滤掉了导致这些区域全是无效值年均值计算时被周围像素填充。解决检查QC掩膜阈值是否过严冬季青藏高原和东北的LST可能被标为QC2低质量如果全部剔除这些区域就没有数据。可以放宽到QC2但要在论文里说明。5.5 现象Python脚本处理到一半内存溢出原因全年366天、每天全国范围的Float32栅格堆叠后占用内存约366 × 5000 × 5000 × 4字节 ≈ 36GB。解决不要一次性读入所有天数改用增量计算——维护一个sum和count数组逐天累加最后相除。或者按季度分批处理再合成年均。# 增量计算年均值避免内存溢出 sum_arr None count_arr None for f in tif_files: arr, gt, proj lst_to_celsius(f) arr[arr -9999] np.nan if sum_arr is None: sum_arr np.zeros_like(arr) count_arr np.zeros_like(arr) valid ~np.isnan(arr) sum_arr[valid] arr[valid] count_arr[valid] 1 annual_mean np.where(count_arr 0, sum_arr / count_arr, -9999)这个写法内存占用只有两个数组适合普通工作站。6. 进阶用2020年LST做城市热岛强度分级的一个具体技巧拿到年均LST栅格后最直接的进阶用法是做城市热岛强度分级。常规做法是用城区均值减郊区均值但全国尺度下“城区”和“郊区”的界定很模糊。我自己的习惯是用MODIS土地覆盖产品MCD12Q1的IGBP分类把13类城市建成区提取出来然后在每个城市周边缓冲50km作为郊区分别统计LST均值。这样比用行政边界更客观。具体操作分三步。第一步用MCD12Q1的2020年数据提取城市类IGBP第13类重采样到1km与LST栅格对齐。第二步对每个地级市以城市类像素为中心向外做50km缓冲区缓冲区内非城市类像素作为郊区样本。第三步计算每个城市的LST均值差按差值大小分五级0°C为无热岛0–1°C为弱1–2°C为中等2–3°C为强3°C为极强。from scipy import ndimage def urban_heat_island_intensity(lst_arr, urban_mask, buffer_km50): 计算每个连通城市区域的热岛强度 # 标记连通的城市斑块 labeled, num_features ndimage.label(urban_mask) results [] for i in range(1, num_features 1): city_mask (labeled i) # 膨胀50km作为郊区范围1km分辨率下约50像素 buffer_mask ndimage.binary_dilation(city_mask, iterations50) suburb_mask buffer_mask (~city_mask) (~np.isnan(lst_arr)) city_lst np.nanmean(lst_arr[city_mask]) suburb_lst np.nanmean(lst_arr[suburb_mask]) if not np.isnan(city_lst) and not np.isnan(suburb_lst): results.append({ city_id: i, city_lst: city_lst, suburb_lst: suburb_lst, uhi_intensity: city_lst - suburb_lst }) return resultsndimage.label标记连通区域binary_dilation做缓冲区膨胀iterations50对应50km。这个方法的边界条件是大城市连片时会被标为一个斑块比如长三角城市群这时候需要按行政边界拆分或者改用滑动窗口法。分级结果可以叠加到中国地图上用分级设色出图。2020年的数据跑下来长江中下游城市群的热岛强度普遍在2–3°C珠三角部分城市超过3°C而西北城市如兰州、乌鲁木齐的热岛强度反而偏低因为周边荒漠的LST本身很高城郊温差被压缩。这个反直觉的结果提醒我做全国对比时不能只看绝对值还要考虑气候背景。最后说一个我踩过的坑一开始我用年均LST直接算热岛强度结果发现北方城市冬季热岛很强、夏季很弱年均值把季节信号抹平了。后来改成按季节分别计算才发现夏季热岛才是真正需要关注的问题。所以如果你要做热岛分析建议至少分夏季6–8月和冬季12–2月两套结果别偷懒只用年均。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/12 20:10:36 · 阅读 1.2k MODIS 2020年中国1km地表温度数据集：下载、处理与城市热岛分析全流程
简介这份数据集提供2020年中国区域1km空间分辨率的地表温度LST栅格成果面向遥感、地理信息、气候与生态等方向的研究人员和学生可用于城市热岛、地表能量平衡、干旱监测及长时间序列温度变化分析等场景。数据源自NASA MODIS MOD11A2产品经提取子数据集、拼接、投影转换、单位换算与裁剪后按8天合成再平均得到年度LST采用Albers等面积投影与WGS84椭球中央经线105°、标准纬线25°与47°。压缩包共11个文件约66.16MB包含2个tif栅格开氏与摄氏两套温度值、2个tfw坐标文件、4个xml元数据、2个txt说明及1个ovr金字塔文件便于直接加载与快速渲染。目前已有3545人学习下载适合需要现成中国区年度LST底图、希望省去预处理环节并快速开展空间分析的用户。1. MODIS 2020年中国1km地表温度数据集从下载到出图的完整路径2020年夏季长江中下游持续高温很多做城市热岛、农业干旱、电力负荷分析的同行都在找同一份底图覆盖全国、1km分辨率、能直接算年均和季节均值的2020年地表温度栅格。MODIS LST 产品恰好卡在这个需求点上——它每天过境两次空间分辨率1km时间跨度从2000年至今2020年这一年的数据完整性和质量都相当好。但真正上手会发现从NASA下载的原始HDF文件到一张能放进报告里的中国区LST分布图中间隔着投影转换、无效值处理、云掩膜、拼接裁剪、单位换算五道坎。这篇笔记就按我自己跑过一遍的流程把每一步的命令、参数和翻车点讲清楚适合做遥感应用、GIS分析、气候相关研究的从业者照着复现。2. MODIS LST 产品选型与2020年数据获取2.1 MOD11A1 还是 MYD11A1先搞清楚你要哪颗星MODIS搭载在Terra和Aqua两颗卫星上地表温度产品对应两套Terra的MOD11系列和Aqua的MYD11系列。两者算法基本一致但过境时间不同——Terra地方时约10:30和22:30过境Aqua约13:30和01:30过境。做城市热岛分析午后温度更有代表性通常选Aqua的MYD11A1做干旱监测或蒸散发估算上午的Terra数据更常用。2020年中国区1km LST数据集绝大多数论文和业务系统用的是MOD11A1因为它的时间序列更长、与历史数据衔接更好。产品层级上MOD11A1是L3级日产品每天一景包含LST_Day_1km和LST_Night_1km两个主要波段以及质量控制QC波段。如果你要的是年均或月均需要自己从日产品合成如果不想做合成可以直接用MOD11A28天合成或MOD11C3月合成但空间分辨率会降到0.05度。1km分辨率这个硬指标决定了你必须从MOD11A1或MYD11A1出发。2.2 用 Earthdata 账号批量下载 2020 年全年数据下载入口是NASA的LAADS DAAC。先注册Earthdata账号然后在Find Data里搜MOD11A1时间范围选2020-01-01到2020-12-31空间范围用中国区经纬度框选大约73°E–135°E18°N–54°N。网页端一次最多下几百景全年366天、每天4个瓦片覆盖中国区手动点不现实。常见做法是用wget脚本批量拉取。# 先在LAADS DAAC生成下载链接列表保存为urls.txt # 每行格式类似 # https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MOD11A1/2020/001/MOD11A1.A2020001.h26v05.061.2021003191946.hdf # 用wget批量下载带Earthdata认证 while read url; do wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies \ --keep-session-cookies --auth-no-challengeon \ --user你的用户名 --password你的密码 \ -P ./MOD11A1_2020/ $url done urls.txt这段脚本的核心是--load-cookies和--save-cookies配合让wget在多次请求间保持Earthdata的会话状态。--auth-no-challengeon避免服务器返回401时wget反复重试。-P指定下载目录建议按年份建文件夹。参数上如果你只做中国区可以在生成urls.txt时只保留h26v05、h27v05、h28v04等覆盖中国的瓦片编号能把数据量从几十GB压到几GB。提示LAADS DAAC的下载链接有时效性生成后尽快下载。如果中断重新生成链接比续传更省事。2.3 瓦片编号与中国区覆盖范围对照MODIS采用正弦投影瓦片划分中国区主要涉及以下瓦片瓦片编号大致覆盖区域是否必选h26v05华东、华中是h27v05华南、西南是h28v04西北、青藏是h25v04华北北部是h26v04华北、东北南部是h27v04西南、藏南是h23v04新疆北部视范围h24v04新疆中部视范围实际做全国图建议把h23–h28、v04–v05全部下载避免边缘缺失。每个瓦片约10–20MB全年下来大约20–40GB取决于你保留的波段数。3. 从HDF到中国区LST栅格投影、拼接与裁剪3.1 用GDAL读取HDF子数据集并转投影MODIS HDF文件里包含多个子数据集GDAL用subdataset语法访问。先看一个文件里有什么gdalinfo MOD11A1.A2020001.h26v05.061.2021003191946.hdf输出里会列出SUBDATASET_1_NAME到SUBDATASET_12_NAME其中LST_Day_1km通常是第1个QC_Day是第3个左右。读取并转成GeoTIFF# 读取LST_Day_1km子数据集转成GeoTIFF gdal_translate -of GTiff \ HDF4_EOS:EOS_GRID:MOD11A1.A2020001.h26v05.061.2021003191946.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km \ LST_Day_2020001_h26v05.tif # 重投影到WGS84经纬度 gdalwarp -t_srs EPSG:4326 -r bilinear \ LST_Day_2020001_h26v05.tif \ LST_Day_2020001_h26v05_wgs84.tifgdal_translate的-of GTiff指定输出格式subdataset路径必须用双引号包住否则冒号会被shell解析。gdalwarp的-t_srs EPSG:4326把正弦投影转成经纬度-r bilinear用双线性插值对连续型温度数据比最近邻更合适。注意MODIS LST的缩放因子是0.02单位开尔文转成摄氏度需要value * 0.02 - 273.15这一步在后续用Python处理。3.2 按天拼接瓦片并裁剪中国区边界每天有多个瓦片需要先拼接再裁剪。用GDAL的VRT虚拟拼接# 生成VRT拼接文件 gdalbuildvrt LST_Day_2020001.vrt \ LST_Day_2020001_h26v05_wgs84.tif \ LST_Day_2020001_h27v05_wgs84.tif \ LST_Day_2020001_h28v04_wgs84.tif \ LST_Day_2020001_h25v04_wgs84.tif # 用中国区shapefile裁剪 gdalwarp -cutline china_boundary.shp -crop_to_cutline \ -dstnodata -9999 LST_Day_2020001.vrt \ LST_Day_2020001_china.tifgdalbuildvrt不实际复制像素只建索引速度快。-cutline指定裁剪边界-crop_to_cutline让输出范围紧贴边界-dstnodata -9999设无效值。中国区边界shapefile可以从公开地理数据平台获取注意用WGS84坐标系。3.3 用Python批量处理全年数据并合成年均值逐日处理366天手动跑命令不现实。用Python脚本串起来import os import glob import numpy as np from osgeo import gdal def hdf_to_lst(hdf_path, out_path): 读取HDF中的LST_Day_1km转投影并裁剪中国区 subdataset fHDF4_EOS:EOS_GRID:{hdf_path}:MODIS_Grid_Daily_1km_LST:LST_Day_1km # 转GeoTIFF temp_tif out_path.replace(.tif, _temp.tif) gdal.Translate(temp_tif, subdataset, formatGTiff) # 重投影裁剪 gdal.Warp(out_path, temp_tif, dstSRSEPSG:4326, cutlineDSNamechina_boundary.shp, cropToCutlineTrue, dstNodata-9999, resampleAlggdal.GRA_Bilinear) os.remove(temp_tif) def lst_to_celsius(tif_path): 读取LST栅格转摄氏度无效值保持-9999 ds gdal.Open(tif_path) band ds.GetRasterBand(1) arr band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() # MODIS LST缩放因子0.02单位K arr arr * 0.02 - 273.15 arr[arr -100] -9999 # 异常值过滤 return arr, ds.GetGeoTransform(), ds.GetProjection() # 批量处理 hdf_list sorted(glob.glob(./MOD11A1_2020/*.hdf)) for hdf in hdf_list: date_str os.path.basename(hdf).split(.)[1] # A2020001 out_tif f./LST_china/LST_{date_str}.tif hdf_to_lst(hdf, out_tif)gdal.Translate和gdal.Warp是GDAL Python绑定的核心函数参数与命令行一一对应。arr * 0.02 - 273.15完成开尔文到摄氏度的换算arr[arr -100] -9999过滤掉换算后明显异常的值比如原始填充值0换算后是-273.15。批量循环里按日期命名输出文件方便后续按时间索引。年均值合成用numpy的nanmean# 读取全年栅格堆叠后求年均 tif_files sorted(glob.glob(./LST_china/LST_*.tif)) stack [] for f in tif_files: arr, gt, proj lst_to_celsius(f) arr[arr -9999] np.nan stack.append(arr) stack np.array(stack) annual_mean np.nanmean(stack, axis0) annual_mean[np.isnan(annual_mean)] -9999 # 写出GeoTIFF driver gdal.GetDriverByName(GTiff) out_ds driver.Create(./LST_china_annual_2020.tif, annual_mean.shape[1], annual_mean.shape[0], 1, gdal.GDT_Float32) out_ds.SetGeoTransform(gt) out_ds.SetProjection(proj) out_ds.GetRasterBand(1).WriteArray(annual_mean) out_ds.GetRasterBand(1).SetNoDataValue(-9999) out_ds Nonenp.nanmean忽略NaN适合处理有云缺失的日数据。如果某天整景都是无效值该天不参与均值计算这是合理的。输出时把NaN还原成-9999保持GIS软件兼容。4. 质量控制与无效值处理别让云像素污染你的均值4.1 读懂QC波段哪些像素能用哪些必须扔MOD11A1的QC_Day波段是一个8位整数不同位段编码不同信息。关键位是bit 0–1LST质量和bit 2–3LST误差。常见做法是只保留QC值0高质量误差1K和1次高质量误差2K其余全部当无效值。用Python解析def mask_by_qc(lst_arr, qc_arr): 根据QC波段掩膜LST只保留QC 0和1的像素 # QC bit 0-1: 00高质量, 01次高质量, 10低质量, 11未生产 qc_quality qc_arr 0b11 mask (qc_quality 0) | (qc_quality 1) lst_masked lst_arr.copy() lst_masked[~mask] np.nan return lst_maskedqc_arr 0b11取出最低两位判断是否在0或1。这一步必须在转摄氏度之前做因为QC波段和LST波段是独立子数据集需要分别读取后对齐。4.2 云掩膜与异常值过滤的实操参数除了QC还可以用LST值本身做粗过滤。2020年中国区地表温度合理范围大约在-50°C到70°C之间超出这个范围的像素大概率是云污染或算法异常。另外水体、冰雪的LST与陆地差异大做全国均值时如果不区分下垫面结果会被水体拉低。常见做法是叠加土地覆盖数据如MCD12Q1做分层统计或者至少用NDVI阈值把水体剔掉。# 异常值过滤 lst_masked[lst_masked -50] np.nan lst_masked[lst_masked 70] np.nan # 可选用NDVI掩膜水体 # ndvi 0 通常是水体 ndvi_arr ... # 读取同期NDVI lst_masked[ndvi_arr 0] np.nan参数上-50和70是经验阈值可以根据研究区调整。青藏高原冬季LST可能低于-50如果研究区包含高原下限要放宽到-60。4.3 无效值比例统计判断某天数据能不能用处理完每天的数据后统计无效值占比。如果某天全国范围内无效值超过70%这天基本不可用建议从年均值计算中剔除。用Python快速统计valid_ratio np.sum(~np.isnan(lst_masked)) / lst_masked.size if valid_ratio 0.3: print(f{date_str} 有效像素不足30%建议剔除)2020年夏季长江流域梅雨季偏长6–7月云覆盖严重那段时间的日数据有效像素可能只有40%–50%。年均值合成时这些天的权重自然降低但如果做月均值需要额外注意。5. 避坑与排查MODIS LST 处理中最容易翻车的 5 个地方5.1 现象转投影后图像整体偏移几十公里原因MODIS HDF文件自带正弦投影参数但某些GDAL版本读取HDF4时投影信息解析不完整导致gdalwarp用了错误的源投影。解决先用gdalinfo确认源文件的投影字符串如果显示为空白或异常手动指定-s_srs projsinu R6371007.181 nadgridsnull wktext no_defs。这个正弦投影参数是MODIS的标准定义R6371007.181是球体半径。5.2 现象LST值全是0或异常大原因读取子数据集时路径写错GDAL返回了填充值而不是实际数据。MODIS HDF的填充值是0换算后是-273.15°C。解决检查subdataset路径中的MODIS_Grid_Daily_1km_LST是否与实际文件一致不同版本061 vs 006的网格名可能不同。用gdalinfo列出所有SUBDATASET确认。5.3 现象拼接后瓦片接缝处有明显条带原因不同瓦片的LST值在重叠区不一致gdalbuildvrt默认取最后一个瓦片的值导致接缝。解决用gdalwarp的-r bilinear在拼接时做插值融合或者用gdal_merge.py的-n -9999参数让无效值不覆盖有效值。更好的做法是在VRT阶段用-resolution highest统一分辨率。5.4 现象年均值图上有大片异常低温区原因冬季高纬度地区LST本身很低但QC掩膜把有效低温也滤掉了导致这些区域全是无效值年均值计算时被周围像素填充。解决检查QC掩膜阈值是否过严冬季青藏高原和东北的LST可能被标为QC2低质量如果全部剔除这些区域就没有数据。可以放宽到QC2但要在论文里说明。5.5 现象Python脚本处理到一半内存溢出原因全年366天、每天全国范围的Float32栅格堆叠后占用内存约366 × 5000 × 5000 × 4字节 ≈ 36GB。解决不要一次性读入所有天数改用增量计算——维护一个sum和count数组逐天累加最后相除。或者按季度分批处理再合成年均。# 增量计算年均值避免内存溢出 sum_arr None count_arr None for f in tif_files: arr, gt, proj lst_to_celsius(f) arr[arr -9999] np.nan if sum_arr is None: sum_arr np.zeros_like(arr) count_arr np.zeros_like(arr) valid ~np.isnan(arr) sum_arr[valid] arr[valid] count_arr[valid] 1 annual_mean np.where(count_arr 0, sum_arr / count_arr, -9999)这个写法内存占用只有两个数组适合普通工作站。6. 进阶用2020年LST做城市热岛强度分级的一个具体技巧拿到年均LST栅格后最直接的进阶用法是做城市热岛强度分级。常规做法是用城区均值减郊区均值但全国尺度下“城区”和“郊区”的界定很模糊。我自己的习惯是用MODIS土地覆盖产品MCD12Q1的IGBP分类把13类城市建成区提取出来然后在每个城市周边缓冲50km作为郊区分别统计LST均值。这样比用行政边界更客观。具体操作分三步。第一步用MCD12Q1的2020年数据提取城市类IGBP第13类重采样到1km与LST栅格对齐。第二步对每个地级市以城市类像素为中心向外做50km缓冲区缓冲区内非城市类像素作为郊区样本。第三步计算每个城市的LST均值差按差值大小分五级0°C为无热岛0–1°C为弱1–2°C为中等2–3°C为强3°C为极强。from scipy import ndimage def urban_heat_island_intensity(lst_arr, urban_mask, buffer_km50): 计算每个连通城市区域的热岛强度 # 标记连通的城市斑块 labeled, num_features ndimage.label(urban_mask) results [] for i in range(1, num_features 1): city_mask (labeled i) # 膨胀50km作为郊区范围1km分辨率下约50像素 buffer_mask ndimage.binary_dilation(city_mask, iterations50) suburb_mask buffer_mask (~city_mask) (~np.isnan(lst_arr)) city_lst np.nanmean(lst_arr[city_mask]) suburb_lst np.nanmean(lst_arr[suburb_mask]) if not np.isnan(city_lst) and not np.isnan(suburb_lst): results.append({ city_id: i, city_lst: city_lst, suburb_lst: suburb_lst, uhi_intensity: city_lst - suburb_lst }) return resultsndimage.label标记连通区域binary_dilation做缓冲区膨胀iterations50对应50km。这个方法的边界条件是大城市连片时会被标为一个斑块比如长三角城市群这时候需要按行政边界拆分或者改用滑动窗口法。分级结果可以叠加到中国地图上用分级设色出图。2020年的数据跑下来长江中下游城市群的热岛强度普遍在2–3°C珠三角部分城市超过3°C而西北城市如兰州、乌鲁木齐的热岛强度反而偏低因为周边荒漠的LST本身很高城郊温差被压缩。这个反直觉的结果提醒我做全国对比时不能只看绝对值还要考虑气候背景。最后说一个我踩过的坑一开始我用年均LST直接算热岛强度结果发现北方城市冬季热岛很强、夏季很弱年均值把季节信号抹平了。后来改成按季节分别计算才发现夏季热岛才是真正需要关注的问题。所以如果你要做热岛分析建议至少分夏季6–8月和冬季12–2月两套结果别偷懒只用年均。希望帮到你。本文还有配套的精品资源点击获取