1. 从数据文件到第一张图为什么是xarray和netCDF如果你刚开始接触气象、海洋或者地球科学的数据分析打开数据文件夹看到一堆以.nc结尾的文件可能会有点懵。这些就是 netCDF 文件是这个领域的“普通话”。我第一次处理这类数据时也试图用传统的pandas去读取结果要么报错要么读出来的数据结构完全不对路。后来才明白处理这种自带维度、坐标和多变量的科学数据你需要一个更趁手的工具——xarray。简单来说netCDFNetwork Common Data Form是一种为存储多维科学数据如温度、气压、风速随经度、纬度、高度和时间变化而设计的文件格式。它不仅仅是存数字还把数据的维度信息比如经度数组、纬度数组、时间序列、变量属性单位、长名称以及全局属性数据来源、作者都打包在一起。这就好比一个精心整理的仓库不仅货物数据值摆放整齐还有完整的货架标签维度和货物说明书属性。而xarray就是为操作这种“带标签的多维数组”而生的 Python 库。它深受pandas的影响但将核心数据结构从一维的Series和二维的DataFrame扩展到了 N 维的DataArray和Dataset。一个DataArray可以看作是一个带有坐标标签的多维numpy数组而一个Dataset则是多个共享相同坐标的DataArray的集合。这种设计让基于维度名称的选择、切片、计算和绘图变得异常直观。所以这个系列的第一篇我们就从最基础的环节开始如何用xarray打开一个 netCDF 文件并快速画出第一张科学图表。这个过程远不止是open加plot那么简单里面涉及到数据结构的理解、坐标系的处理以及如何避免初次绘图时常见的“坑”。我会结合一个真实的数据示例把每一步为什么这么做讲清楚。2. 环境准备与数据获取搭建你的分析工作台在写第一行代码之前我们需要一个合适的环境。我强烈建议使用conda来管理科学计算的环境它能很好地处理各个库特别是涉及地理信息处理的库之间的依赖关系。2.1 创建并激活专用环境打开你的终端或 Anaconda Prompt执行以下命令来创建一个名为geo_viz的新环境并安装核心包# 创建新环境指定Python版本如3.9 conda create -n geo_viz python3.9 # 激活环境 conda activate geo_viz # 安装核心数据分析与可视化库 conda install -c conda-forge xarray netcdf4 dask # 安装绘图库。cartopy是地理绘图神器但安装稍复杂我们先装基础的。 conda install -c conda-forge matplotlib numpy这里有几个关键点为什么用conda-forge通道conda-forge社区维护的软件包通常更新更快版本兼容性更好特别是对于xarray,cartopy这类科学计算栈的包。netcdf4是必须的xarray本身不直接读取文件它需要后端引擎。netcdf4是读取 netCDF 文件最常用、最稳定的引擎。dask是可选的但推荐对于大型数组比如全球高分辨率气候模式数据dask可以实现惰性计算和并行处理避免内存爆炸。我们在初学阶段可能用不到它的高级功能但先装上无害。如果你后续需要做地理投影的绘图比如把数据画到地图上那么cartopy是必不可少的。但它依赖PROJ和GEOS等地理空间库在 Windows 上通过pip安装容易出错。因此最稳妥的方式还是通过conda安装conda install -c conda-forge cartopy安装过程可能会慢一些请耐心等待。2.2 获取示例数据为了演示我们需要一个 netCDF 文件。你可以从许多气候数据中心免费下载比如NASA GES DISC、ECMWF或NOAA。这里我推荐一个对新手非常友好的来源Pangeo Gallery的示例数据集。我们使用xarray内置的教程数据集或者从一个公开的 OPeNDAP 服务器加载一个小型数据集。但为了完全模拟本地文件操作我们可以下载一个小的示例文件。这里我们使用pooch库来下载一个示例文件首先安装它conda install -c conda-forge pooch然后在 Python 脚本或 Jupyter Notebook 中我们可以这样获取一个全球海表温度SST的示例数据import pooch import xarray as xr # 定义一个文件下载器 downloader pooch.create( pathpooch.os_cache(xarray_tutorial), base_urlhttps://github.com/pangeo-data/weather-bench/raw/main/data/, registry{ sst_1990_2018_1.40625deg.nc: md5:abc123..., # 此处应为真实MD5仅为示例 } ) # 由于真实MD5校验较复杂我们换一种更直接的方式使用xarray自带的示例数据集 # 但xarray的tutorial数据集可能不在本地我们改用一个更稳定的在线小数据集 # 例如读取一个ERA5再分析数据的单时间片样例通过OPeNDAP try: # 这是一个公开的测试数据集URL可能随时间失效但常用于演示 data_url https://psl.noaa.gov/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/air.sig995.1948.nc ds xr.open_dataset(data_url) print(成功从网络加载示例数据集。) except Exception as e: print(f网络加载失败: {e}) print(我们将使用xarray内置的示例数据集。) # 加载xarray内置的一个非常小的示例数据集 ds xr.tutorial.open_dataset(rasm)为了本教程的稳定性和可重复性我们直接使用xarray内置的rasm数据集。这是一个模拟的北极区域气候模式输出包含了多个变量如气温、降水以及经纬度坐标。import xarray as xr import matplotlib.pyplot as plt # 加载内置示例数据集 ds xr.tutorial.open_dataset(rasm) print(ds)运行上述代码你会看到类似下面的输出它展示了xarray.Dataset的结构xarray.Dataset Dimensions: (time: 36, y: 205, x: 275) Coordinates: * time (time) object 1980-09-16 12:00:00 ... 1983-08-17 12:00:00 * x (x) float64 0.0 1.0 2.0 3.0 4.0 ... 271.0 272.0 273.0 274.0 * y (y) float64 0.0 1.0 2.0 3.0 4.0 ... 200.0 201.0 202.0 203.0 204.0 Data variables: Tair (time, y, x) float64 ... Attributes: title: example Rasm data institution: U.W. source: RACM Rasm example output_frequency: daily ...这个输出就是理解 netCDF 数据的钥匙。它告诉我们Dimensions维度: 数据有三个维度time36个时间点、y205个点、x275个点。这通常对应着时间、纬度、经度。Coordinates坐标: 每个维度对应的具体坐标值。例如time坐标是36个日期时间对象。Data variables数据变量: 这里显示了一个变量Tair近地表气温它的形状是(time, y, x)数据类型是浮点数。Attributes属性: 数据集级别的全局信息如标题、机构等。3. 解剖数据集理解xarray的数据结构在动手画图之前我们必须弄清楚手里数据的“长相”。盲目画图很可能得到一张坐标轴奇怪、单位不明的图片。3.1 探索数据集内容让我们深入查看一下这个数据集# 查看所有变量 print(ds.data_vars) # 查看单个变量例如气温 Tair tair ds[Tair] print(tair) print(f变量单位: {tair.attrs.get(units, N/A)}) print(f变量长名: {tair.attrs.get(long_name, N/A)}) # 查看坐标 print(ds.coords) # 查看全局属性 print(ds.attrs)关键点在于attrs属性。科学数据中units单位和long_name长名称是绘图的灵魂它们会自动被用于美化坐标轴标签。例如Tair的单位可能是K开尔文长名称可能是Air temperature。3.2 理解坐标与投影细心的你可能发现了这个数据集的坐标叫x和y而不是longitude和latitude。这在区域气候模式输出中很常见数据可能存储在一个特定的地图投影如 Lambert Conformal Conic的笛卡尔网格上。真正的经纬度信息通常作为附加的二维坐标变量存在。# 检查是否有经纬度坐标 if lon in ds.coords or longitude in ds.coords: print(找到经度坐标) if lat in ds.coords or latitude in ds.coords: print(找到纬度坐标) # 更常见的是经纬度是依赖于 (y, x) 的二维变量 print(ds[xc]) print(ds[yc])你会发现ds[xc]和ds[yc]也是DataArray它们的维度是(y, x)。这意味着网格上每个(y, x)点都对应一个经度值xc和一个纬度值yc。这种坐标称为“二维笛卡尔坐标”或“投影坐标”。在绘图时我们需要用xc和yc来定位数据而不是用x和y的索引值。3.3 选择与切片数据我们很少需要一次性画出所有时间和空间的数据。xarray的强大之处在于其基于标签的索引。# 选择第一个时间点 tair_first ds[Tair].isel(time0) # isel 代表 integer select用索引选择 # 或者用更直观的标签选择如果坐标是datetime类型 # tair_first ds[Tair].sel(time1980-09-16) # 选择空间子区域例如y索引从50到150x索引从100到200 tair_sub tair_first.isel(yslice(50, 150), xslice(100, 200)) # 基于实际经纬度范围选择这需要二维坐标 # 假设我们想选取经度-160到-140纬度60到75的大致范围 # 我们需要找到对应的索引这通常需要一点技巧因为坐标是二维的。 # 一个近似的方法是 lon_min, lon_max -160, -140 lat_min, lat_max 60, 75 mask_lon (ds[xc] lon_min) (ds[xc] lon_max) mask_lat (ds[yc] lat_min) (ds[yc] lat_max) mask mask_lon mask_lat # 这个mask是二维布尔数组要应用到数据上需要结合.where()方法 tair_region tair_first.where(mask, dropTrue) # dropTrue会丢弃所有为NaN的点sel和isel是xarray中最常用的两个方法务必熟练掌握。sel用于按坐标值选择isel用于按整数索引选择。4. 绘制第一张科学图表从简单到实用现在数据已经在手理解也到位了是时候生成可视化结果了。4.1 最基础的快速绘图xarray.DataArray自带一个.plot()方法它基于matplotlib对于快速查看数据极其方便。# 绘制第一个时间点的气温场 tair_first ds[Tair].isel(time0) tair_first.plot() plt.show() # 在脚本中需要调用show()来显示图形就这么简单一行.plot()一张带有颜色条colorbar、自动使用变量名和单位作为标签的图片就生成了。但是你可能会发现横纵坐标是x和y而不是我们习惯的经纬度。4.2 使用二维坐标进行地理绘图为了画出更符合认知的地图我们需要告诉绘图函数每个数据点对应的经纬度是什么。这可以通过x和y参数来实现。fig, ax plt.subplots(figsize(10, 8)) # 使用 pcolormesh 绘制并传入二维的经纬度坐标 # 注意.T 是转置因为通常我们期望经度是x轴纬度是y轴。 # 但这里数据存储的维度顺序是 (y, x)而 pcolormesh 要求 (X, Y) 是二维且维度与数据一致。 # 更安全的做法是直接使用 xarray 的 .plot.pcolormesh mesh tair_first.plot.pcolormesh(xxc, yyc, axax, add_colorbarTrue) # 添加一些美化 ax.set_title(fSurface Air Temperature - {str(tair_first.time.values)[:10]}) # xarray 通常会自动从属性中设置标签但我们也可以手动覆盖 ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.tight_layout() plt.show()这里的关键参数是xxc和yyc。.plot.pcolormesh方法会自动从数据集ds中寻找名为xc和yc的坐标数组并用它们来定位每个网格单元。这样画出来的图坐标轴就变成了有实际意义的经纬度。注意.plot()默认使用.plot.imshow()它假设数据是均匀网格绘图速度很快。.plot.pcolormesh()则能处理非均匀网格比如我们的二维投影坐标并且能正确绘制每个网格单元但速度稍慢。对于规则经纬度网格两者区别不大对于投影坐标或非均匀网格必须使用pcolormesh。4.3 处理时间序列与多子图我们经常需要查看数据随时间的变化或者对比不同时间点的空间格局。# 选择前4个时间点 tair_subset ds[Tair].isel(timeslice(0, 4)) # 创建一个2x2的子图 fig, axes plt.subplots(nrows2, ncols2, figsize(14, 10)) axes axes.flatten() # 将二维的axes数组展平成一维方便循环 # 循环绘制每个时间点 for i, time_idx in enumerate(range(4)): tair_slice tair_subset.isel(timetime_idx) # 在每个子图上绘图 mesh tair_slice.plot.pcolormesh(axaxes[i], xxc, yyc, add_colorbarTrue, cmapRdBu_r) # 设置子图标题为时间 axes[i].set_title(fTime: {str(tair_slice.time.values)[:10]}) plt.tight_layout() plt.show()这里我们引入了cmapRdBu_r参数来更改颜色映射。RdBu_r是一个红蓝渐变色常用于表示温度异常冷蓝热红_r表示反转色带。选择合适的色带对于准确传达科学信息至关重要。4.4 自定义颜色条和图形美化默认的绘图可能不符合论文或报告的要求我们需要进行深度定制。import matplotlib.pyplot as plt import numpy as np tair_first ds[Tair].isel(time0) fig, ax plt.subplots(figsize(12, 8)) # 绘图并返回一个“可映射对象”的集合包含图像、颜色条等 # 我们暂时不自动添加颜色条 im tair_first.plot.pcolormesh(axax, xxc, yyc, add_colorbarFalse, cmapviridis, robustTrue) # 手动添加颜色条并控制其位置和大小 cbar plt.colorbar(im, axax, orientationvertical, pad0.03, shrink0.8) # 设置颜色条标签优先使用变量的 long_name 和 units cbar_label f{tair_first.attrs.get(long_name, Tair)} [{tair_first.attrs.get(units, )}] cbar.set_label(cbar_label, fontsize12) # 设置标题和坐标轴标签 ax.set_title(Regional Climate Model Output - Near Surface Air Temperature, fontsize14, fontweightbold) ax.set_xlabel(Longitude (degrees_east), fontsize12) ax.set_ylabel(Latitude (degrees_north), fontsize12) # 设置坐标轴刻度密度和格式 ax.xaxis.set_major_locator(plt.MaxNLocator(6)) ax.yaxis.set_major_locator(plt.MaxNLocator(6)) # 添加网格线虚线浅灰色 ax.grid(True, linestyle--, linewidth0.5, alpha0.7, colorgray) # 自动调整布局防止标签重叠 plt.tight_layout() # 保存图形到文件设置高DPI以保证印刷质量 plt.savefig(first_plot.png, dpi300, bbox_inchestight) plt.show()这段代码中包含了几个重要的技巧robustTrue这个参数非常实用。它会自动根据数据的第2和第98百分位数而非最小最大值来设定颜色条的范围。这能有效避免个别异常值如缺省值、极端值把整个颜色条范围拉得很宽导致主要数据区域的颜色对比度丧失。手动控制颜色条通过plt.colorbar()手动添加可以更灵活地控制颜色条的位置 (pad)、大小 (shrink)、方向等。动态生成标签使用.attrs.get(long_name, Tair)来获取属性如果属性不存在则提供一个默认值。这使代码对不同的数据集更具鲁棒性。图形保存plt.savefig()的bbox_inchestight参数可以自动裁剪图形周围的空白区域dpi300确保输出图片有足够的分辨率。5. 实战中的常见“坑”与解决之道纸上得来终觉浅绝知此事要躬行。在实际操作中你会遇到各种预料之外的情况。下面是我踩过的一些坑和总结的经验。5.1 坑一内存不足与懒加载netCDF 文件尤其是再分析数据或气候模式输出动辄几十GB。直接用xr.open_dataset(huge_file.nc)可能会瞬间耗尽内存。xarray的默认行为是“急加载”eager loading即打开时就把所有变量的数据读入内存。解决方案懒加载Lazy Loading使用chunks参数进行分块结合dask后端。# 使用 chunks 参数打开大型数据集 # ‘auto’ 让 dask 自动决定块大小也可以指定如 chunks{time: 10, lat: 100, lon: 100} ds_lazy xr.open_dataset(large_file.nc, chunksauto) print(ds_lazy)这时ds_lazy中的DataArray存储的不再是numpy数组而是dask.array。只有当你真正需要数据时如调用.compute()或进行.plot()等触发计算的操作dask才会按块读取和计算。这让你可以在普通笔记本电脑上处理远超内存大小的数据。注意事项懒加载状态下一些操作如.values会触发计算。在开发阶段可以先对数据的一个子集如ds_lazy.isel(time0)进行计算和绘图确认无误后再处理全量数据。5.2 坑二时间坐标的混乱处理科学数据中的时间坐标格式千奇百怪有的是“days since 1900-01-01”有的是“seconds since 1970-01-01 00:00:00”还有的用整数表示年月如202301。xarray在打开 netCDF 时如果识别出units属性符合 CF 公约标准会自动将其解码为datetime对象。但并非所有数据都那么规范。# 检查时间坐标类型 print(ds.time.dtype) print(ds.time.attrs) # 如果时间坐标是整数或浮点数且带有‘units’和‘calendar’属性xarray通常能自动解码。 # 如果解码失败可以手动转换 import pandas as pd import cftime # 假设时间坐标是‘days since 1950-01-01’ # ds[time] 可能还是数值类型 if not np.issubdtype(ds[time].dtype, np.datetime64): try: # 使用 cftime 库进行解码处理非标准日历如360_day times cftime.num2date(ds[time].values, unitsds[time].attrs[units], calendards[time].attrs.get(calendar, standard)) # 将解码后的时间赋值回数据集 ds[time] (time, pd.DatetimeIndex(times)) # 转换为 pandas DatetimeIndex except Exception as e: print(f自动解码时间失败: {e}) # 备选方案使用 pandas 手动解析如果知道格式的话 # ds[time] pd.to_datetime(ds[time].astype(str), format%Y%m)处理时间坐标是科学数据可视化的基石错误的时间会导致序列分析、气候态计算等全部出错。务必在绘图前确认ds.time的数据类型是datetime64[ns]或类似的日期时间类型。5.3 坑三缺失值与无效值的处理netCDF 数据中常用一个特定的填充值如-9999.0,1e20来表示缺失值或陆地/海洋掩膜。xarray在读取时如果检测到_FillValue或missing_value属性会自动将这些值替换为NaN。# 检查变量是否有填充值属性 print(tair.attrs) # 如果数据没有自动处理或者你想用不同的阈值可以手动处理 # 例如将所有小于 -100 的值设为 NaN tair_clean tair.where(tair -100) # 在绘图时NaN 区域默认是透明的。如果你想用特定颜色填充如陆地用灰色可以使用 .plot() 的 add_colorbar 和 extend 参数并配合设置 bad color cmap plt.cm.viridis cmap.set_bad(gray, 1.0) # 将 NaN 区域设置为灰色 tair_clean.plot(cmapcmap)5.4 坑四投影转换与地图绘制我们的示例数据rasm使用的是投影坐标。如果你想把它画到常见的地理底图上比如带有海岸线、国界线的地图就需要用到cartopy。这是一个比单纯传入x、y坐标更高级但也更复杂的话题。import cartopy.crs as ccrs import cartopy.feature as cfeature # 1. 创建带有地理投影的图形和轴 # 我们需要知道原始数据的投影。rasm数据使用了一个旋转极地立体投影但为简化我们假设它是普通经纬度。 # 实际上对于未知投影的数据用其原生坐标绘图如前所述是最安全的。 # 这里演示如何将数据画到 PlateCarree普通经纬度投影的地图上。 proj ccrs.PlateCarree() # 目标投影 fig, ax plt.subplots(figsize(12, 8), subplot_kw{projection: proj}) # 2. 添加地理特征 ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS, linewidth0.5) # 3. 绘图。关键必须指定 transform 参数告诉 cartopy 数据本身的坐标系统。 # 如果数据是经纬度transformccrs.PlateCarree()。 # 如果数据是其他投影需要找到对应的 cartopy.crs 对象。 # 对于 rasm 数据其 xc, yc 是投影坐标但这里我们不知道具体投影所以这是一个难点。 # 一个常见做法如果数据文件提供了投影信息如 grid_mapping 属性可以据此创建 crs。 # 本例中我们退一步仅用其二维坐标绘图不叠加地理要素。 # 假设我们强行将其视为经纬度仅用于演示不科学 im tair_first.plot.pcolormesh(axax, xxc, yyc, transformccrs.PlateCarree(), add_colorbarFalse, cmapviridis, robustTrue) # 4. 设置图形范围 ax.set_extent([-180, 180, 40, 90], crsccrs.PlateCarree()) # 大致北极区域 plt.colorbar(im, axax, orientationvertical, shrink0.8) ax.set_title(Air Temperature with Cartopy Basemap) plt.show()使用cartopy的核心是理解“数据投影”transform参数和“地图投影”subplot_kw{projection: ...}参数的区别。数据必须从其原生投影转换到地图显示的投影上cartopy会在内部完成这个插值计算。如果数据本身没有明确定义的投影信息很多 netCDF 文件确实如此那么使用cartopy绘制地理底图就会非常棘手。在这种情况下像我们之前那样直接用二维坐标(xc, yc)绘制填充图并手动添加经纬度刻度往往是更实际的选择。6. 代码封装与可重复工作流当你掌握了基本流程后应该考虑将常用操作封装成函数以提高效率并减少错误。def plot_netcdf_var(nc_file, var_name, time_idx0, level_idxNone, regionNone, output_pathNone): 绘制netCDF文件中指定变量的空间分布图。 参数 ---------- nc_file : str netCDF文件路径。 var_name : str 要绘制的变量名。 time_idx : int, 可选 时间维度的索引默认为0。 level_idx : int, 可选 层次维度的索引如气压层默认为None不选择。 region : dict, 可选 空间区域限制格式为 {lon_min:, lon_max:, lat_min:, lat_max:}。 需要数据有标准的 lon/lat 一维坐标。 output_path : str, 可选 图片保存路径默认为None不保存。 返回 ------- fig, ax : matplotlib图形和坐标轴对象 # 1. 打开数据集 ds xr.open_dataset(nc_file) # 2. 选择变量和数据切片 da ds[var_name] if time in da.dims: da da.isel(timetime_idx) if level_idx is not None and level in da.dims: da da.isel(levellevel_idx) # 3. 区域选择简化版仅适用于规则经纬度网格 if region: # 假设坐标名是 lon 和 lat da da.sel(lonslice(region[lon_min], region[lon_max]), latslice(region[lat_min], region[lat_max])) # 4. 创建图形 fig, ax plt.subplots(figsize(12, 8)) # 5. 绘图 - 尝试使用经纬度坐标如果存在的话 plot_kwargs {ax: ax, robust: True, add_colorbar: True} if lon in da.coords and lat in da.coords: # 规则经纬度网格 im da.plot.pcolormesh(xlon, ylat, **plot_kwargs) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) elif xc in da.coords and yc in da.coords: # 投影坐标网格 im da.plot.pcolormesh(xxc, yyc, **plot_kwargs) ax.set_xlabel(X Coordinate) ax.set_ylabel(Y Coordinate) else: # 回退到索引 im da.plot.pcolormesh(**plot_kwargs) # 6. 美化标题 title f{var_name} if long_name in da.attrs: title da.attrs[long_name] if units in da.attrs: title f [{da.attrs[units]}] if time in da.coords: title f\nTime: {str(da.time.values)[:19]} ax.set_title(title, fontsize14) # 7. 保存图形 if output_path: plt.savefig(output_path, dpi300, bbox_inchestight) print(f图形已保存至: {output_path}) return fig, ax # 使用函数 fig, ax plot_netcdf_var(your_data.nc, temperature, time_idx5, output_pathtemp_plot.png) plt.show()这个函数只是一个起点你可以根据需求扩展它比如添加对更多维度高度、成员的支持、集成cartopy地图、自定义颜色映射等。将重复性工作函数化是迈向高效数据分析的关键一步。从双击一个陌生的.nc文件到生成一张信息丰富、可用于分析的图表这个过程涉及了对数据格式的理解、对工具的熟练运用以及对科学可视化原则的把握。xarray极大地简化了 netCDF 数据的操作但其真正的威力在于将数据维度、坐标和属性作为一等公民对待的思维方式。记住在画图之前花时间print(ds)和探索ds.coords、ds.attrs永远是值得的。它不仅能帮你避免很多低级错误还能让你更深刻地理解手中的数据所讲述的科学故事。在接下来的系列文章中我们会深入更多专题比如时间序列分析、垂直剖面绘制、多变量对比以及更高级的地图定制技巧。