ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

ERA5逐小时数据聚合为日数据的三种方法:CDO、NCL与Python实战指南

ERA5逐小时数据聚合为日数据的三种方法:CDO、NCL与Python实战指南 1. 项目概述从小时到天数据聚合的必经之路处理气象再分析数据尤其是像ERA5这样的高分辨率全球数据集是很多气象、气候乃至水文、生态领域研究者的日常。当你从Copernicus Climate Data StoreCDS辛辛苦苦下载了TB级别的逐小时数据后面临的第一个现实问题往往就是如何高效、准确地将这些高频数据聚合成更易管理的逐日数据无论是为了节省存储空间、匹配其他日尺度数据还是直接进行气候态分析这个“降尺度”在时间维度上的操作都是绕不开的一步。我自己在几年前刚开始用ERA5做极端降水分析时就曾在这个环节上踩过不少坑。比如直接用编程循环去读取成千上万个NetCDF文件做日平均结果内存爆了又比如没搞清楚累积量和瞬时量的区别把“总降水”的日总和算成了小时平均值的24倍闹了笑话。所以今天我想结合自己的实战经验系统聊聊从ERA5逐小时数据获取逐日数据的三种主流方法使用气候数据运算符CDO、利用NCAR命令语言NCL以及用Pythonxarray进行灵活处理。这三种方法各有优劣适合不同的工作流和需求场景我会把它们的核心命令、背后的计算逻辑、我踩过的坑以及如何避坑都详细拆解给你。无论你是刚接触气象数据处理的新手还是想优化现有流程的老手这篇文章都能给你提供可直接“抄作业”的方案和值得注意的细节。我们的目标很明确安全、准确、高效地把原始的小时数据变成我们科研或业务中真正可用的日数据。2. 核心概念与数据准备理解你的“原料”在动手之前我们必须先搞清楚要处理的是什么“食材”。ERA5数据变量繁多但大致可以分为两类这对后续的聚合操作至关重要。2.1 ERA5变量类型辨析瞬时量与累积量这是整个处理流程的基石一旦搞错结果全盘皆输。瞬时量指在某个特定时刻如UTC 00:00, 01:00...的大气状态。例如2米气温 (2t)海平面气压 (msl)10米风场的U/V分量 (10u,10v)相对湿度 (r)对于瞬时量获取日数据通常意味着计算日平均值。即将一天内24个时刻的值相加后除以24。这代表了这一天气象要素的平均状态。累积量指从前一个时刻到当前时刻的累积总量。在ERA5中这类数据通常带有“_accumulated”的描述。最典型的就是总降水 (tp)地表净太阳辐射 (ssr)地表净热辐射 (str)对于累积量获取日数据意味着计算日总和。更精确地说由于ERA5的逐小时累积量给出的是从“time-1 hour”到“time”这一小时的累积值因此一天UTC 00:00到次日00:00的总量应该是当日01:00到次日00:00这24个小时的累积值之和。这里有一个关键陷阱如果你下载的数据是“逐小时”的那么每个文件中的tp已经是小时累积量。但如果你下载的是“月度”数据其中可能包含一个名为tp的日总量变量那就不要再重复累加了。注意在CDS下载时务必确认你选择的“时间聚合”类型。对于需要日总量的变量应选择“逐小时”数据来自行聚合这样灵活性最高。2.2 数据下载与结构检查假设我们已经从CDS下载了2023年1月全球的逐小时2米气温(2t)和总降水(tp)数据。文件可能被分割成多个NetCDF例如era5_2t_tp_20230101-20230115.nc和era5_2t_tp_20230116-20230131.nc。在开始聚合前强烈建议先用轻量级工具检查一下文件结构# 使用ncdump查看文件头信息了解变量和维度 ncdump -h era5_2t_tp_20230101-20230115.nc # 或者使用CDO的sinfo命令信息更规整 cdo sinfo era5_2t_tp_20230101-20230115.nc这个步骤能帮你确认时间维度的单位是否是hours since ...、变量名、是否有缺失值、数据的填充方式等。我曾经遇到过时间坐标不连续因下载失败导致某小时数据缺失的情况提前检查可以避免后续聚合出错。3. 方法一使用CDOClimate Data Operators—— 瑞士军刀法CDO是气象海洋领域数据处理的“瑞士军刀”它通过命令行提供了一系列强大的统计和操作功能。其最大优点是高效、内存友好因为它通常是流式处理数据无需将整个数据集加载到内存中。3.1 基础聚合命令与原理对于瞬时量如2t求日平均cdo daymean input_hourly.nc output_daily_mean.nc这个命令背后CDO会按照时间维度识别出同一天的所有时间步通常为24个。对每个网格点分别计算这24个值的算术平均值。输出一个新的NetCDF文件时间维度从“小时”变为“天”。对于累积量如tp求日总和cdo daysum input_hourly.nc output_daily_sum.nc同理daysum操作会对每个网格点将一天内的24个小时累积值进行求和得到该日的总降水量。3.2 多文件处理与批量操作实战我们很少只处理一个文件。更常见的场景是处理数个月甚至数年的数据它们被分成多个NetCDF文件。方案A使用cat或mergetime合并后再聚合这是最直观的方法但适用于文件数量不多、总大小可控的情况。# 首先按时间顺序合并所有小时文件 cdo mergetime era5_*.nc merged_hourly_all.nc # 然后进行日聚合 cdo daymean merged_hourly_all.nc final_daily_mean.nc缺点合并成一个巨大的文件可能占用大量磁盘空间且如果中间步骤出错需要重头再来。方案B流式管道聚合推荐这是更优雅和高效的方式尤其适合处理大量文件。CDO可以直接读取文件列表进行聚合。# 一次性完成按时间合并并计算日平均 cdo daymean -mergetime era5_20230101-20230115.nc era5_20230116-20230131.nc output_daily.nc # 或者使用通配符但要确保文件顺序正确按文件名排序 cdo daymean -mergetime era5_*.nc output_daily.nc这里的-符号是CDO的管道符表示将前一个操作mergetime的输出直接作为后一个操作daymean的输入不生成中间文件。方案C循环处理再合并如果你的数据是按月存放并且希望保持月度文件的组织方式可以采用循环for month in {01..12}; do cdo daymean era5_2023${month}_hourly.nc era5_2023${month}_daily.nc done # 最后将12个日文件合并 cdo mergetime era5_2023*_daily.nc era5_2023_daily.nc3.3 CDO实战心得与避坑指南时间坐标的陷阱确保你的输入文件时间坐标是连续的且没有重复。使用cdo showdate input.nc可以快速检查。如果存在缺失daymean会基于实际存在的小时数计算平均可能导致错误。建议先用cdo setmissval,1e20 input.nc fixed.nc处理缺测值。内存与磁盘空间流式操作-能极大节省磁盘IO。对于超大型数据可以结合-f nc4选项指定输出NetCDF4格式它支持压缩如-z zip_6能显著减少输出文件体积。变量选择如果文件中有多个变量但你只想处理其中一个可以使用selvar命令先选择变量避免不必要的IO。cdo daymean -selvar,tp merged_hourly_all.nc daily_tp_sum.nc性能调优在拥有多核的服务器上可以为CDO编译OpenMP支持并通过设置环境变量export OMP_NUM_THREADS4来加速部分计算密集型操作。4. 方法二使用NCLNCAR Command Language—— 传统脚本法NCL是一门为气象数据可视化与分析而生的解释型语言其数组语法和丰富的内置函数对处理网格数据非常友好。虽然NCAR已宣布对其进入维护模式并推荐转向Python但大量遗留脚本和其强大的文件读写能力使其仍是许多业务系统或老牌科研团队的选择。4.1 NCL脚本编写核心步骤下面是一个完整的NCL脚本示例用于计算2米气温的日平均和总降水的日总和。load $NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl begin ; 1. 设置输入输出文件路径 fname_in era5_2t_tp_20230101-20230115.nc fname_out era5_daily_20230101-20230115.nc ; 2. 打开文件读取变量和时间 fin addfile(fname_in, r) t2m_hour fin-t2m ; 假设变量名是t2m请根据实际情况修改 tp_hour fin-tp time_hour fin-time lat fin-latitude lon fin-longitude ; 3. 获取时间信息并创建日维度 utc_date cd_calendar(time_hour, 0) ; 将时间转换为年月日时分秒数组 year tointeger(utc_date(:,0)) month tointeger(utc_date(:,1)) day tointeger(utc_date(:,2)) ; 构建一个唯一的“日标签”例如20230101 date_str sprinti(%0.4i, year) sprinti(%0.2i, month) sprinti(%0.2i, day) uniq_dates get_unique_values(date_str) ; 需要 contributed.ncl 中的这个函数 ndays dimsizes(uniq_dates) ; 4. 预定义日尺度数组 dims dimsizes(t2m_hour) nlat dims(1) nlon dims(2) t2m_daily new((/ndays, nlat, nlon/), typeof(t2m_hour), t2m_hour_FillValue) tp_daily new((/ndays, nlat, nlon/), typeof(tp_hour), tp_hour_FillValue) copy_VarAtts(t2m_hour, t2m_daily) ; 复制属性 copy_VarAtts(tp_hour, tp_daily) t2m_daily!0 time t2m_daily!1 latitude t2m_daily!2 longitude tp_daily!0 time tp_daily!1 latitude tp_daily!2 longitude ; 5. 核心循环按天聚合 do d 0, ndays-1 ; 找到属于这一天的所有小时索引 day_indices ind(date_str .eq. uniq_dates(d)) if(.not. any(ismissing(day_indices))) then ; 计算日平均气温 t2m_daily(d, :, :) dim_avg_n_Wrap(t2m_hour(day_indices, :, :), 0) ; 计算日总降水 tp_daily(d, :, :) dim_sum_n_Wrap(tp_hour(day_indices, :, :), 0) end if delete(day_indices) ; 及时清理内存 end do ; 6. 创建新的时间坐标通常取每日的00:00或12:00作为代表 ; 这里简化处理取每天第一个小时的时间作为日时间坐标 daily_time_ind new(ndays, integer) do d0, ndays-1 daily_time_ind(d) min(ind(date_str .eq. uniq_dates(d))) end do time_daily time_hour(daily_time_ind) time_dailyunits time_hourunits time_dailycalendar time_hourcalendar ; 7. 关联维度坐标 t2m_dailytime time_daily t2m_dailylatitude lat t2m_dailylongitude lon tp_dailytime time_daily tp_dailylatitude lat tp_dailylongitude lon ; 8. 写入新文件 system(rm -f fname_out) ; 强制删除旧文件 fout addfile(fname_out, c) ; 定义文件维度 dim_names (/time, latitude, longitude/) dim_sizes (/ndays, nlat, nlon/) dim_unlim (/True, False, False/) filedimdef(fout, dim_names, dim_sizes, dim_unlim) ; 定义变量 filevardef(fout, time, typeof(time_daily), time) filevardef(fout, latitude, typeof(lat), latitude) filevardef(fout, longitude, typeof(lon), longitude) filevardef(fout, t2m_daily, typeof(t2m_daily), dim_names) filevardef(fout, tp_daily, typeof(tp_daily), dim_names) ; 复制变量属性 filevarattdef(fout, time, time_daily) filevarattdef(fout, latitude, lat) filevarattdef(fout, longitude, lon) filevarattdef(fout, t2m_daily, t2m_daily) filevarattdef(fout, tp_daily, tp_daily) ; 写入数据 fout-time time_daily fout-latitude lat fout-longitude lon fout-t2m_daily t2m_daily fout-tp_daily tp_daily print(处理完成输出文件 fname_out) end4.2 NCL方法优缺点与适用场景优点控制力极强你可以精确控制聚合的每一个逻辑例如处理非24小时数据如3小时数据、自定义聚合函数如日最高/最低温、或处理不规则时间序列。内存管理直观通过分块读取或索引操作可以处理大于内存的数据集。与图形化无缝衔接聚合后的数据可以立即用NCL强大的绘图功能进行可视化检查。缺点代码冗长相比于CDO的一行命令NCL需要编写数十行脚本。学习曲线需要熟悉NCL的语法、数组索引和文件读写API。生态衰退官方支持减弱新项目更推荐Python。适用场景当你的聚合逻辑非常复杂例如需要根据另一变量的阈值来有条件地聚合或者你的整个数据处理分析流程已经建立在NCL生态中时使用NCL是合理的选择。5. 方法三使用Pythonxarray dask—— 现代灵活法这是目前科研和业务中增长最快的方法。xarray库完美地封装了NetCDF数据模型而dask库则提供了并行计算能力两者结合使得在Python中处理大型气象数据变得既直观又高效。5.1 基于xarray的核心聚合操作首先确保安装环境pip install xarray netcdf4 dask基础聚合脚本示例import xarray as xr import numpy as np # 1. 打开数据集 # 使用open_mfdataset可以轻松打开多个文件并自动按时间合并 ds_hourly xr.open_mfdataset(era5_*.nc, combineby_coords, chunks{time: 240}) # 使用chunks参数启用dask惰性加载这里将时间维度分块每块240个时间步 # 2. 查看数据 print(ds_hourly) print(ds_hourly[t2m].attrs) # 查看变量属性确认是瞬时量还是累积量 # 3. 执行聚合操作 # 对于瞬时量如t2m计算日平均 ds_daily_mean ds_hourly[t2m].resample(time1D).mean(dimtime) # 对于累积量如tp计算日总和 ds_daily_sum ds_hourly[tp].resample(time1D).sum(dimtime) # 4. 将结果合并到一个新的Dataset中 ds_daily xr.Dataset({t2m_daily: ds_daily_mean, tp_daily: ds_daily_sum}) # 5. 写入新文件 # 设置编码进行压缩节省空间 encoding {var: {zlib: True, complevel: 5} for var in ds_daily.data_vars} ds_daily.to_netcdf(era5_daily_output.nc, encodingencoding) # 6. 关闭文件句柄使用with语句更佳 ds_hourly.close()resample(1D)是核心方法它先将数据按“1天”的频率重新分组然后接.mean()或.sum()进行聚合。xarray会自动处理时间坐标。5.2 利用Dask处理超大规模数据当数据量远超单机内存时Dask的优势就体现出来了。上面的open_mfdataset中已经使用了chunks参数这表示数据没有被立即加载到内存而是被分成了多个块chunk。所有的聚合操作resample,mean,sum都是惰性计算只有在执行compute()或to_netcdf()时才会真正触发计算。# 惰性计算示例构建复杂的处理链 lazy_daily_mean ds_hourly[t2m].resample(time1D).mean() lazy_daily_sum ds_hourly[tp].resample(time1D).sum() # 此时没有实际计算发生 print(type(lazy_daily_mean)) # 输出class xarray.core.dataarray.DataArray # 方案A直接写入文件触发计算 lazy_daily_mean.to_netcdf(t2m_daily.nc, computeTrue) # 方案B显式计算并保存到内存适用于后续还需多次操作的结果 daily_mean_computed lazy_daily_mean.compute() # 现在daily_mean_computed是一个普通的numpy数组在内存中你可以通过配置Dask的分布式调度器将计算任务分发到多核CPU甚至计算集群上从而大幅提升处理TB级数据的速度。5.3 Python方法进阶技巧与常见问题时间轴对齐问题ERA5数据的时间坐标通常是“hours since 1900-01-01”。xarray在resample时能很好地理解这一点。但如果你发现聚合后的时间点不是你期望的比如日数据的时间戳是00:00还是12:00可以使用resample(time1D, labelleft)或labelright来控制标签位置使用closedleft或right来控制区间闭合端。处理缺失值resample操作默认会跳过全为NaN的组。如果你希望保留NaN组需要设置skipnaFalse。但要注意对于累积量求和如果某天有部分小时数据缺失skipnaFalse会导致该日总和为NaN。自定义聚合函数resample后不仅可以接mean,sum还可以接min,max,std或者使用apply方法传入自定义函数灵活性极高。# 计算日最高温和最低温 daily_tmax ds_hourly[t2m].resample(time1D).max() daily_tmin ds_hourly[t2m].resample(time1D).min() # 自定义函数计算日温差 def daily_range(da): return da.max() - da.min() daily_trange ds_hourly[t2m].resample(time1D).apply(daily_range)性能优化chunks的大小设置是关键。通常一个chunk的大小应适合内存如100MB-1GB。对于时间序列按时间分块是高效的。可以通过ds_hourly.chunks查看当前分块情况并使用rechunk()方法调整。6. 三种方法对比与选型建议为了更直观地对比我将三种方法的核心特点总结如下表特性维度CDO (气候数据运算符)NCL (NCAR命令语言)Python (xarray dask)上手速度极快一行命令完成核心操作较慢需学习语法和API中等需基础Python和库知识处理效率非常高流式处理内存占用低中等取决于脚本优化程度高结合Dask可并行但内存管理需注意功能灵活性中等受限于预定义运算符极高可编写任意复杂逻辑极高可无缝集成Python生态大数据支持优秀原生支持流式处理需手动分块处理优秀原生集成Dask并行框架代码可读性命令行简洁但功能隐式脚本较长过程式编程脚本清晰声明式编程易于理解生态与未来稳定气象领域专用工具维护中不推荐新项目活跃主流趋势社区庞大调试便利性较差错误信息有时晦涩中等可逐行调试好可使用Python调试工具典型适用场景常规、批量的日/月/年统计复杂、非标准的聚合需求或遗留系统现代数据分析流程需与机器学习、可视化等深度集成选型建议如果你是新手或只想快速、稳定地完成标准的日聚合任务无脑选择CDO。它的学习成本最低效率最高一行命令就能得到可靠结果。把时间花在分析结果上而不是数据处理上。如果你的聚合逻辑非常特殊例如只聚合降水大于0.1mm的那些小时或需要复杂的时空条件判断或者你所在的团队/项目已有成熟的NCL代码库可以考虑使用NCL。它给你最大的控制权。如果你正在进行一项新的研究处理的数据量巨大TB级或者后续的分析、可视化、机器学习都计划在Python生态中进行那么Python (xarray)是你的不二之选。它代表了当前地球科学数据计算的主流方向灵活性和扩展性无可比拟。从长期来看投资学习xarray的回报率最高。在实际工作中我经常混合使用这些工具。比如用CDO进行数据的初步裁剪、格式转换和快速检查用xarraydask在Jupyter Notebook里进行交互式的探索性分析和复杂计算而一些特定的、已成型的分析图表可能还是会用NCL脚本来批量生成。工具是为人服务的选择最适合当前任务和自身技能栈的那一个。7. 常见问题排查与实战技巧实录即使知道了方法在实际操作中还是会遇到各种“坑”。下面是我和同事们总结的一些典型问题及解决方法。7.1 时间处理相关错误问题1聚合后时间坐标错乱或不是整日。现象日数据文件里的时间显示为“2023-01-01 12:00:00”而不是“2023-01-01 00:00:00”。原因与解决CDOCDO的daymean默认输出时间戳是当天所有时间点的平均值。如果你希望代表日的是00:00可以使用shifttime命令调整cdo shifttime,-12hour daily_mean.nc daily_mean_00UTC.nc。xarray使用resample(time1D, labelleft, closedleft)。labelleft表示使用区间左边界00:00作为标签closedleft表示区间是左闭右开[00:00, 次日00:00)这符合大多数气象日界的定义。核心理解你需要的“日”是日历日00-24 UTC还是业务日如12-12 UTC并根据需求调整聚合逻辑和时间标签。问题2遇到闰秒或不连续时间轴导致聚合失败。现象CDO报错“Grids have different size”或xarray报错“Cannot resample”。解决首先用cdo showtimestamp或xr.decode_cf()检查时间坐标是否连续、等间隔。对于ERA5通常很规整。如果数据源本身有问题可能需要先用cdo setgridtime或xarray的reindex进行时间轴修复。7.2 数据精度与单位转换问题降水单位从“米”到“毫米”的转换。说明ERA5中的总降水(tp)单位是“米”这是国际单位制。但气象学中常用“毫米/天”。操作在聚合之后进行单位转换。CDO:cdo mulc,1000 daily_tp_sum_m.nc daily_tp_sum_mm.ncxarray:ds_daily[tp_daily] ds_daily[tp_daily] * 1000切记同时更新变量的units属性例如tp_daily.attrs[units] mm day-1。7.3 内存不足与性能优化问题处理全球多年数据时内存爆炸OOM。CDO方案CDO本身流式处理内存压力小。但如果文件太多mergetime可能产生巨大的中间文件。最佳实践是使用管道或分时段处理。# 不好的做法先合并全年数据巨大文件再聚合 cdo mergetime *.nc huge.nc cdo daymean huge.nc out.nc # 好的做法流式管道一次处理一个月 for year in {2020..2023}; do for month in {01..12}; do cdo daymean -mergetime era5_${year}${month}*.nc daily_${year}${month}.nc done done cdo mergetime daily_*.nc final_daily.ncxarrayDask方案这是其优势所在。关键在于合理设置chunks。# 错误的打开方式直接加载到内存 ds xr.open_mfdataset(era5_*.nc) # 如果文件很大立刻OOM # 正确的打开方式使用分块和惰性加载 ds xr.open_mfdataset(era5_*.nc, chunks{time: 240, latitude: 100, longitude: 100}) # 分块原则使每个块的大小在10MB-100MB左右便于在内存中操作。 # 可以通过 ds.nbytes / 1e6 估算数据总大小通过 ds.chunks 查看分块情况。如果单机内存依然不足可以考虑使用Dask的分布式调度器连接到一个计算集群。7.4 结果验证如何确保你的日数据是对的这是最后也是最关键的一步。不要假设程序运行没报错结果就是对的。抽样检查选择一个你知道天气过程的区域和时间点。例如你知道2023年7月20日北京下了一场大雨。用Panoply、ncview或Pythonmatplotlib快速可视化该日的降水场看空间分布是否合理。提取该格点的小时降水序列手动计算日总和与程序输出的日值对比。气候态检查计算一段时间如一个月的日数据空间平均绘制时间序列。看看是否出现异常的跳变或恒值这可能是缺测值处理不当或单位错误。与已知数据源对比如果可能将你处理得到的日平均气温、日总降水与公开的ERA5日数据产品如ERA5-Land或站点观测数据进行粗略对比检查量级和变化是否一致。使用CDO的diff命令如果你用两种方法如CDO和Python处理了同一份数据可以用cdo diff比较两个输出文件看是否在机器精度内一致。cdo diff method1_daily.nc method2_daily.nc数据处理是一个需要耐心和细心的工作。每次处理新数据或更换方法后养成交叉验证的习惯能帮你节省大量后续调试和论文返工的时间。
返回列表