做气候变化、农业区划、生态模型这类项目的人手头最缺的往往不是算法而是一份能直接用、精度过得去的长时间序列气温数据。2005到2025年整整21年每天一张全国范围的平均气温栅格图这就是一套能支撑很多研究的基础底图。我是做GIS和气象数据处理出身的这几年帮好几个团队搭建过类似的数据集今天就把从数据源处理、空间插值到质量检验的完整思路和实操流程拆开讲。这篇东西适合搞气候分析的研究生、做灾害评估的工程师也包括那些已经买了站点观测数据、但不知道怎么能转成栅格产品的人按照文中的流程走一遍基本能做出可用的逐日平均气温栅格数据集。1. 数据需求与方案设计先把用途搞明白再做技术选型1.1 站点数据到栅格数据本质是从“点”到“面”的跨越气象站观测得到的是单个点的气温数值但气候模型、作物模型、生态模型需要的往往是空间连续的气温场。栅格数据把研究区划分为规则网格每个像元存储一个气温值这样才能和土地利用、NDVI、人口密度等其它栅格数据做逐像元运算也才能按流域、行政区、任意多边形去做统计汇总。我经常打一个比方站点数据像是十几个温度计挂在房间不同角落栅格数据则是给整个房间铺了一张连续的“温度地图”。没有这张地图你只能回答“某某站今天多少度”有了它才能回答“某个区域平均多少度、哪里出现了低温冷害”。2005到2025年这个时间窗口选得很有讲究21年跨度的逐日数据既能用来算气候态平均值也能分析极端事件的频率变化还能支撑农业积温、采暖制冷度日等应用。相比只做典型年份比如2008年雪灾、2013年高温长序列逐日数据最大的优势在于可以灵活聚合——你要月平均就从逐日平均你要生长期积温就按日累加统计口径完全在自己手里。1.2 日均温的统计口径直接决定数据能不能复用做逐日栅格之前必须先把“日平均气温”的定义讲清楚。国内气象业务上通常用一天中02时、08时、14时、20时四次定时观测的算术平均值但很多科研用户拿到的资料里只有日最高和日最低气温这时会采用简化的二时次平均也就是最高气温与最低气温的平均值Tmean(TmaxTmin)/2。这两种口径的结果并不完全等价四次观测平均更接近真实日均温而最高最低平均在锋面过境、强对流天气下误差会明显偏大山区夏季午后甚至能差出2到3摄氏度。所以方案设计的第一步不是选插值方法而是定口径、写文档。我的做法是无论最终用哪种口径都会在数据说明文件里明确标注并且在文件名里体现比如使用最高最低平均的可以叫Tmaxmin使用四次平均的可以叫T4obs。否则数据传递两轮之后使用者根本不知道手里的日均温是怎么算出来的后续做趋势分析或者积温计算都会埋下隐患。2. 数据源选型与预处理地基不牢后面全是白搭2.1 数据源选择观测资料、再分析资料如何取舍构建逐日气温栅格数据源大致有三类国家气象站逐日观测数据、区域自动站加密观测数据、再分析资料如ERA5。国家级站点空间分布相对均匀但密度不高全国大约2400多个站西部稀疏、东部密集用于全国尺度插值时精度基本够用区域自动站密度高但历史一致性差2010年前的站点数量远少于现在而且观测仪器、维护水平参差不齐早期数据有明显噪声。再分析资料如ERA5的好处是空间连续、时间序列完整但它是模式同化产物近地面气温在复杂地形、极端天气事件中的偏差可达4摄氏度以上直接拿来当“真值”用容易出问题。我的建议是长时序气候网格产品以国家级站观测数据为骨架区域自动站只做验证和局地修正。不要不同来源混着插否则站点密度在时间上剧烈变化插值结果会出现虚假的年际波动。站点的密度不均匀问题也要正视青藏高原、新疆南部站点稀疏插值结果在这些区域的不确定性就是大这不是换插值方法能解决的需要在成果报告中写明“哪些区域可信、哪些区域仅供参考”。2.2 质量控制处理缺测、异常值、台站迁移必须逐一排查站点观测数据从来不是拿来就能用的质量控制是整套流程里最花时间、也最容易被低估的一步。我处理过不少数据集第一个遇到的问题就是缺测。单站连续缺测不超过5天的可以用相邻站回归或气候值插补连续缺测超过5天的我倾向于在插值当天直接剔除该站而不是强行用统计方法填补——长时间缺测后“脑补”出来的值会以虚假的空间平滑掩盖真实气候事件。异常值检测我通常用双保险第一道是百分位法计算该站历史同期比如过去21年的同一天前后各15天的0.01至99.99百分位超出范围的标记为可疑第二道是空间一致性检测将待检站值与周围100公里内站点的值做对比如果偏差超过了同期空间标准差的三倍就标记为异常交给人工判断。台站迁移的问题也很隐蔽2005到2025期间不少台站搬迁过迁站后海拔、周围环境都变了气温序列会出现系统性台阶。查询台站元数据、标记断点、必要时分段处理这一步偷懒后面插值出来的空间分布就会在断点年份出现一圈古怪的“突变”。3. 空间插值核心细节如何在复杂地形下把精度挤出来3.1 插值方法对比IDW、克里金与回归残差法空间插值方法的选择直接决定栅格数据在站点稀少区域的表现。反距离权重IDW原理简单、计算快但它有一个天生的毛病——“牛眼效应”也就是站点周围会出现同心圆状的数值伪影在站点分布不均时会非常难看。普通克里金通过半变异函数建模空间自相关性能给出插值误差估计理论上比IDW更优雅但逐日数据量巨大每天重新拟合半变异函数很耗时而且克里金对非平稳过程比如气温随海拔的确定性变化并不敏感如果研究区横跨青藏高原和华北平原直接普通克里金插值会把山地和平原混为一谈。我自己最推荐的是“多元回归残差插值”的组合方案先用气温与经纬度、海拔建立回归模型拟合出大尺度气候趋势面这部分抓住了“随纬度升高变冷、随海拔升高变冷”的确定性规律再对回归残差做IDW或克里金插值这部分抓住局地小尺度波动最后两者相加。这个方案在复杂地形区的表现实测下来比单纯的IDW或克里金有明显优势而且计算量可控。简单说就是把“确定性趋势”和“随机剩余”分开处理各有各的方法而不是一锅烩。3.2 气温直减率复杂地形插值的胜负手全国尺度逐日气温插值最难的不是算法本身而是如何把地形的影响“塞”进去。对流层大气中气温随海拔升高而降低平均直减率大约是0.6摄氏度每100米但这个值并非固定不变——夏季、午后、湿润条件下直减率偏小冬季、夜间、干燥条件下偏大山区甚至会出现逆温层。用固定直减率做校正夏季山区午后误差能到5摄氏度以上。实操中我一般这样处理把全年分组每组内建立气温与经纬度、海拔的逐步回归方程把海拔系数作为“当日直减率的经验估计”。比如按候5天分组一年73组每组拟合一次回归就能捕捉直减率的季节变化。然后对残差插值再到高分辨率DEM上逐像元使用当日回归方程计算趋势面。这样一来每个像元的气温都包含了海拔校正山脉走向、河谷盆地的局部气候特征都能呈现出来。这个方法比固定直减率略麻烦但精度提升非常明显尤其是做积温、低温冷害这类对温度绝对值敏感的应用值得多花这一步。3.3 分辨率和坐标系统的选择分辨率不是越高越好。全国范围的逐日栅格如果做到1公里分辨率数据量约为960万个像元Float32单精度存储一天大约38MB21年累计超过280GB如果做到5公里一天大约6MB21年总计大约50GB。对大多数气候统计和区域评估研究来说5公里分辨率已经足够但如果要结合高分辨率土地利用做精细化农业区划1公里版本会更合适。我的建议是交付两种规格一套5公里用于快速分析和长期趋势统计一套1公里用于精细化应用两套数据用同一条处理流程生成保证一致性。坐标系统方面全国范围我推荐使用Albers等积圆锥投影或Lambert等角圆锥投影做面积计算和空间统计发布时同时提供WGS84地理坐标的GeoTIFF版本方便在GIS软件里直接叠加其它数据。有一个细节容易被忽略栅格数据必须确保像元对齐不同时期生成的两张栅格如果投影参数或像元左上角坐标有微小差异后续逐像元运算时就会出现错位表现为时间序列上明显的“锯齿”。我在批处理中会把投影参数写死在一个配置文件里每次生成都从同一个模板复制地理变换信息这样就不会出现对齐问题。4. 批处理工作流实操从原始站点数据到逐日栅格4.1 全流程步骤拆解整个批处理流程我分成六个环节每个环节都有明确的输入输出方便中途检查数据清洗标准化统一站点编号、日期格式、字段命名把各来源数据整理成统一CSV或Parquet格式。计算日平均气温按既定口径计算逐日平均温同时保留最高、最低气温字段备用。分组拟合回归趋势面按候或月分组对每组建气温与经纬度、海拔的回归方程。残差空间插值对每个站点当日观测值与回归趋势面的差值做IDW或克里金插值生成残差栅格。叠加生成最终栅格用回归方程在DEM上计算趋势面栅格加上残差栅格得到当日的最终气温栅格。质量检查与归档检查极值范围、空间分布、与站点实测值的差异生成质量报告再归档保存。图示这个流程其实就像做一道菜清洗是备菜回归趋势面是调底味残差插值是加 garnish最后叠加起锅。每一步都有独立的中间产物哪一步出了问题都能回溯定位。4.2 Python实现思路与关键代码我主要用Python做整套流程核心库是GDAL、NumPy和PyKrige插值前的数据管理用Pandas。下面给一个残差插值加趋势面叠加的核心代码片段抛砖引玉import numpy as np import pandas as pd from osgeo import gdal from scipy.interpolate import griddata def interpolate_daily_tmaxmin(stations_df, dem_path, out_path, date_str): stations_df: 包含 lon, lat, elev, tmean 的DataFrame dem_path: 与目标栅格同投影的DEM路径 out_path: 输出GeoTIFF路径 # 读取DEM ds gdal.Open(dem_path) geotransform ds.GetGeoTransform() elev ds.ReadAsArray().astype(np.float32) rows, cols elev.shape # 网格坐标 lon0, lat0 geotransform[0], geotransform[3] xres, yres geotransform[1], geotransform[5] lon_grid lon0 (np.arange(cols) 0.5) * xres lat_grid lat0 (np.arange(rows) 0.5) * yres lon_mesh, lat_mesh np.meshgrid(lon_grid, lat_grid) # 1. 拟合回归趋势面tmean ~ lon lat elev X np.column_stack([stations_df[lon], stations_df[lat], stations_df[elev]]) y stations_df[tmean].values X_design np.column_stack([np.ones(len(X)), X]) coef, _, _, _ np.linalg.lstsq(X_design, y, rcondNone) # 2. 计算残差 fitted X_design coef residual y - fitted # 3. 残差空间插值IDW也可以用克里金 points stations_df[[lon, lat]].values # 注意全国范围用投影坐标更好这里简化为经纬度 grid_res griddata(points, residual, (lon_mesh, lat_mesh), methodlinear) # 4. 计算趋势面并叠加残差 trend coef[0] coef[1]*lon_mesh coef[2]*lat_mesh coef[3]*elev result trend grid_res # 5. 写出GeoTIFF driver gdal.GetDriverByName(GTiff) out_ds driver.Create(out_path, cols, rows, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(geotransform) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(result) out_band.SetNoDataValue(-9999.0) out_ds.FlushCache() print(f{date_str} 完成) # 示意调用 # for date_str in date_list: # stations_df load_stations(date_str) # interpolate_daily_tmaxmin(stations_df, dem_path, out_path, date_str)这段代码有几个细节需要注意。第一是网格坐标的构建用(np.arange(cols)0.5)*xres是取像元中心而不是像元边界这个细节很多人会忽略但它直接影响栅格与站点的空间匹配关系。第二是lstsq求解回归系数用最小二乘一次算完比用sklearn的LinearRegression更省内存在逐日循环里能快不少。第三是残差插值用scipy的griddata方法选linear它在边界处会外插产生NaN后续要用最近邻或者填-9999处理。4.3 大数据量下的并行与存储方案2005到2025年大约7660天逐日插值如果单线程跑每个文件算它3到5秒一天天跑下来也要七八个小时加上回归拟合和IO整体耗时差不多一两天。实际生产环境中我一般用Python的multiprocessing或者直接上分布式调度把7660天的任务按年份拆给多个进程16核机器一晚上就能跑完一年半载的数据。注意并行时要小心临时文件夹的写入冲突给每个进程分配独立的中间目录最后再由一个汇总任务做质量检查和归档。存储层面逐日GeoTIFF文件数量会达到七千多个管理起来要规范目录结构。我的习惯是year/month/date_tmaxmin_1km.tif这样分层并同时生成NetCDF格式的集合文件把所有天的数据打包到一个NetCDF里维度为(时间, 纬度, 经度)方便用xarray直接做时间维度的切片和统计。NetCDF还有一个好处是能压缩存储一日气温数据用NetCDF压缩后只有GeoTIFF的1/3到1/5适合长期保存。5. 数据质量检验与常见问题排查5.1 交叉验证怎么做才靠谱栅格数据生成之后不能看一眼“颜色好看”就完事必须做定量验证。最常用的方法是交叉验证把站点分成训练集和验证集用训练集建模型插值在验证集站点上比较预测值与实测值。我一般用十折交叉验证计算平均绝对误差MAE、均方根误差RMSE和决定系数R2。对这个全国逐日气温数据集理想情况下MAE在1摄氏度以内、RMSE在1.5摄氏度以内复杂地形区如横断山区和站点稀疏区如藏北高原的误差会明显偏大这要在数据文档里单独说明。除了数值指标空间分布合理性检查也很重要。我会随机抽查几十天把栅格结果和站点实测值叠加显示重点看三点一是等温线是否沿山脉走向、海岸线走向延伸二是盆地、谷地是否出现明显的局地高温或低温中心三是高海拔山区是否比周边明显偏冷。这几年检查中我发现最常见的问题出在“牛眼效应”上尤其在站点密集但数据质量差的地区解决方法是改用克里金或增加趋势面回归的惩罚减弱孤立点的影响。5.2 时间序列一致性最容易被忽略的坑逐日数据最怕的是时间不一致。我们做的是21年长序列任何某个时间段内的插值参数、站点数量、质量控制的改变都可能在年际对比中制造虚假信号。比如2010年后区域自动站大量增加如果把这些站混入插值2010年后的插值结果会因为站点密度上升而出现“空间图案变化”容易被误读为气候变化。所以我在设计流程时会固定一个“站点集合”——即使是区域自动站参与插值也只使用一套在2005到2025年间始终保持观测的固定站点这样至少在空间代表性上是稳定的。但对于国家级站点来说这部分站点密度已经固定用它们做长序列分析不用担心站点变迁带来的伪变化。另外检查逐年平均值、逐年极端值是否平滑过渡、没有异常跳变这个步骤必不可少。5.3 数据使用中的常见问题速查整理几个使用这套数据时最常遇到的问题做成了速查表问题现象可能原因处理建议河谷、盆地区域栅格值异常偏高站点多位于河谷插值时高密度站点主导趋势面在地形相对高程大的区域使用回归残差法并加大海拔权重局部等温线呈同心圆状IDW插值产生的牛眼效应改用克里金或对站点做聚类抽稀后再插值同一天不同版本数据值不同日平均气温统计口径不同核对数据文件名和元数据统一口径后再比较时间序列上某年出现整体跳变站点集合变化或台站迁移检查该年份站点元数据和插值参数是否变更与站点实测值相比栅格值在山区偏差大站点稀疏、直减率估计不准使用当日分组回归的直减率而不是固定值栅格间配准错位投影参数或像元对齐不一致统一从同一模板复制地理变换参数这些坑我前前后后都踩过尤其是牛眼效应和口径混乱耗费了不少时间。好的数据集必须自带元数据文档把处理口径、站点来源、误差情况讲清楚这份文档的价值有时比数据本身还高。6. 栅格数据的通用扩展应用场景6.1 与其它栅格产品的叠加分析既然做出来的是标准栅格数据它和市面上常见的栅格数据产品天然兼容比如热搜词里提到的“全国城市形态栅格数据集”、地下水位栅格数据这类产品本质都是同一套空间数据范式。气温栅格和城市形态栅格叠加就能分析城市热岛强度与建成区密度的关系和地下水位栅格叠加就能评估气候变暖对地下水补给的影响。栅格数据的最大优势就在于格式统一、坐标对齐后可以直接做逐像元运算而不需要复杂的空间匹配和拓扑运算。这类分析我做过不少主要经验是数据在分析之前先检查坐标系、分辨率、范围是否一致不一致的用重采样或裁剪统一到同一网格然后再进入计算。很多人拿着气温是WGS84、城市形态是Albers投影的数据就直接做栅格计算器结果因为投影不一致导致结果整体偏移了几十公里这是很冤的错误。6.2 从逐日到多时间尺度产品的聚合逐日栅格数据最灵活的地方在于能够按需聚合。可以聚合成逐月平均、逐季平均、逐年平均做气候态分析也可以计算年积温、极端高温日数、霜冻日数等指标还可以提取任意时间段的距平和异常做气象灾害监测。我自己的经验是做聚合时先做好“按像元累加”的中间结果再除以有效天数这样能处理缺测像元的问题——在某个像元上如果某几天没有值累加时按有效天数计数最后除以计数即可避免把缺测当成0值算进平均。对绝大多数研究者来说拿到标准逐日气温栅格相当于获得了一个可以自由加工的原材料而不是一个固定不可变的成品。这也是我坚持把整套方法论写出来的原因——授人以鱼不如授人以渔有了流程和数据后续的所有扩展分析都成了顺理成章的填空题。7. 一点经验之谈数据集的构建流程说完了最后聊一点跟技术无关、但又比技术更磨人的体会。做2005到2025年这个跨度的时间序列栅格数据真正的难点从来不是插值算法或者Python代码而是耐心——耐心地处理每一个异常值、耐心地核对每一个年份的站位变化、耐心地给每一步操作写文档。我见过太多人花了大力气把数据跑出来却因为不愿意写元数据说明最后自己半年后再看都不确定当时处理的细节。另外一个很深的体会是先做小范围验证再铺开到全国。刚开始做的时候我直接在Python里循环跑全年结果跑到第五天发现趋势面系数异常返工浪费了半天时间。后来改成先挑两个月做全流程验证把误差、耗时、存储量全部确认一遍再批量跑剩下来的7600多天效率反而高出一大截。做这类工程化的数据产品稳定压倒一切慢就是快。最后一句话送给大家做气候变化空间分析数据质量永远大于分析方法。一套经过严格质量控制、口径清晰、文档完备的气温栅格数据集哪怕只用简单的趋势分析也能得出可信的结论反过来即使算法再花哨垃圾进垃圾出结果也毫无意义。希望这篇思路和流程的梳理能帮到正在跟气温数据较劲的你。