简介:全球海水表面温度与海冰浓度数据集(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_cf=True, mask_and_scale=True) print(ds["sst"].attrs) print(ds["sea_ice_conc"].attrs if "sea_ice_conc" in ds else ds)这段代码里decode_cf=True会让 xarray 自动识别scale_factor和add_offset,把 NetCDF 里常见的有符号整型(如int16或short)还原成真实物理值;mask_and_scale=True则会把_FillValue自动转成NaN。打印出来的units能直接告诉你两件事:SST 需不需要加减 273.15,海冰浓度需不需要除以 100。如果属性里写的是units = %,那么在插值前先:
sic_da = ds["sea_ice_conc"] / 100.0这类数据集的变量形态大致可以按下面的表对照:
| 字段 | 常见变量名 | 单位 | 典型范围 | 说明 |
|---|---|---|---|---|
| 海表温度 | sst/sea_surface_temperature | K / degC | -2 ~ 32 | 冰区可能缺失或已被填充 |
| 海冰浓度 | siconc/sea_ice_conc/sea_ice_area_fraction | 1 / % | 0 ~ 1 或 0~100 | 插值前后务必统一为 0~1 |
| 经度 | lon/longitude | degrees_east | 0~360 或 -180~180 | 决定了后期裁切方式 |
| 纬度 | lat/latitude | degrees_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_cf=True, mask_and_scale=True) # 看整体维度 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_scale=True。如果你之前发现某个变量全是很奇怪的整数,比如 27600,那多半是上一次读取时没有应用scale_factor。这个参数在open_dataset里默认就是开启的,但如果你用engine="h5netcdf"或者手动拼接多个文件时,偶尔会丢,所以我习惯显式写出来。
输出解读也简单:(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(time=slice("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 python=3.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(lat=slice(10, 50), lon=slice(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,它把每个源格点值原样搬到距离最近的目标格点,不会跨边界平滑,代价是空间分辨率略有损失。
| 插值方法 | 适用变量 | 优点 | 风险 |
|---|---|---|---|
| bilinear | SST | 平滑、连续 | 会穿过陆海边界产生虚假值 |
| conservative | SIC、通量 | 守恒性好 | 需要网格 bounds,慢 |
| nearest_s2d | SIC、掩膜 | 不产生中间值 | 分辨率受损、可能出现块状 |
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", engine="netcdf4", encoding={ "sst": {"zlib": True, "complevel": 4}, "sea_ice_conc": {"zlib": True, "complevel": 4}, }, )强制float32是因为海洋模型读驱动场时,大部分代码用单精度,float64会增加文件体积和 I/O 时间。zlib配合complevel=4是净收益比较高的组合,压缩率通常在 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(min=0, max=40) # 海水冰点经验公式: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(time="2020-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", engine="netcdf4", 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(min=0, max=1)如果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(lon=lon_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,而部分国产数据产品按北京时间(UTC+8)输出,decode_cf=True会按字符串里的基准来,不会自动加时区。
解决:先确认维度再决定要不要展开:
if ds.time.ndim == 0: ds = ds.expand_dims("time") print(ds.time.attrs.get("units"))如果是 UTC 基准但业务上需要北京时,手动加偏移并写清attrs:
ds = ds.assign_coords(time=ds.time + np.timedelta64(8, "h")) ds.time.attrs["units"] = "hours since 2020-01-01 00:00:00+08"这里最容易翻车的是你对“偏移 8 小时”的理解。先确认源数据文档里写的是 UTC 还是本地时,再决定加不加。盲目加偏移等于让模型初始时刻整体平移,对强强迫问题的结果影响不大,对海冰日变化敏感的问题会直接影响日循环相位。
5.5 现象五:岸线附近重插出“碎冰带”,海冰蔓延到内陆网格
现象:海冰浓度场在格陵兰岛东岸和加拿大北岸出现大量离散的碎冰格点,甚至越过海岸线出现在陆地上。
原因:插值器沿陆海边界处理时,把陆地缺失值当成 0 参与加权平均,或者conservative插值法把陆侧格点的面积权重也算进了海冰浓度。源数据的陆海掩膜和模型陆海掩膜如果分辨率不一致,问题更严重。
解决:海冰浓度建议优先用nearest_s2d,这类插值不会产生跨边界的中间值。如果必须用conservative,那就先把源数据的陆地格点设成NaN:
sic_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(time=slice("2020-07-01", "2020-09-30")) - oisst["sst"] print(diff.mean(dim=["lat", "lon"]).values) print(diff.sel(lat=slice(-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, sigma=1.5, mode="nearest") sst_smooth = xr.DataArray(sst_smooth_np, coords=ds_out["sst"].coords, attrs=ds_out["sst"].attrs) # 平滑会改变全场平均海温,简单加回一个常量做守恒修正 sst_smooth = sst_smooth + (ds_out["sst"].mean() - sst_smooth.mean())sigma=1.5意味着在 0.25 度网格上影响半径约 0.375 度,是一个很轻的平滑;如果网格更粗,比如 1 度,我会把sigma降到 1。mode="nearest"是为了避免边界外推产生异常值。海冰浓度不要做这种平滑,它的物理边界比温度更锐利,平滑之后阈值判定会失真。
这套方案我前后用了快两年,最大的体会是:海冰浓度永远比 SST 先处理,填缺失值永远比插值先处理。每次“看起来都正常”的输出,最后翻车都翻在那些没打印过的属性和没有核对过的坐标范围上。希望帮到你。
本文还有配套的精品资源,点击获取