简介全球海水表面温度与海冰浓度数据集2020a专用源自 Met Office Hadley Centre 观测数据集包含覆盖全球海域的海表温度和海冰浓度要素是海洋气候研究中常用的基础数据资源适合需要处理 NetCDF 格式但尚未熟悉的新手。压缩包采用7z格式整体约167.98MB共包含3个文件——两个.nc数据文件分别存储全球海表温度与海冰浓度一个.py脚本提供入门级的数据读取与结构检查代码。脚本以简单易懂的方式演示了从打开NetCDF文件到查看变量名、维度、坐标、单位及全局属性的全过程帮助用户快速理清数据构造和变量情况省去自行摸索格式的时间在此基础上使用者可顺利开展海温变化统计、海冰范围监测或基础可视化分析。由于脚本具备通用性也适用于后续其他HadISST系列数据的快速预览与预处理便于扩展分析范围。该资源已有4928人学习下载是研究生、工程师和气候爱好者快速上手实测海洋数据的低门槛选择。1. 全球海水表面温度和海冰浓度数据集2020a 专用它解决什么场景的问题如果你跑过海洋-海冰耦合模型一定对初始场和强迫场的“拼图”过程有印象海水表面温度SST要连续海冰浓度要和 SST 在冰区保持一致时间和网格还要跟模型严格对齐。这套标注“2020a 专用”的全球海水表面温度和海冰浓度数据集本质是一份面向 2020 年模拟时段的再分析后处理产品把 SST 和海冰密集度统一到同一套网格和时间分辨率上省去了从多个数据源拼接和订正的环节。做北极海冰预报、区域海洋环流模拟、以及模型冷启动实验的从业者拿它做侧边界或初始场最合适即便你是刚入手海洋数据处理的新手也能借这份标准 NetCDF 样例把读、查、插、写的完整链路走通。2. 把 2020a 数据集的变量、单位和网格拆开看先看懂再动手拿到数据集先别急着画图或插值我一般会先花十分钟把变量属性、坐标范围和网格类型看清楚。多数“后期翻车”都不是模型问题而是数据的前处理没对齐。下面这几步是固定动作顺序最好不要反。2.1 变量名与单位为什么海冰浓度用 0~1 而不是百分比这类数据集在 NetCDF 里的变量名通常有两套习惯。SST 字段常见的是sst或sea_surface_temperature单位可能是开尔文K也可能是摄氏度degC海冰浓度字段常见的是siconc、sea_ice_conc或sea_ice_area_fraction单位在 CF 标准里一般写成1表示 0~1 之间的面积占比但有些再分析产品直接输出 0~100 的百分数。这两者在插值前后混用是第一个容易踩的坑。我拿到文件后的第一件事是打印属性而不是直接看数值import xarray as xr ds xr.open_dataset(GLB_SST_SIC_2020a.nc, decode_cfTrue, mask_and_scaleTrue) print(ds[sst].attrs) print(ds[sea_ice_conc].attrs if sea_ice_conc in ds else ds)这段代码里decode_cfTrue会让 xarray 自动识别scale_factor和add_offset把 NetCDF 里常见的有符号整型如int16或short还原成真实物理值mask_and_scaleTrue则会把_FillValue自动转成NaN。打印出来的units能直接告诉你两件事SST 需不需要加减 273.15海冰浓度需不需要除以 100。如果属性里写的是units %那么在插值前先sic_da ds[sea_ice_conc] / 100.0这类数据集的变量形态大致可以按下面的表对照字段常见变量名单位典型范围说明海表温度sst/sea_surface_temperatureK / degC-2 ~ 32冰区可能缺失或已被填充海冰浓度siconc/sea_ice_conc/sea_ice_area_fraction1 / %0 ~ 1 或 0~100插值前后务必统一为 0~1经度lon/longitudedegrees_east0~360 或 -180~180决定了后期裁切方式纬度lat/latitudedegrees_north-90~90一般为单调递增如果变量名和你预期不一致不要凭经验硬编码直接打印ds的完整结构看数据集的维度名dims和坐标名。很多时候上游处理脚本换过变量名比如把siconc改成sic你按旧名取数据就会得到KeyError。2.2 网格坐标系规则经纬网格与旋转/位移网格怎么认网格类型直接决定你选哪种重插方法。2020a 这类面向模型驱动的产品最常见的是规则经纬网格经度从 0 到 360 或者从 -180 到 180纬度均匀递增。这种网格最简单用xarray的sel就能按范围切用xESMF做重插也比较稳。但也有些产品源自海洋环流模型输出会带有位移网格staggered grid特征比如变量定义在 B 网格或 C 网格的角点/边上lat和lon是二维数组而不是一维坐标。判断办法很简单打印ds.lat.dims如果是(y, x)而不是(lat,)这就是曲网格这块内容超出了今天的话题不过遇到的情况也不多。更常见的是经度起点不统一——同是 0.25 度网格一个文件从 0 开始另一个从 -180 开始。这不是同一网格重插后会产生半个网格的偏移。在真正动手前我习惯打印坐标范围把网格特征落到纸面上lon_name [c for c in ds.coords if c.lower() in (lon, longitude)][0] lat_name [c for c in ds.coords if c.lower() in (lat, latitude)][0] print(float(ds[lon_name].min()), float(ds[lon_name].max())) print(float(ds[lat_name].min()), float(ds[lat_name].max())) print(ds[lon_name].values[:5])这段输出会让你立刻知道经度是 0~360 还是 -180~180纬度是不是从南极到北极。很多区域模型只关注北半球或中纬度但全局重插前这个信息直接关系到后续裁剪的边界会自动拆成两段的问题。边界断裂通常就出在跨 0 度经线或跨日期变更线的地方。2.3 用 xarray 读取 2020a 数据的最小脚本与输出解读把读取、属性检查和坐标检查放在一起形成每次开文件都会跑的固定模板import xarray as xr f GLB_SST_SIC_2020a.nc ds xr.open_dataset(f, decode_cfTrue, mask_and_scaleTrue) # 看整体维度 print(ds) # 看变量属性和时间坐标 print(ds[sst].attrs) print(ds.time.dtype, ds.time.attrs.get(units)) # 看坐标范围 lon ds.lon if lon in ds.coords else ds.longitude lat ds.lat if lat in ds.coords else ds.latitude print(lon:, float(lon.min()), float(lon.max())) print(lat:, float(lat.min()), float(lat.max()))这里有一个容易被忽略的参数mask_and_scaleTrue。如果你之前发现某个变量全是很奇怪的整数比如 27600那多半是上一次读取时没有应用scale_factor。这个参数在open_dataset里默认就是开启的但如果你用engineh5netcdf或者手动拼接多个文件时偶尔会丢所以我习惯显式写出来。输出解读也简单(time, lat, lon)是标准三维结构SST 和 SIC 应该都是这个形状如果两者维度顺序不一致比如一个是(time, lat, lon)另一个是(lat, lon, time)在后期xarray自动对齐时也能对上但写 NetCDF 时最好统一。时间维度不是必选项——如果文件只给了某一时刻time可能是标量坐标这时后续做时间序列就要先expand_dims我在避坑章节再展开。3. 从全球场到模型驱动文件预处理链路与参数选择数据本身质量再高也必须经过“时间裁剪、空间重插、变量规范化”三道工序才能变成模型认得的驱动文件。这一章的每一步我都给可执行脚本和参数选择依据。3.1 时间子集与去重2020a 时间坐标的常见形态2020a 专用数据集通常覆盖完整年份但模型积分窗口可能只取某个季节。先用slice裁剪时间再处理重复和排序ds ds.sel(timeslice(2020-06-01, 2020-08-31)) # 按时间排序避免文件里记录乱序 ds ds.sortby(time) # 去掉重复时刻数据拼接时偶尔会出现同一时刻两条记录 ds ds.drop_duplicates(time) print(ds.time.values[:3], ds.time.values[-3:])这里slice的参数是字符串xarray 会按时间坐标的基准自动解析。要注意的是如果时间坐标是cftime对象比如日历用了noleap或360_day字符串解析通常是安全的但打印出来的dtype可能是object而不是datetime64[ns]。这种情况不需要立刻转换sel一样能用。为什么先把drop_duplicates放前面因为如果重复的是 7 月 1 日 0 点那一条而它恰好是模型启动时刻重复记录会导致同化或强迫数据在同一time索引上出现两个不同值后期写文件并不报错但模型读进去后时间循环会卡死。3.2 空间裁剪与重插到模型网格xESMF 插值方法怎么选空间处理是整条链路里最需要判断的一步。如果任务是区域模拟我建议先裁剪再重插而不是先全局重插再裁剪。原因是重插算法在全网格上会沿着陆地边界或数据边界额外产生一些怪异值裁剪之后再插能把误差限制在目标范围内。环境准备可以一条命令搞定conda create -n sst2020 python3.11 xarray netcdf4 xesmf esmpy -c conda-forge -y注意xesmf依赖esmpy二者版本最好一起装分开pip install经常出现 ESMF 库版本不匹配的问题。裁剪与重插的典型过程如下import numpy as np import xarray as xr import xesmf as xe ds xr.open_dataset(GLB_SST_SIC_2020a.nc) # 1. 先裁剪到研究区 region ds.sel(latslice(10, 50), lonslice(90, 150)) # 2. 构建目标网格注意 lon 范围必须与模型约定一致 target_grid xr.Dataset( { lat: ([lat], np.arange(10, 51, 0.25)), lon: ([lon], np.arange(90, 150.5, 0.25)), } ) # 3. SST 用双线性海冰浓度用保守插值 regrid_sst xe.Regridder(region[sst], target_grid, bilinear) regrid_sic xe.Regridder(region[sea_ice_conc], target_grid, conservative) sst_rg regrid_sst(region[sst]) sic_rg regrid_sic(region[sea_ice_conc])这里有两个参数需要格外说明。第一bilinear适合连续场海表温度空间梯度本身就平滑用双线性不会引入明显畸变第二海冰浓度是 0~1 的有限值有清晰的“有冰/无冰”边界用conservative能保留总面积但这种方法要求输入网格有bounds如果原始文件没有bounds变量运行会报错。这个时候我一般改成nearest_s2d它把每个源格点值原样搬到距离最近的目标格点不会跨边界平滑代价是空间分辨率略有损失。插值方法适用变量优点风险bilinearSST平滑、连续会穿过陆海边界产生虚假值conservativeSIC、通量守恒性好需要网格 bounds慢nearest_s2dSIC、掩膜不产生中间值分辨率受损、可能出现块状3.3 写出驱动 NetCDF变量顺序、fill_value 与压缩参数重插结果不能直接to_netcdf了事特别是给数值模型用的时候变量类型和缺失值要先规范化。下面是一套推荐的写文件方式ds_out xr.Dataset( { sst: ([time, lat, lon], sst_rg.astype(float32)), sea_ice_conc: ([time, lat, lon], sic_rg.astype(float32)), }, coords{ time: sst_rg.time, lat: sst_rg.lat, lon: sst_rg.lon, }, ) ds_out[sst].attrs[units] degC ds_out[sea_ice_conc].attrs[units] 1 ds_out.to_netcdf( SST_SIC_2020a_region_ready.nc, enginenetcdf4, encoding{ sst: {zlib: True, complevel: 4}, sea_ice_conc: {zlib: True, complevel: 4}, }, )强制float32是因为海洋模型读驱动场时大部分代码用单精度float64会增加文件体积和 I/O 时间。zlib配合complevel4是净收益比较高的组合压缩率通常在 3~5 倍CPU 开销不大。这里没有显式设置_FillValue因为源数据里被mask_and_scale转成NaN的格点会自动在写出时保留为NaN。要确认缺测标记是否一致可以写完后重新打开一次chk xr.open_dataset(SST_SIC_2020a_region_ready.nc) print(chk.sst.encoding.get(_FillValue)) print(chk.sea_ice_conc.encoding.get(_FillValue))这一步不能省模型读到NaN或-9999时表现完全不同。有些模型把NaN当成有效值参与计算第一个时间步就会出现整片区域水温异常。4. 海冰浓度与 SST 的衔接处理掩膜、冰点温度和填充顺序很多人在数据预处理里已经把 SST 和海冰浓度单独处理得“看起来正常”但放到一起就出问题海冰边缘的 SST 高达十几度或者冰下 SST 是明显的陆地填充值。这一章解决的就是“海冰-海温一致”这件事。4.1 冰下 SST 缺失值填充直接填 -1.8 还是临近插值2020a 这类再分析产品里海冰覆盖格点的 SST 有两种可能存在一是真的缺测二是数据源把海洋表层温度算到接近冰点。如果缺测直接填-1.8是最省事的办法但也最粗暴因为不同盐度下的冰点并不一样。更稳妥的做法是先看数据里有没有盐度场有的话按海水冰点经验公式计算逐格点冰点import numpy as np sst_da region[sst].copy() sic_da region[sea_ice_conc].copy() ice_mask sic_da 0.15 # 常见海冰阈值可按模型配置改 sst_missing sst_da.isnull() # 如果数据集附带盐度就用盐度算冰点 if salinity in region: S region[salinity].clip(min0, max40) # 海水冰点经验公式Tf -0.0575*S 1.710523e-3*S^1.5 - 2.154996e-4*S^2 tf -0.0575 * S 1.710523e-3 * S**1.5 - 2.154996e-4 * S**2 sst_da sst_da.where(~(sst_missing ice_mask), tf) else: sst_da sst_da.where(~(sst_missing ice_mask), -1.8)这段代码里最关键的是sst_missing ice_mask它限定了只填充“海冰覆盖且 SST 缺测”的格点而不是把所有NaN都填成冰点。之前有同事图省事直接fillna(-1.8)结果把陆地掩膜格点也填成了海温模型海岸线附近冷得一塌糊涂。需要警惕的是ice_mask的阈值不是固定的。气候模式常用 0.15海冰预报模型有时用 0.1天气尺度强迫场用 0.5 的也有。这个阈值不是数据集的属性而是模型物理方案的设置宁可多检查一遍模型手册也别照着别人配置抄。4.2 海冰密集度阈值与 SST 掩膜的先后顺序另一个常见争议是先做掩膜还是先做插值。我推荐的处理顺序是源数据掩膜 → 插值 → 以海冰阈值重设 SST → 再用模型陆海掩膜裁剪。如果先按模型掩膜裁剪再插值源数据的边界会和模型掩膜边界错位海冰浓度很容易沿着海岸线渗进陆地格点。实际操作时我会把陆海掩膜和格点面积权重一起处理# 假设模型掩膜 1海洋, 0陆地 landsea_model xr.open_dataset(model_landsea.nc)[mask] # 1. 保护源数据掩膜把陆地先设置为 NaN source_ocean np.isfinite(sic_rg) sic_rg_masked sic_rg.where(source_ocean) # 2. 海冰覆盖格点SST 重设为冰点 sst_final sst_rg.where(~(sic_rg_masked 0.15), sst_filled) # 3. 用模型陆海掩膜做最后裁剪 sst_final sst_final.where(landsea_model 0.5) sic_final sic_rg_masked.where(landsea_model 0.5)为什么要插值之后再设 SST因为重插过程中海冰浓度会被平滑原本 0.13 的格点可能变 0.16也可能原本 0.18 的格点变 0.14。如果在插值之前就按源数据阈值把 SST 填成冰点插值之后 SST 会再次被邻近暖水污染前面等于白做。4.3 生成一张“海冰-海温一致”驱动场的标准流程示例把以上思路串成一个完整示例适合直接改成你的任务脚本import numpy as np import xarray as xr import xesmf as xe # 读源数据 ds xr.open_dataset(GLB_SST_SIC_2020a.nc).sel(time2020-07-15) # 重插到目标网格 target xr.Dataset({ lat: ([lat], np.arange(60, 90, 0.25)), lon: ([lon], np.arange(-180, 180, 0.25)), }) regridder_sst xe.Regridder(ds[sst], target, bilinear) regridder_sic xe.Regridder(ds[sea_ice_conc], target, nearest_s2d) sst regridder_sst(ds[sst]).astype(float32) sic regridder_sic(ds[sea_ice_conc]).astype(float32) # 海冰阈值判定 ice_threshold 0.15 ice_mask sic ice_threshold # 冰区 SST 统一设为 -1.8非冰区保留原值 sst sst.where(~ice_mask, -1.8) # 模型海冰浓度低于阈值的按 0 处理 sic sic.where(sic ice_threshold, 0.0) # 写文件 out xr.Dataset( {sst: sst, sea_ice_conc: sic}, coords{lat: target.lat, lon: target.lon}, ) out.to_netcdf(arctic_sst_sic_2020a.nc, enginenetcdf4, encoding{sst: {zlib: True, complevel: 4}, sea_ice_conc: {zlib: True, complevel: 4}})这段流程里最容易被新手忽略的是最后一步sic.where(sic ice_threshold, 0.0)。如果不做海冰浓度在 0.05~0.14 的格点会保留一个“亚阈值冰量”模型的海冰热力过程会把这些薄冰格点当成有效冰面反射太阳短波辐射热量收支一下子就偏了。这个值到底设多少取决于你对模型物理方案的理解但处理流程本身必须要有。5. 常见问题避坑处理 2020a 海冰/SST 数据时遇到的 5 条真实踩坑记录这一章写给我自己在内的后来人。每一条都是真实会发生的现象按“现象 → 原因 → 解决”展开希望能帮你少做几次无用功。5.1 现象一海冰浓度插值后出现负值和超过 100%现象从 0~1 的海冰浓度场做双线性插值后数组最小值到了 -0.3最大值到了 1.2。原因海冰浓度是一个有明确上下界的物理量但双线性插值本质是加权平均不会主动约束边界。当源网格里有一侧是无冰区、邻侧是满冰区插值窗口横跨两个格点时受缺失值和网格边缘影响可以出现超界。解决不依赖插值器自动处理在写出前显式约束sic_rg sic_rg.clip(min0, max1)如果nearest_s2d仍然出现超界通常是因为源数据本身有异常值。这时要先print(sic_rg.min(), sic_rg.max())确认源头而不是直接clip否则问题会被掩盖到后面。5.2 现象二海冰覆盖格点上出现 25℃ 的 SST现象重插后的 SST 场里海冰浓度达到 0.9 的格点海温显示 25℃。原因源 SST 在海冰覆盖区本来就是缺测或无效值插值时这些格点附近的暖水被平均进来或者源数据里的海温来自“无冰期气候态”导致数值偏高。再往深一步可能是源数据集的海温与海冰掩膜不是同步生成的。解决不要在插值后手工挑异常点而是用海冰浓度做强制修正sst_final sst_rg.where(~(sic_rg 0.15), -1.8)关键点是海冰浓度必须先于 SST 处理顺序不能反。如果再叠上一层“SST 低于 -2.5℃ 视为无效”的过滤条件也能兜住一部分问题但治标不治本。5.3 现象三写出的文件经度范围 0~360模型读取后多出一条裂痕现象区域模型运行后在太平洋中部出现一条南北向的异常边界流场沿这条线断开。原因源数据经度是 0~360而模型网格内部约定是 -180~180。重插到目标网格时我把目标lon设成了np.arange(-180, 180, 0.25)但源数据经度大于 180 的部分没有先做转换导致数据在 180 度处出现跳跃。解决在裁剪与插值前统一经度范围lon ds[lon] if float(lon.max()) 180: # 0~360 - -180~180 lon_adj ((lon 180) % 360) - 180 ds ds.assign_coords(lonlon_adj).sortby(lon)注意这里用了取模和减法而不是lon - 360因为源数据经度可能到 360 附近直接减会弄成负数范围。转换后一定要sortby否则经度序列乱序插值器会按原始顺序处理结果是另一套错位。5.4 现象四时间坐标读取出来全是同一日期或整体偏移了 8 小时现象ds.time.values打印出来所有值都一样或者 2020 年 1 月 1 日 00:00 的场次被读成了 08:00。原因前一种情况多半是时间坐标在文件里被定义为标量没有作为维度展开或者 NetCDF 变量time的维度是time但长度为 1没有沿着时间轴展开后一种情况是时间单位里包含了时区基准比如hours since 2020-01-01 00:00:00用的是 UTC而部分国产数据产品按北京时间UTC8输出decode_cfTrue会按字符串里的基准来不会自动加时区。解决先确认维度再决定要不要展开if ds.time.ndim 0: ds ds.expand_dims(time) print(ds.time.attrs.get(units))如果是 UTC 基准但业务上需要北京时手动加偏移并写清attrsds ds.assign_coords(timeds.time np.timedelta64(8, h)) ds.time.attrs[units] hours since 2020-01-01 00:00:0008这里最容易翻车的是你对“偏移 8 小时”的理解。先确认源数据文档里写的是 UTC 还是本地时再决定加不加。盲目加偏移等于让模型初始时刻整体平移对强强迫问题的结果影响不大对海冰日变化敏感的问题会直接影响日循环相位。5.5 现象五岸线附近重插出“碎冰带”海冰蔓延到内陆网格现象海冰浓度场在格陵兰岛东岸和加拿大北岸出现大量离散的碎冰格点甚至越过海岸线出现在陆地上。原因插值器沿陆海边界处理时把陆地缺失值当成 0 参与加权平均或者conservative插值法把陆侧格点的面积权重也算进了海冰浓度。源数据的陆海掩膜和模型陆海掩膜如果分辨率不一致问题更严重。解决海冰浓度建议优先用nearest_s2d这类插值不会产生跨边界的中间值。如果必须用conservative那就先把源数据的陆地格点设成NaNsic_source ds[sea_ice_conc].where(ds[sea_ice_conc] 0) regridder_sic xe.Regridder(sic_source, target, conservative) sic_rg regridder_sic(sic_source)注意where(sea_ice_conc 0)会把陆地上的负值缺测标记滤除。插值结束后还要跟模型掩膜再乘一次才能保证最终文件里陆地格点的海冰浓度严格为 0。6. 进阶用法用独立观测验证字段再给初始场做一次平滑数据处理完只是第一步真正让我放心的验证方式是把结果和独立观测放在一起算偏差。对于 2020a 专用数据集最常见的独立参照是 NOAA OISST 月平均海温和浮标/Argo 剖面。这里给一套足够轻量的验证脚本你不需要完整的观测网只要有一份 OISST 就能跑通。6.1 快速验证与 OISST 月平均做偏差统计oisst xr.open_dataset(OISST_2020.nc)[sst] # 先把两个场对齐到同一套坐标这里省略重插细节 diff ds_out[sst].sel(timeslice(2020-07-01, 2020-09-30)) - oisst[sst] print(diff.mean(dim[lat, lon]).values) print(diff.sel(latslice(-60, 60)).std(dim[time, lat, lon]).values)第一行输出是平均偏差正负只说明系统高低第二行是标准差如果超过 1.5℃ 就要回头检查是不是经度范围或者海冰掩膜处理出了问题。对北极区域标准差的警戒线通常会更高一点因为海冰边缘本身就是高变率区但均值偏差大于 1℃ 仍然值得警惕。6.2 初始场平滑冷启动压力波的抑制最后再追加一个技巧。直接把再分析场塞进海洋模型做冷启动初始场在陆架陡坡区域容易出现“压力波”头一两天模拟结果有明显的高频振荡。常见做法是给初始 SST 做一次浅层平滑让海洋表面不再呈现格点尺度的锯齿。from scipy.ndimage import gaussian_filter import numpy as np sst_np ds_out[sst].values sst_smooth_np gaussian_filter(sst_np, sigma1.5, modenearest) sst_smooth xr.DataArray(sst_smooth_np, coordsds_out[sst].coords, attrsds_out[sst].attrs) # 平滑会改变全场平均海温简单加回一个常量做守恒修正 sst_smooth sst_smooth (ds_out[sst].mean() - sst_smooth.mean())sigma1.5意味着在 0.25 度网格上影响半径约 0.375 度是一个很轻的平滑如果网格更粗比如 1 度我会把sigma降到 1。modenearest是为了避免边界外推产生异常值。海冰浓度不要做这种平滑它的物理边界比温度更锐利平滑之后阈值判定会失真。这套方案我前后用了快两年最大的体会是海冰浓度永远比 SST 先处理填缺失值永远比插值先处理。每次“看起来都正常”的输出最后翻车都翻在那些没打印过的属性和没有核对过的坐标范围上。希望帮到你。本文还有配套的精品资源点击获取 海温海冰数据预处理实战：海洋-海冰模型驱动场构建指南
简介全球海水表面温度与海冰浓度数据集2020a专用源自 Met Office Hadley Centre 观测数据集包含覆盖全球海域的海表温度和海冰浓度要素是海洋气候研究中常用的基础数据资源适合需要处理 NetCDF 格式但尚未熟悉的新手。压缩包采用7z格式整体约167.98MB共包含3个文件——两个.nc数据文件分别存储全球海表温度与海冰浓度一个.py脚本提供入门级的数据读取与结构检查代码。脚本以简单易懂的方式演示了从打开NetCDF文件到查看变量名、维度、坐标、单位及全局属性的全过程帮助用户快速理清数据构造和变量情况省去自行摸索格式的时间在此基础上使用者可顺利开展海温变化统计、海冰范围监测或基础可视化分析。由于脚本具备通用性也适用于后续其他HadISST系列数据的快速预览与预处理便于扩展分析范围。该资源已有4928人学习下载是研究生、工程师和气候爱好者快速上手实测海洋数据的低门槛选择。1. 全球海水表面温度和海冰浓度数据集2020a 专用它解决什么场景的问题如果你跑过海洋-海冰耦合模型一定对初始场和强迫场的“拼图”过程有印象海水表面温度SST要连续海冰浓度要和 SST 在冰区保持一致时间和网格还要跟模型严格对齐。这套标注“2020a 专用”的全球海水表面温度和海冰浓度数据集本质是一份面向 2020 年模拟时段的再分析后处理产品把 SST 和海冰密集度统一到同一套网格和时间分辨率上省去了从多个数据源拼接和订正的环节。做北极海冰预报、区域海洋环流模拟、以及模型冷启动实验的从业者拿它做侧边界或初始场最合适即便你是刚入手海洋数据处理的新手也能借这份标准 NetCDF 样例把读、查、插、写的完整链路走通。2. 把 2020a 数据集的变量、单位和网格拆开看先看懂再动手拿到数据集先别急着画图或插值我一般会先花十分钟把变量属性、坐标范围和网格类型看清楚。多数“后期翻车”都不是模型问题而是数据的前处理没对齐。下面这几步是固定动作顺序最好不要反。2.1 变量名与单位为什么海冰浓度用 0~1 而不是百分比这类数据集在 NetCDF 里的变量名通常有两套习惯。SST 字段常见的是sst或sea_surface_temperature单位可能是开尔文K也可能是摄氏度degC海冰浓度字段常见的是siconc、sea_ice_conc或sea_ice_area_fraction单位在 CF 标准里一般写成1表示 0~1 之间的面积占比但有些再分析产品直接输出 0~100 的百分数。这两者在插值前后混用是第一个容易踩的坑。我拿到文件后的第一件事是打印属性而不是直接看数值import xarray as xr ds xr.open_dataset(GLB_SST_SIC_2020a.nc, decode_cfTrue, mask_and_scaleTrue) print(ds[sst].attrs) print(ds[sea_ice_conc].attrs if sea_ice_conc in ds else ds)这段代码里decode_cfTrue会让 xarray 自动识别scale_factor和add_offset把 NetCDF 里常见的有符号整型如int16或short还原成真实物理值mask_and_scaleTrue则会把_FillValue自动转成NaN。打印出来的units能直接告诉你两件事SST 需不需要加减 273.15海冰浓度需不需要除以 100。如果属性里写的是units %那么在插值前先sic_da ds[sea_ice_conc] / 100.0这类数据集的变量形态大致可以按下面的表对照字段常见变量名单位典型范围说明海表温度sst/sea_surface_temperatureK / degC-2 ~ 32冰区可能缺失或已被填充海冰浓度siconc/sea_ice_conc/sea_ice_area_fraction1 / %0 ~ 1 或 0~100插值前后务必统一为 0~1经度lon/longitudedegrees_east0~360 或 -180~180决定了后期裁切方式纬度lat/latitudedegrees_north-90~90一般为单调递增如果变量名和你预期不一致不要凭经验硬编码直接打印ds的完整结构看数据集的维度名dims和坐标名。很多时候上游处理脚本换过变量名比如把siconc改成sic你按旧名取数据就会得到KeyError。2.2 网格坐标系规则经纬网格与旋转/位移网格怎么认网格类型直接决定你选哪种重插方法。2020a 这类面向模型驱动的产品最常见的是规则经纬网格经度从 0 到 360 或者从 -180 到 180纬度均匀递增。这种网格最简单用xarray的sel就能按范围切用xESMF做重插也比较稳。但也有些产品源自海洋环流模型输出会带有位移网格staggered grid特征比如变量定义在 B 网格或 C 网格的角点/边上lat和lon是二维数组而不是一维坐标。判断办法很简单打印ds.lat.dims如果是(y, x)而不是(lat,)这就是曲网格这块内容超出了今天的话题不过遇到的情况也不多。更常见的是经度起点不统一——同是 0.25 度网格一个文件从 0 开始另一个从 -180 开始。这不是同一网格重插后会产生半个网格的偏移。在真正动手前我习惯打印坐标范围把网格特征落到纸面上lon_name [c for c in ds.coords if c.lower() in (lon, longitude)][0] lat_name [c for c in ds.coords if c.lower() in (lat, latitude)][0] print(float(ds[lon_name].min()), float(ds[lon_name].max())) print(float(ds[lat_name].min()), float(ds[lat_name].max())) print(ds[lon_name].values[:5])这段输出会让你立刻知道经度是 0~360 还是 -180~180纬度是不是从南极到北极。很多区域模型只关注北半球或中纬度但全局重插前这个信息直接关系到后续裁剪的边界会自动拆成两段的问题。边界断裂通常就出在跨 0 度经线或跨日期变更线的地方。2.3 用 xarray 读取 2020a 数据的最小脚本与输出解读把读取、属性检查和坐标检查放在一起形成每次开文件都会跑的固定模板import xarray as xr f GLB_SST_SIC_2020a.nc ds xr.open_dataset(f, decode_cfTrue, mask_and_scaleTrue) # 看整体维度 print(ds) # 看变量属性和时间坐标 print(ds[sst].attrs) print(ds.time.dtype, ds.time.attrs.get(units)) # 看坐标范围 lon ds.lon if lon in ds.coords else ds.longitude lat ds.lat if lat in ds.coords else ds.latitude print(lon:, float(lon.min()), float(lon.max())) print(lat:, float(lat.min()), float(lat.max()))这里有一个容易被忽略的参数mask_and_scaleTrue。如果你之前发现某个变量全是很奇怪的整数比如 27600那多半是上一次读取时没有应用scale_factor。这个参数在open_dataset里默认就是开启的但如果你用engineh5netcdf或者手动拼接多个文件时偶尔会丢所以我习惯显式写出来。输出解读也简单(time, lat, lon)是标准三维结构SST 和 SIC 应该都是这个形状如果两者维度顺序不一致比如一个是(time, lat, lon)另一个是(lat, lon, time)在后期xarray自动对齐时也能对上但写 NetCDF 时最好统一。时间维度不是必选项——如果文件只给了某一时刻time可能是标量坐标这时后续做时间序列就要先expand_dims我在避坑章节再展开。3. 从全球场到模型驱动文件预处理链路与参数选择数据本身质量再高也必须经过“时间裁剪、空间重插、变量规范化”三道工序才能变成模型认得的驱动文件。这一章的每一步我都给可执行脚本和参数选择依据。3.1 时间子集与去重2020a 时间坐标的常见形态2020a 专用数据集通常覆盖完整年份但模型积分窗口可能只取某个季节。先用slice裁剪时间再处理重复和排序ds ds.sel(timeslice(2020-06-01, 2020-08-31)) # 按时间排序避免文件里记录乱序 ds ds.sortby(time) # 去掉重复时刻数据拼接时偶尔会出现同一时刻两条记录 ds ds.drop_duplicates(time) print(ds.time.values[:3], ds.time.values[-3:])这里slice的参数是字符串xarray 会按时间坐标的基准自动解析。要注意的是如果时间坐标是cftime对象比如日历用了noleap或360_day字符串解析通常是安全的但打印出来的dtype可能是object而不是datetime64[ns]。这种情况不需要立刻转换sel一样能用。为什么先把drop_duplicates放前面因为如果重复的是 7 月 1 日 0 点那一条而它恰好是模型启动时刻重复记录会导致同化或强迫数据在同一time索引上出现两个不同值后期写文件并不报错但模型读进去后时间循环会卡死。3.2 空间裁剪与重插到模型网格xESMF 插值方法怎么选空间处理是整条链路里最需要判断的一步。如果任务是区域模拟我建议先裁剪再重插而不是先全局重插再裁剪。原因是重插算法在全网格上会沿着陆地边界或数据边界额外产生一些怪异值裁剪之后再插能把误差限制在目标范围内。环境准备可以一条命令搞定conda create -n sst2020 python3.11 xarray netcdf4 xesmf esmpy -c conda-forge -y注意xesmf依赖esmpy二者版本最好一起装分开pip install经常出现 ESMF 库版本不匹配的问题。裁剪与重插的典型过程如下import numpy as np import xarray as xr import xesmf as xe ds xr.open_dataset(GLB_SST_SIC_2020a.nc) # 1. 先裁剪到研究区 region ds.sel(latslice(10, 50), lonslice(90, 150)) # 2. 构建目标网格注意 lon 范围必须与模型约定一致 target_grid xr.Dataset( { lat: ([lat], np.arange(10, 51, 0.25)), lon: ([lon], np.arange(90, 150.5, 0.25)), } ) # 3. SST 用双线性海冰浓度用保守插值 regrid_sst xe.Regridder(region[sst], target_grid, bilinear) regrid_sic xe.Regridder(region[sea_ice_conc], target_grid, conservative) sst_rg regrid_sst(region[sst]) sic_rg regrid_sic(region[sea_ice_conc])这里有两个参数需要格外说明。第一bilinear适合连续场海表温度空间梯度本身就平滑用双线性不会引入明显畸变第二海冰浓度是 0~1 的有限值有清晰的“有冰/无冰”边界用conservative能保留总面积但这种方法要求输入网格有bounds如果原始文件没有bounds变量运行会报错。这个时候我一般改成nearest_s2d它把每个源格点值原样搬到距离最近的目标格点不会跨边界平滑代价是空间分辨率略有损失。插值方法适用变量优点风险bilinearSST平滑、连续会穿过陆海边界产生虚假值conservativeSIC、通量守恒性好需要网格 bounds慢nearest_s2dSIC、掩膜不产生中间值分辨率受损、可能出现块状3.3 写出驱动 NetCDF变量顺序、fill_value 与压缩参数重插结果不能直接to_netcdf了事特别是给数值模型用的时候变量类型和缺失值要先规范化。下面是一套推荐的写文件方式ds_out xr.Dataset( { sst: ([time, lat, lon], sst_rg.astype(float32)), sea_ice_conc: ([time, lat, lon], sic_rg.astype(float32)), }, coords{ time: sst_rg.time, lat: sst_rg.lat, lon: sst_rg.lon, }, ) ds_out[sst].attrs[units] degC ds_out[sea_ice_conc].attrs[units] 1 ds_out.to_netcdf( SST_SIC_2020a_region_ready.nc, enginenetcdf4, encoding{ sst: {zlib: True, complevel: 4}, sea_ice_conc: {zlib: True, complevel: 4}, }, )强制float32是因为海洋模型读驱动场时大部分代码用单精度float64会增加文件体积和 I/O 时间。zlib配合complevel4是净收益比较高的组合压缩率通常在 3~5 倍CPU 开销不大。这里没有显式设置_FillValue因为源数据里被mask_and_scale转成NaN的格点会自动在写出时保留为NaN。要确认缺测标记是否一致可以写完后重新打开一次chk xr.open_dataset(SST_SIC_2020a_region_ready.nc) print(chk.sst.encoding.get(_FillValue)) print(chk.sea_ice_conc.encoding.get(_FillValue))这一步不能省模型读到NaN或-9999时表现完全不同。有些模型把NaN当成有效值参与计算第一个时间步就会出现整片区域水温异常。4. 海冰浓度与 SST 的衔接处理掩膜、冰点温度和填充顺序很多人在数据预处理里已经把 SST 和海冰浓度单独处理得“看起来正常”但放到一起就出问题海冰边缘的 SST 高达十几度或者冰下 SST 是明显的陆地填充值。这一章解决的就是“海冰-海温一致”这件事。4.1 冰下 SST 缺失值填充直接填 -1.8 还是临近插值2020a 这类再分析产品里海冰覆盖格点的 SST 有两种可能存在一是真的缺测二是数据源把海洋表层温度算到接近冰点。如果缺测直接填-1.8是最省事的办法但也最粗暴因为不同盐度下的冰点并不一样。更稳妥的做法是先看数据里有没有盐度场有的话按海水冰点经验公式计算逐格点冰点import numpy as np sst_da region[sst].copy() sic_da region[sea_ice_conc].copy() ice_mask sic_da 0.15 # 常见海冰阈值可按模型配置改 sst_missing sst_da.isnull() # 如果数据集附带盐度就用盐度算冰点 if salinity in region: S region[salinity].clip(min0, max40) # 海水冰点经验公式Tf -0.0575*S 1.710523e-3*S^1.5 - 2.154996e-4*S^2 tf -0.0575 * S 1.710523e-3 * S**1.5 - 2.154996e-4 * S**2 sst_da sst_da.where(~(sst_missing ice_mask), tf) else: sst_da sst_da.where(~(sst_missing ice_mask), -1.8)这段代码里最关键的是sst_missing ice_mask它限定了只填充“海冰覆盖且 SST 缺测”的格点而不是把所有NaN都填成冰点。之前有同事图省事直接fillna(-1.8)结果把陆地掩膜格点也填成了海温模型海岸线附近冷得一塌糊涂。需要警惕的是ice_mask的阈值不是固定的。气候模式常用 0.15海冰预报模型有时用 0.1天气尺度强迫场用 0.5 的也有。这个阈值不是数据集的属性而是模型物理方案的设置宁可多检查一遍模型手册也别照着别人配置抄。4.2 海冰密集度阈值与 SST 掩膜的先后顺序另一个常见争议是先做掩膜还是先做插值。我推荐的处理顺序是源数据掩膜 → 插值 → 以海冰阈值重设 SST → 再用模型陆海掩膜裁剪。如果先按模型掩膜裁剪再插值源数据的边界会和模型掩膜边界错位海冰浓度很容易沿着海岸线渗进陆地格点。实际操作时我会把陆海掩膜和格点面积权重一起处理# 假设模型掩膜 1海洋, 0陆地 landsea_model xr.open_dataset(model_landsea.nc)[mask] # 1. 保护源数据掩膜把陆地先设置为 NaN source_ocean np.isfinite(sic_rg) sic_rg_masked sic_rg.where(source_ocean) # 2. 海冰覆盖格点SST 重设为冰点 sst_final sst_rg.where(~(sic_rg_masked 0.15), sst_filled) # 3. 用模型陆海掩膜做最后裁剪 sst_final sst_final.where(landsea_model 0.5) sic_final sic_rg_masked.where(landsea_model 0.5)为什么要插值之后再设 SST因为重插过程中海冰浓度会被平滑原本 0.13 的格点可能变 0.16也可能原本 0.18 的格点变 0.14。如果在插值之前就按源数据阈值把 SST 填成冰点插值之后 SST 会再次被邻近暖水污染前面等于白做。4.3 生成一张“海冰-海温一致”驱动场的标准流程示例把以上思路串成一个完整示例适合直接改成你的任务脚本import numpy as np import xarray as xr import xesmf as xe # 读源数据 ds xr.open_dataset(GLB_SST_SIC_2020a.nc).sel(time2020-07-15) # 重插到目标网格 target xr.Dataset({ lat: ([lat], np.arange(60, 90, 0.25)), lon: ([lon], np.arange(-180, 180, 0.25)), }) regridder_sst xe.Regridder(ds[sst], target, bilinear) regridder_sic xe.Regridder(ds[sea_ice_conc], target, nearest_s2d) sst regridder_sst(ds[sst]).astype(float32) sic regridder_sic(ds[sea_ice_conc]).astype(float32) # 海冰阈值判定 ice_threshold 0.15 ice_mask sic ice_threshold # 冰区 SST 统一设为 -1.8非冰区保留原值 sst sst.where(~ice_mask, -1.8) # 模型海冰浓度低于阈值的按 0 处理 sic sic.where(sic ice_threshold, 0.0) # 写文件 out xr.Dataset( {sst: sst, sea_ice_conc: sic}, coords{lat: target.lat, lon: target.lon}, ) out.to_netcdf(arctic_sst_sic_2020a.nc, enginenetcdf4, encoding{sst: {zlib: True, complevel: 4}, sea_ice_conc: {zlib: True, complevel: 4}})这段流程里最容易被新手忽略的是最后一步sic.where(sic ice_threshold, 0.0)。如果不做海冰浓度在 0.05~0.14 的格点会保留一个“亚阈值冰量”模型的海冰热力过程会把这些薄冰格点当成有效冰面反射太阳短波辐射热量收支一下子就偏了。这个值到底设多少取决于你对模型物理方案的理解但处理流程本身必须要有。5. 常见问题避坑处理 2020a 海冰/SST 数据时遇到的 5 条真实踩坑记录这一章写给我自己在内的后来人。每一条都是真实会发生的现象按“现象 → 原因 → 解决”展开希望能帮你少做几次无用功。5.1 现象一海冰浓度插值后出现负值和超过 100%现象从 0~1 的海冰浓度场做双线性插值后数组最小值到了 -0.3最大值到了 1.2。原因海冰浓度是一个有明确上下界的物理量但双线性插值本质是加权平均不会主动约束边界。当源网格里有一侧是无冰区、邻侧是满冰区插值窗口横跨两个格点时受缺失值和网格边缘影响可以出现超界。解决不依赖插值器自动处理在写出前显式约束sic_rg sic_rg.clip(min0, max1)如果nearest_s2d仍然出现超界通常是因为源数据本身有异常值。这时要先print(sic_rg.min(), sic_rg.max())确认源头而不是直接clip否则问题会被掩盖到后面。5.2 现象二海冰覆盖格点上出现 25℃ 的 SST现象重插后的 SST 场里海冰浓度达到 0.9 的格点海温显示 25℃。原因源 SST 在海冰覆盖区本来就是缺测或无效值插值时这些格点附近的暖水被平均进来或者源数据里的海温来自“无冰期气候态”导致数值偏高。再往深一步可能是源数据集的海温与海冰掩膜不是同步生成的。解决不要在插值后手工挑异常点而是用海冰浓度做强制修正sst_final sst_rg.where(~(sic_rg 0.15), -1.8)关键点是海冰浓度必须先于 SST 处理顺序不能反。如果再叠上一层“SST 低于 -2.5℃ 视为无效”的过滤条件也能兜住一部分问题但治标不治本。5.3 现象三写出的文件经度范围 0~360模型读取后多出一条裂痕现象区域模型运行后在太平洋中部出现一条南北向的异常边界流场沿这条线断开。原因源数据经度是 0~360而模型网格内部约定是 -180~180。重插到目标网格时我把目标lon设成了np.arange(-180, 180, 0.25)但源数据经度大于 180 的部分没有先做转换导致数据在 180 度处出现跳跃。解决在裁剪与插值前统一经度范围lon ds[lon] if float(lon.max()) 180: # 0~360 - -180~180 lon_adj ((lon 180) % 360) - 180 ds ds.assign_coords(lonlon_adj).sortby(lon)注意这里用了取模和减法而不是lon - 360因为源数据经度可能到 360 附近直接减会弄成负数范围。转换后一定要sortby否则经度序列乱序插值器会按原始顺序处理结果是另一套错位。5.4 现象四时间坐标读取出来全是同一日期或整体偏移了 8 小时现象ds.time.values打印出来所有值都一样或者 2020 年 1 月 1 日 00:00 的场次被读成了 08:00。原因前一种情况多半是时间坐标在文件里被定义为标量没有作为维度展开或者 NetCDF 变量time的维度是time但长度为 1没有沿着时间轴展开后一种情况是时间单位里包含了时区基准比如hours since 2020-01-01 00:00:00用的是 UTC而部分国产数据产品按北京时间UTC8输出decode_cfTrue会按字符串里的基准来不会自动加时区。解决先确认维度再决定要不要展开if ds.time.ndim 0: ds ds.expand_dims(time) print(ds.time.attrs.get(units))如果是 UTC 基准但业务上需要北京时手动加偏移并写清attrsds ds.assign_coords(timeds.time np.timedelta64(8, h)) ds.time.attrs[units] hours since 2020-01-01 00:00:0008这里最容易翻车的是你对“偏移 8 小时”的理解。先确认源数据文档里写的是 UTC 还是本地时再决定加不加。盲目加偏移等于让模型初始时刻整体平移对强强迫问题的结果影响不大对海冰日变化敏感的问题会直接影响日循环相位。5.5 现象五岸线附近重插出“碎冰带”海冰蔓延到内陆网格现象海冰浓度场在格陵兰岛东岸和加拿大北岸出现大量离散的碎冰格点甚至越过海岸线出现在陆地上。原因插值器沿陆海边界处理时把陆地缺失值当成 0 参与加权平均或者conservative插值法把陆侧格点的面积权重也算进了海冰浓度。源数据的陆海掩膜和模型陆海掩膜如果分辨率不一致问题更严重。解决海冰浓度建议优先用nearest_s2d这类插值不会产生跨边界的中间值。如果必须用conservative那就先把源数据的陆地格点设成NaNsic_source ds[sea_ice_conc].where(ds[sea_ice_conc] 0) regridder_sic xe.Regridder(sic_source, target, conservative) sic_rg regridder_sic(sic_source)注意where(sea_ice_conc 0)会把陆地上的负值缺测标记滤除。插值结束后还要跟模型掩膜再乘一次才能保证最终文件里陆地格点的海冰浓度严格为 0。6. 进阶用法用独立观测验证字段再给初始场做一次平滑数据处理完只是第一步真正让我放心的验证方式是把结果和独立观测放在一起算偏差。对于 2020a 专用数据集最常见的独立参照是 NOAA OISST 月平均海温和浮标/Argo 剖面。这里给一套足够轻量的验证脚本你不需要完整的观测网只要有一份 OISST 就能跑通。6.1 快速验证与 OISST 月平均做偏差统计oisst xr.open_dataset(OISST_2020.nc)[sst] # 先把两个场对齐到同一套坐标这里省略重插细节 diff ds_out[sst].sel(timeslice(2020-07-01, 2020-09-30)) - oisst[sst] print(diff.mean(dim[lat, lon]).values) print(diff.sel(latslice(-60, 60)).std(dim[time, lat, lon]).values)第一行输出是平均偏差正负只说明系统高低第二行是标准差如果超过 1.5℃ 就要回头检查是不是经度范围或者海冰掩膜处理出了问题。对北极区域标准差的警戒线通常会更高一点因为海冰边缘本身就是高变率区但均值偏差大于 1℃ 仍然值得警惕。6.2 初始场平滑冷启动压力波的抑制最后再追加一个技巧。直接把再分析场塞进海洋模型做冷启动初始场在陆架陡坡区域容易出现“压力波”头一两天模拟结果有明显的高频振荡。常见做法是给初始 SST 做一次浅层平滑让海洋表面不再呈现格点尺度的锯齿。from scipy.ndimage import gaussian_filter import numpy as np sst_np ds_out[sst].values sst_smooth_np gaussian_filter(sst_np, sigma1.5, modenearest) sst_smooth xr.DataArray(sst_smooth_np, coordsds_out[sst].coords, attrsds_out[sst].attrs) # 平滑会改变全场平均海温简单加回一个常量做守恒修正 sst_smooth sst_smooth (ds_out[sst].mean() - sst_smooth.mean())sigma1.5意味着在 0.25 度网格上影响半径约 0.375 度是一个很轻的平滑如果网格更粗比如 1 度我会把sigma降到 1。modenearest是为了避免边界外推产生异常值。海冰浓度不要做这种平滑它的物理边界比温度更锐利平滑之后阈值判定会失真。这套方案我前后用了快两年最大的体会是海冰浓度永远比 SST 先处理填缺失值永远比插值先处理。每次“看起来都正常”的输出最后翻车都翻在那些没打印过的属性和没有核对过的坐标范围上。希望帮到你。本文还有配套的精品资源点击获取