去年做农业气候区划客户要求把全国积温分带落到县级精度。我一开始图省事直接用了月平均气温栅格数据结果边界一细化就露馅最冷月零线跟站点实测差了将近两百公里项目负责人盯着我屏幕问“你这等值线是怎么画出来的”。被将了一军之后我才老老实实去找日尺度的产品最后换成了2000-2020年中国4km分辨率逐日2米平均温度栅格数据把积温、霜冻日数、极端温度指标全部重算了一遍才把精度兜住。这套数据说简单也简单就是给中国全境大约每4公里一个格网提供21年里每一天的2米高度平均气温。对做农业气象、生态遥感、城市化研究、健康风险区划的人来说它的价值是实打实的。这篇博文我把它从数据原理、文件读取、裁剪处理、指标计算到常见的坑完整梳理一遍适合刚接触栅格数据的GIS/遥感学生也适合已经跑过不少气温数据但没细究过逐日序列处理细节的从业者。1. 为什么说4km逐日气温栅格是“刚需”数据1.1 2米气温和地表温度是两回事先强调一个特别容易被搞混的概念。标题里写的是“2米平均温度”这是气象站百叶箱里测的那个气温代表近地面空气的热力状态。很多做过遥感的人习惯用卫星反演的地表温度LST这两个东西完全不是一回事。夏季午后裸土表面温度能到60℃以上但2米气温通常也就35℃上下反过来晴夜地表强烈辐射冷却地表温度可能比2米气温低很多。如果你做积温、作物生长模型、人体舒适度必须用2米气温用LST会得到一堆离谱结果。这套栅格数据瞄准的就是站点观测意义上的气温空间上按规则网格铺开。它把散落的气象站点观测值转换成了连续分布的栅格场这就是它作为栅格数据的核心价值你可以在任意一个经纬度位置提取温度序列而不必依赖距离最近的站点。1.2 4km分辨率在空间尺度上是“甜点位”空间分辨率不是越高越好得看数据体量和实际精度。25km/0.25°的产品跑全国尺度很快但中国地形太复杂天山、秦岭、横断山这些地方的温度空间梯度很大25km完全压不住地形信号。1km甚至500m的产品当然更细腻但全国范围逐日序列的数据量会膨胀到难以日常处理而且站点插值在1km尺度上未必真能把精度做上去有时反而会制造出虚假的细节。4km这个档位大约对应0.04°的经纬度网格在全国尺度刚好能捕捉主要山脉走向、盆地和河谷的温度分异同时数据量可控。我拿四川盆地和川西高原对比过4km数据里成都平原和邛崃山区的温度差能清晰体现出来站点点位上的验证误差也控制得不错。对中国这种多山国家来说4km是区域尺度和县级区划之间的一个很好的折中点。1.3 逐日时间尺度从“看趋势”升级到“数天数”很多公开气温产品是月尺度一个月一个值能做气候倾向率、距平分析但做不了过程类指标。农业上算积温要逐日累加生态上判断物候起止日期要精确到天极端高温、霜冻事件本身就是以“日”为单位定义的。有了逐日数据你可以自己聚合出任意时间尺度的任意指标月、季、年、多年平均都行这是月尺度产品做不到的灵活性。2000年到2020年一共21年其中2000、2004、2008、2012、2016、2020这6个闰年总天数7671天。也就是说这套数据至少包含7671个时次每个时次一张全国栅格。体量上看数据参数一般如下表所示拿到手先对照一遍避免后面对不上。项目常见参数时间范围2000-01-01 至 2020-12-31逐日空间范围中国全境及周边空间分辨率约4km0.04°左右数据内容逐日2米平均气温常见格式GeoTIFF逐日单文件或 NetCDF单文件多时次常见投影WGS84地理坐标EPSG:4326常见单位摄氏度℃少数产品可能是开尔文2. 从站点到栅格这套气温数据是怎么生产出来的2.1 台站观测与质量控制数据的源头气温栅格不是凭空算出来的它的原始输入是全国几千个国家级气象站的逐日观测记录。中国气象站网的密度分布很不均匀东部平原差不多每一两百公里有一个站青藏高原和新疆腹地则稀疏得多。做插值之前原始观测数据要过质量控制关包括气候极值检查比如某站某天温度超过该站历史极值、时间一致性检查逐日温度不会出现毫无道理的突变、空间一致性检查邻站对比发现明显孤点等。如果源头数据有错后面所有插值都是白搭。实际上很多公开数据集的算法还会把台站周边的土地利用、人口密度信息放进质量控制流程。这一步做得好不好直接决定数据集在极端天气事件的记录上准不准。我记得检查过几次极端低温事件个别数据集在寒潮过程里会出现站间温度衔接不上的现象就是质量控制环节没处理干净。2.2 空间插值用什么方法把散点变成连续面从站点到栅格核心是空间插值。常见的路线有三种。第一种是反距离加权IDW简单粗暴但完全不考虑地形山区效果很差现在已经很少用于正式产品。第二种是克里金和薄板样条它们考虑空间自相关结构能把大尺度温度场的变化趋势拟合得比较好学术圈常用的ANUSPLIN软件就是薄板样条的代表。第三种是在插值过程中加入高程协变量因为气温随海拔升高而降低是明确的物理规律平均而言每上升100米气温下降约0.6℃叫气温直减率。把DEM高程作为协变量带入插值是4km分辨率数据能做好的关键。如果没有高程项一个平坦的插值面会把太行山两侧的温度混成一锅粥有了高程约束太行山东麓和西麓就能拉出合理的温差。很多产品用的就是“薄板样条高程协变量”的组合对温度这种与地形强相关的要素非常有效。2.3 再分析产品做背景场另一种技术路线还有一类数据集不走纯插值路线而是把全球再分析资料比如ERA5-Land作为背景场再降尺度到4km并用站点观测做偏差订正。再分析资料的优点是有完整的物理模型支撑逐日序列连续且没有站点缺测的问题缺点是原始分辨率粗、地形细节不足降尺度后需要用高程、坡度这些信息把地形效应补回来。偏差订正通常按月分开做因为不同季节的系统性误差差别很大。两条路线做出来的产品各有千秋。纯插值产品在站点密集区更贴近观测但站点稀疏区容易产生不自然的边界再分析降尺度产品空间上更平滑、物理上更自洽但局部极值往往被平滑掉。你在使用某一套具体数据之前最好先看它的技术文档搞清楚它是哪个路线这样遇到异常值的时候能判断是数据问题还是真实天气过程。2.4 不确定性数据不是真理尤其在这三个地方任何栅格产品都有不确定性搞数据的人必须心里有数。第一是站点稀疏区比如青藏高原中西部、塔克拉玛干沙漠腹地站点密度极低插值结果更多依赖模型和协变量误差可能达到3-4℃。第二是复杂地形区山地峡谷里一个格网内的实际温度差异可能超过几度栅格值代表的是格网平均用它去代表一个具体的小流域或谷底要谨慎。第三是日尺度误差会被过程指标放大。这一点很多人忽略逐日数据当天的误差可能不大但你把它累加成积温、算作极端事件天数时系统性偏差会累积。比如某个格网常年平均偏低1℃算21年积温就少了几百度拿去做农业区划可能把整个种植界线画错。所以在使用任何栅格数据之前先在你有实测站点的区域做验证这一步省不得。3. 上手第一步文件读取与海量数据的组织方式3.1 GeoTIFF还是NetCDF先搞清楚自己拿到的是什么这类数据的组织方式通常有两种。一种是逐日一个GeoTIFF文件文件名带着日期21年7671个文件另一种是打包成一个NetCDF文件里面一个三维数组维度是时间、纬度、经度。两者各有侧重点GeoTIFF适合按天切片、做单日分析和传统的GIS软件操作NetCDF适合做时间序列统计因为它把整个时间维都放在一个结构里用xarray处理非常顺手。我建议你把两种格式都稍微了解因为实际使用中你会遇到混合情况有些平台下发的就是一组GeoTIFF有些则是一个NetCDF还有的会同时给。最稳妥的做法是全部统一成DataArray对象后面所有分析都基于xarray来做。3.2 用xarray读取NetCDF一行的东西别拆成一百行如果你拿到的是NetCDF处理会非常舒服。打开一个文件就能看到所有时次import xarray as xr # 打开单个NetCDF内含2000-2020全部逐日数据 ds xr.open_dataset(China_4km_daily_tavg_2000_2020.nc) print(ds) print(ds.tavg.shape) # 比如 (7671, 900, 1550) print(ds.tavg.attrs) # 查看单位、填充值说明等元数据打开之后先做三件事看维度顺序看坐标范围看变量属性。维度顺序直接决定你后面写代码时的缩放下标坐标范围能帮你快速确认经纬度边界对不对属性里通常会写清楚单位、填充值和数据说明这是最容易帮你避坑的地方。立刻检查有没有缺失值print(ds.tavg.isnull().sum())如果缺失值比例很低后面分析就省心如果某些地区常年缺失就要留意是不是数据本身覆盖不完整。3.3 GeoTIFF批量读取目录结构和命名规则要提前设计如果拿到的是逐日GeoTIFF千万不要把7671个文件全塞进一个文件夹文件管理器光列目录都要卡半天。建议按年份分目录文件名里带上标准日期格式这样排序不会乱tavg_2000/ tavg_20000101.tif tavg_20000102.tif ... tavg_2001/ tavg_20010101.tif ...用Python批量读取时最怕的就是文件名排序错乱。日期用yyyymmdd格式可以保证字符串排序等于时间排序如果用别的格式一定要用正则或日期解析排序不然时间序列就全乱了。import rasterio import numpy as np from pathlib import Path files sorted(Path(./tavg_daily).rglob(*.tif)) print(len(files)) # 如果是完整序列应该是7671 # 读取第一天 with rasterio.open(files[0]) as src: data src.read(1) profile src.profile print(data.shape, profile[crs], profile[nodata])这里的profile里面有crs和nodata信息nodata就是无效值标记后面统计时一定要处理。建议先把单日文件读通再套进循环或xarray的open_mfdataset里做时间堆叠。4. 从全国到局部裁剪与重投影实操4.1 按行政边界裁剪研究区域限定在哪数据就留哪全国数据用起来毕竟体量大大多数研究只需要某个省、某个流域或某个县域。按行政边界裁剪最方便的工具链是geopandas加rioxarray。前提是你的数据本身带着地理坐标信息NetCDF里一般有crs变量GeoTIFF则天然带投影信息。import geopandas as gpd import rioxarray ds xr.open_dataset(China_4km_daily_tavg_2000_2020.nc) # 如果坐标没有绑定CRS先显式声明 ds ds.rio.write_crs(EPSG:4326) # 读取省级边界比如四川省 shp gpd.read_file(sichuan_boundary.shp) ds_sichuan ds.rio.clip(shp.geometry, dropTrue)clip的原理是把边界外的像元全部遮挡掉保留边界内的网格。注意shp的坐标系必须和数据一致都是WGS84才不用先做投影转换。如果你要裁剪的边界是投影坐标记得先to_crs(“EPSG:4326”)。裁剪之后先看一眼格网大小确认维度和原来期望一致。4.2 经纬度范围裁剪没有边界文件的时候用这招如果你的研究区边界就是一个矩形范围比如“东经103°到105°北纬28°到32°”根本不需要边界文件直接按坐标切片就行速度极快ds_sub ds.sel( lonslice(103, 105), latslice(32, 28) # 注意纬度从大到小 )这里有个新手最容易翻车的地方纬度坐标经常是降序排列的即从高纬度到低纬度。如果你写成slice(28, 32)会得到一个空的数组或者倒序的结果。解决办法是先打印ds.lat.values看起止顺序或者用排序后的坐标再切片。养成先看坐标再切片的习惯能省掉很多莫名其妙的空数组排查时间。4.3 投影坐标系统一面积统计之前必须做的事WGS84经纬坐标的栅格在做面积计算、距离计算、或者和某些投影坐标系的矢量数据叠加时会出现问题。经纬度一度在不同纬度对应的物理距离不一样在投影坐标下直接算面积会扭曲。如果你的下游分析涉及面积统计比如算某个县的霜冻风险面积比例建议先把栅格重投影到Albers等积投影或UTM分区投影再做裁剪和统计。# 重投影到Albers等积圆锥中央经线105°E标准纬线25°N和47°N ds_alb ds.rio.reproject( EPSG:102025, res4000, # 目标分辨率4km )重投影会重采样最常用的方法是双线性内插。对于温度这种连续变量双线性没问题如果你哪天换成土地利用分类这种离散变量就得用最近邻法。这个区别一定要记住。4.4 裁剪结果的快速体检裁剪和重投影做完别急着跑分析先做一次体检看数据的空间范围是否和预期一致看有没有整块空洞看边界处有没有异常值。一个简单的办法是把裁剪结果的第一天画出来和行政边界叠加看是否对齐。这一步能发现很多隐形问题裁剪边界差了几公里导致沿海岛屿被切掉重投影时经纬度网格旋转导致边界处出现空白带或者裁剪后nodata填充值变成新值后面统计时直接污染结果。反正画一张图只要十秒值得每次都做。5. 数据分析实战从逐日气温到研究指标5.1 多年平均与气候态基准先立好“标尺”很多分析都需要一个基准期。比如做距平就得先算出某段时期的气候态均值。你可以计算2000年到2020年逐年平均温度再算21年总体平均# 计算逐年平均 annual_tavg ds.tavg.resample(time1Y).mean(time) # 计算多年平均2000-2010作为基准期 clim ds.tavg.sel(timeslice(2000, 2010)).mean(time)为什么基准期选2000-2010气象上常用WMO推荐的30年气候标准期但这套数据只有21年一般会把前期十年作基准或者直接用整个21年做平均。不同基准期会得到不同距平结果论文里要写清楚。5.2 距平与异常年识别把“哪年反常”找出来有了气候态就可以算逐年距平。最直观的做法是先聚合成月或年尺度再用多年平均做差# 逐年温度距平 annual_anomaly annual_tavg - clim # 逐月距平分月计算基准避免季节性干扰 monthly ds.tavg.resample(time1M).mean(time) monthly_clim monthly.sel(timeslice(2000, 2010)).groupby(time.month).mean(time) monthly_anomaly monthly.groupby(time.month) - monthly_clim做出来后你会发现2000-2020年间有几个明显的冷暖异常年份。2007年冬季全国普遍偏暖2010年前后西南地区出现严重干旱伴随高温2013年夏季南方极端高温2020年末的“霸王级”寒潮这些异常在距平图上都会以显著正负距平的形式出现。这其实是数据质量的一个验证手段真实发生的极端气候事件应该在数据里留下信号。5.3 GDD与积温农业上最常用的指标积温的计算逻辑不复杂以某个生物学下限温度为基准把每天高于基准的那部分温度累加起来。比如小麦常用10℃以下的界限做“无效”温度那就用日平均温度减去10℃只保留正值一年内求和就是生长度日GDD。# 计算2000-2020逐年的GDD下限10℃ daily ds.tavg - 10.0 gdd_daily daily.where(daily 0, 0) gdd_annual gdd_daily.resample(time1Y).sum(time)注意这里用了where把低于下限的日期置零而不是直接乘一个布尔掩码因为掩码乘出来的NaN会在求和时变成NaN必须用0填充。这个细节我踩过坑算出来的积温有一半是缺测就是因为当初图省事用了mask乘法。如果你需要按不同作物用不同下限温度把这个函数包起来传参就行。5.4 极端温度指标高温日数、霜冻日数逐日数据另一个大用途是统计极端事件。常见指标有高温日数日平均温或日最高温超过阈值的天数、霜冻日数日最低温低于0℃的天数、热浪持续天数等。注意标题数据是平均温度没有最高和最低所以用平均温做代理指标时要说明口径差异。# 炎热日数统计日平均温≥28℃的天数 hot_days (ds.tavg 28).resample(time1Y).sum(time) # 冷日数统计日平均温0℃的天数近似霜冻指标 cold_days (ds.tavg 0).resample(time1Y).sum(time)平均温度阈值和最高/最低温度阈值刻画的事件不完全一样。日平均温≥28℃代表整日偏热比单纯看午后最高温更严格日平均温0℃也不能完全等同于霜冻因为霜冻主要由夜间最低温决定。如果做正式研究建议同时获取最高/最低温数据如果只是做空间格局和变化趋势的快速摸底这样用也没问题。5.5 趋势分析每个格点都在变暖还是变冷把21年逐日数据聚合成年值后可以做空间趋势分析。最简单的做法是对每个格点的21个年值做线性回归斜率就是年际变化速率。更进一步可以用Mann-Kendall非参数趋势检验它对异常值不敏感是气候研究里常用的方法。输出结果可以叠加成一张趋势空间分布图标出显著的区域。这里特别提醒趋势分析结果对基准期和起止年份敏感。2000-2020这21年如果拿前五年和后五年来对比会得到不同强度的趋势。做结论时不要过度解读短时段的趋势尤其是在年际波动大的地区。6. 这套数据躲不开的坑及排查思路6.1 填充值引发的“全国都是-9999℃”事故我第一次拿这套数据算全国平均温结果西北地区全部出现了-9999左右的值乍一看像极了真实的极寒天气但直觉告诉我这不是温度是填充值漏处理了。逐日GeoTIFF文件里无效像元通常用特定数字填充常见的有-9999、-3.4e38有的文件干脆把这个值写进nodata字段有的没有。排查路径是这样的先打印单个文件的profile看nodata设置再用numpy统计最小值、最大值和唯一值数量看是否存在一个占比过大的“孤值”。确认填充值后统一替换成NaN# 统一把填充值转成NaN data_clean np.where(data -9000, np.nan, data)这个坑的教训是任何时候拿到新栅格第一步永远先看直方图。温度数据的直方图应该是单峰的、连续分布的如果看到某个值的频数异常高十有八九是填充值混进来了。6.2 闰年导致的时间对齐错位我做GDD计算时发现2008年12月31日的结果总是比前后年份少一天数据总是提前结束。排查了一圈发现不是数据缺失而是我在把GeoTIFF文件名转成时间索引时直接用了一个预设的时间序列这个序列没有包含2008年2月29日导致文件与时间错位。更隐蔽的情况恰恰相反数据集本身是按自然日排列的包含了闰年但你的脚本假设每年固定365天于是从2000年往后每过四年就错位一天到2020年已经错位5天。解决办法就是用日期字符串批量解析而不是用索引序号硬推from datetime import datetime dates [datetime.strptime(f.stem, tavg_%Y%m%d) for f in files]构建时间坐标时直接用解析出来的日期数组确保和文件一一对应。这套逻辑虽然多写两行但能避免所有闰年相关的错位问题。6.3 全国格点一次性算直接OOM如果直接用7671天乘全国格点做计算内存很容易爆。粗略算一下全国范围内约1500乘900的网格也就是135万个像元乘以7671个时次后的数据量为10的10次方量级以float32存储需要40GB左右普通电脑根本吃不消。但这不代表你没法用思路是分块处理。最推荐的是按年或按月循环for year in range(2000, 2021): ds_year ds.tavg.sel(timestr(year)) # 做该年的计算比如年均 result_year ds_year.resample(time1Y).mean(time) # 保存或累加更优雅的方案是用dask让xarray延迟计算数据会按块自动调度。但在逐日全国数据上dask的调度开销也不小我自己的经验是优先裁剪到研究区再按年循环既简单又直观内存占用能控制在2GB以内。6.4 沿海与高原的插值边缘异常数据在两类区域容易出现反常识的值一类是海岸带另一类是高大山体边缘。沿海地区有时会出现比同一纬度内陆明显偏高或偏低的孤立像元这往往是插值处理海陆边界时的伪影高原边缘则可能出现一天之内温度剧烈跳变的情况那是插值面在陡峭地形转折处的过冲。应对方法不复杂把异常像元筛出来后和周边站点对比确认是否合理。如果确认是数据伪影直接掩膜掉或者用周边像元均值替换。同时也要注意到这些区域本身观测困难真实日温差大别一概当错误处理。判断标准是看空间连续性单个孤立像元异常大概率是数据问题成片且和地形走势一致的异常则可能是真实信号。6.5 开尔文和摄氏度的深夜翻车单位问题听起来低级但我见过不止一个组在这上面栽过。有些再分析衍生数据集的温度变量单位是开尔文而国内很多脚本默认按摄氏度处理。如果你直接把开尔文温度当摄氏度算积温GDD会大得离谱因为全年的温度都高于“下限”10度。解决方法就一句话拿到数据先打印attrs养成习惯。要是发现单位是K直接减273.15ds[tavg_c] ds.tavg - 273.15我自己现在只要接到新数据第一件事就是看单位而不是看空间范围这个习惯救了好几次。7. 这套温度栅格在真实项目中的典型玩法7.1 农业领域霜冻风险区划和作物种植适宜性逐日气温序列最直接的应用就是农业。拿4km数据做春季最后一次霜冻日期的空间分布可以给经济林果的防冻提供依据计算生长季长度和积温能划定不同熟性玉米、水稻的适宜种植区。比如某个县想评估猕猴桃能否安全越冬用这套数据提取该县极端最低日平均温的历史最小值再叠加地形和坡向分析就能给出比较稳妥的风险区划。农业保险定价也常用这套数据。保险公司需要知道不同区域的极端低温频次和强度逐日栅格数据正好能提供空间上连续的风险概率估计比用散点站数据做费率分区细致得多。7.2 城市方向热岛效应评估与高温健康预警城市热岛研究通常需要城区和郊区的温度对比。有了逐日栅格可以按土地利用数据把城区像元和周边乡村像元分别聚合计算热岛强度的年际变化和夏季高温日数的城乡差异。用这套数据配合人口格网还能估算暴露在高温环境下的人口规模这对健康城市规划和高温预警启动标准的制定很有用。需要注意的仍然是分辨率问题4km在城市尺度只能反映城市群层面的热岛特征对单条街道和单个街区的微气候无能为力。如果你研究的是街区尺度的热环境需要更高的数据源4km适合做区域战略层面的判断。7.3 生态与物候温度是很多生态过程的“总开关”植被物候的起止时间、春季返青的早晚、森林上线的位置都直接受温度控制。把逐日温度数据和NDVI时间序列叠加可以识别春季温度累计量与植被变绿日期之间的关系。比如暖冬之后春季物候提前寒潮之后物候推迟这些过程用逐日温度数据都能定量刻画。再往深一层把逐日温度输入到生态过程模型里可以模拟土壤呼吸、净初级生产力的季节变化。虽然4km分辨率对单点生态过程来说还是粗但对全国尺度的生态区划和碳收支估算这个尺度恰好合适。7.4 叠加其他数据做综合场景分析气温栅格很少单独出场。实际项目里常见的组合包括和降水数据叠加算干旱指数、和风速湿度数据叠加算人体舒适度、和土地利用数据叠加算城市扩张的温度效应、和人口格网叠加算气候风险暴露度。逐日2米气温作为基础变量配合其他要素的空间数据能组合出大量有价值的派生指标。做数据叠加时要先统一空间分辨率。如果降水数据是10km或25km要上采样或下采样到同一网格保持坐标系一致。这一步处理不好后期所有归因分析都可能是错的。用这套数据快三年了我个人的习惯是任何新数据到手先拿一个小范围、短时段做全流程验证再推向全国尺度。具体做法是先选中一个自己熟悉的省份把某一年365天的数据完整跑一遍从读取、裁剪、统计到绘图全部走通确认所有参数没问题再开始处理全部21年。因为全国尺度的逐日分析跑一次往往要几小时一旦中间参数错了重跑的时间成本太高。同时把原始文件的元数据、处理脚本、版本号完整保存下来数据体量大的时候更要严格管理版本。逐日栅格数据的分辨率越高对细节的要求就越苛刻但处理得当之后它带来的分析深度是月尺度数据完全给不了的。