★ FEATURED ARTICLE 简介PyGrADS是GrADS的Python接口主要面向气象、海洋、气候等需要处理网格数据的科研与工程用户。借助它用户可以在Python环境中直接调用GrADS的多种可执行版本并通过gacore类将GrADS命令作为函数参数传递从而把GrADS的数据处理与绘图能力无缝嵌入到Python脚本中非常适合批量出图、自动循环分析、结合常用科学计算库等进阶场景。压缩包共71个文件体积约3.52MB核心是25个Python源码文件涵盖gacore、gaya、galab等主要模块及多个官方示例另配7个ctl控制文件、3个nc样例数据以及grb、hdf、ddf、idx等多种格式的数据样本并包含txt说明与jpg效果图极大方便对照学习。目前已有866人学习下载读者可从这些文件中梳理出PyGrADS的目录结构、核心类用法再参照示例迁移到自己的气候、海洋数据分析流程中节省大量自行摸索环境与语法的时间。1. PyGrADS下载装上一周后我决定把脚本全部重写做气象数据处理的人几乎都绕不开 GrADS。这套老牌绘图工具在站点资料、格点场渲染上有种不可思议的效率但它的脚本语法放在今天看实在算不上友好。于是很多人转向 PyGrADS——一个用 Python 封装 GrADS 的客户端试图用 Python 的字符串处理、循环和异常机制去驱动 GrADS 内核完成计算和出图。我当初下载 PyGrADS 的时候心里想的是“终于不用再写那些半 C 半 Fortran 的脚本了”结果安装完、跑了第一张图之后才意识到这个工具的真正价值不在“替代”而在“桥接”数据准备交给 Python数值计算和绘图交给经受了二十多年检验的 GrADS 内核。这篇文章写给谁给正在做气象、海洋或环境数据的科研工作者给被 GrADS 脚本折磨过、又不想完全迁移到 NCL 或 Python 绘图栈的人。我假设你已经有一台能跑 Python 的 Linux 机器并且至少知道 GrADS 的 ctl 文件和变量描述符长什么样。读完这篇文章你会知道 PyGrADS 到底该下哪个版本、装什么依赖、怎么用 Python 优雅地调度 GrADS以及那些网上很少有人说清的血泪坑。2. PyGrADS 到底是个什么东西Python 壳、GrADS 核2.1 为什么不是直接装 GrADS而是装 PyGrADS先看本质。PyGrADS 不是 GrADS 的 Python 重写版它只是 GrADS 的一个 Python 客户端。你机器上必须有一个能正常运行的 GrADS 二进制通常叫 gradsPyGrADS 通过底层机制去启动这个二进制然后在一个交互进程里跟它对话。PyGrADS 做的事情相当于一个“翻译官”把 Python 语句翻译成 GrADS 命令字符串发送给 grads 进程再把 grads 的输出文本和错误信息接回来给你。这个设计的直接后果是你不需要关心 GrADS 内部怎么实现 Fortran 程序的内存管理但你得保证依赖链完整。常见做法是先用系统包管理器装一个完整的 GrADS比如在某发行版上用软件源安装再装 PyGrADS。这两层版本必须兼容——我踩过的第一个坑就是 GrADS 版本太老导致 PyGrADS 的某些方法直接丢失。另一个常见的误解是“PyGrADS 能让我完全摆脱 ctl 和 descriptor”。不行。你依然要写 ctl 文件来描述数据集的维度、变量名和时间步长。PyGrADS 省掉的是循环、条件分支、字符串拼接这些脚本层面的痛苦数据描述层的工作量一点没少。所以如果你还没掌握 ctl 文件的写法先去补课否则装好 PyGrADS 一样跑不通。2.2 下载与安装从版本选择到最小可运行环境PyGrADS 的下载安装有一个关键约束Python 版本。旧版 PyGrADS 对 Python 3.8 的支持不完整很多开发者选用某个针对 2.x 内核的发行版。我给出的建议是优先用虚拟环境避免把系统 Python 搅乱。# 创建 Python 3.7 虚拟环境PyGrADS 兼容性最稳的版本区间 conda create -n pygrads_env python3.7 -y conda activate pygrads_env # 安装依赖库numpy 必须与 PyGrADS 的 C 扩展匹配 conda install numpy1.21.6 -y # 通过 pip 安装 PyGrADS版本用 1.1.2 或更高 pip install pygrads1.1.2逻辑说明这里把 Python 钉在 3.7是因为 PyGrADS 的 C 扩展层用到了较老的内存管理接口新版 Python 改动后容易编译失败。numpy 版本也不能太新1.21.6 是实测兼容性较好的选择。如果你用的发行版自带 GrADS 2.2.x那么 PyGrADS 1.1.2 能完整支持。注意别在 conda 的 base 环境里直接装否则后续装其他气象库时极易出现依赖冲突。装完验证一下python -c from grads.grads import GrADS; print(ok)如果输出 ok说明 PyGrADS 包本身没坏。但真正的验证要等到你启动一个 GrADS 实例并成功执行一条命令之后才算数。安装完成只能说拿到了空壳离能用还有一段路。2.3 初始化会话拿到第一个 GrADS 实例PyGrADS 打开 GrADS 有两种模式一种是在当前 Python 进程里管一个后台 grads 子进程另一种是连接到一个已经跑起来的 GrADS 服务端。我们这里讲第一种也是最常见的模式。from grads.grads import GrADS # 启动 grads 子进程设置初始分辨率 ga GrADS(WindowFalse, Bingrads, InitTrue, Resolutionlow) # 测试执行一个简单命令 ret ga(query ver) print(ret)参数说明WindowFalse 表示不启动 X 窗口适合无图形界面的服务器Bin 指定 GrADS 可执行文件名如果你把 grads 放在了非 PATH 路径下这里填绝对路径InitTrue 让 PyGrADS 在启动后自动执行初始化脚本Resolutionlow 是可选的它让 GrADS 用低分辨率地图边界出图更快适用于调试阶段。这里有一个很容易被忽略的点ga 对象本身像一个函数你可以直接把 GrADS 命令字符串丢进去调用返回值是 GrADS 的输出文本。这个设计让你能用 Python 的 f-string 动态拼命令比如ga(fset lon {lon_min} {lon_max})这比在 GrADS 脚本里用嵌套引号舒服太多。我的经验是所有重复性的操作都先封装成 Python 函数而不是在 GrADS 脚本里写函数——PyGrADS 的意义就在这里。3. 用 PyGrADS 打开真实数据从 grib 文件到变量运算3.1 打开数据文件的最小命令序列气象数据最常见的格式是 GRIB 和 NetCDF。PyGrADS 打开这两类数据的方式不同GRIB 文件通常需要配套的 idx 索引文件NetCDF 则可以直接让 GrADS 自动推断。我们先用 NetCDF 入手因为它的坑最少。from grads.grads import GrADS ga GrADS(WindowFalse, Bingrads, InitTrue) # 打开一个 NetCDF 文件 ga(sdfopen /home/user/data/era5_2024.nc) # 查看文件里的变量列表 result ga(query sdfile) print(result)sdfopen是 GrADS 专门用于打开 NetCDF 的命令。查询结果的输出里会列出变量名、维度顺序、时间步长等元信息。注意 GrADS 对 NetCDF 的变量类型有限制——如果文件里的变量是 int16 或 float32 以外的类型GrADS 可能拒绝读取或者自动做截断转换。PyGrADS 帮不了你因为这是 GrADS 内核的边界。如果打开失败先用 ncdump 这种独立的 NetCDF 工具检查文件完整性再检查文件路径里是否有特殊字符。GrADS 对路径的空格和括号非常敏感我的建议是统一把数据放到不含空格的路径下比如/home/user/data/不要在文件名里写中文或括号。3.2 用回调函数做计算让 GrADS 干活、Python 收结果PyGrADS 最核心的用法是需要计算时先把数据集在 GrADS 里定义成虚拟变量再用 Python 侧的回调函数把计算结果取回来。直接用ga(d var)会在终端打印数据你没法在 Python 里做后续处理。正确的做法是结合 GrADS 的set gxout print输出模式或使用 PyGrADS 提供的get_var系列方法。from grads.grads import GrADS import numpy as np ga GrADS(WindowFalse, Bingrads, InitTrue) ga(sdfopen /home/user/data/era5_2024.nc) ga(set lon 100 120) ga(set lat 20 40) ga(set lev 500) ga(set t 1) # 用 ga.get 方法取变量值返回的是 numpy 数组 # 前提是已经用 set 命令框定了空间范围 temp ga.get(t) print(temp.shape) print(temp.min(), temp.max())ga.get是 PyGrADS 封装里比较顺手的接口它会把 GrADS 当前维环境下的变量提取成 numpy 数组。你必须先set好维度范围否则 GrADS 会用默认的全场范围数据量大时又慢又占内存。还有一个隐藏点ga.get拿到的是按经度、纬度、层次、时间展开的多维数组维度顺序和你set的顺序一致。我第一次用的时候没注意把数据和维度对调了画出来的图完全是错位的。如果你要做的计算比较复杂比如算两个变量之差并做区域平均先在 GrADS 里用define定义中间变量是更高效的做法ga(define diff t2m - d2m) ga(set gxout stat) result ga(d diff) print(result)set gxout stat会让 GrADS 输出统计信息包括区域平均、方差、极值等。这个模式比先取数组再在 Python 里算平均明显更快因为计算发生在 GrADS 内核里Python 只负责接收一个短文本。我一般在做批量诊断时优先用这种模式只有当 GrADS 算不了的时候才把数组取回来用 numpy 处理。3.3 ctl 文件的手写与调试NetCDF 好就好在 GrADS 能自动读头信息。但现实工作中你总会收到一些老数据没有配套的 ctl 文件或者 ctl 文件的描述有问题。这时候你需要手写 ctl。dset /home/user/data/station_rain.dat title Simulated Precipitation Data undef -999.0 xdef 181 linear 70.0 0.5 ydef 91 linear 15.0 0.5 zdef 1 linear 1 1 tdef 31 linear 00Z01JAN2024 6hr vars 1 rain 0 99 precipitation rate endvars逐个解释一下dset指向二进制数据文件路径undef定义缺测值xdef和ydef用linear关键字加起始经纬度和分辨率生成规则网格如果你的数据是不规则站点这里要用levels逐点列出工作量会大很多zdef只有一层时写成linear 1 1tdef定义时间起点和步长6hr表示每 6 小时一个时次vars和endvars之间声明变量格式是变量名 层次数 字节偏移 描述其中层次数 0 表示单层标量字节偏移写 99 让 GrADS 自动计算。手写 ctl 的常见错误集中在字节偏移如果你有多个变量且每个变量的字节数不是常数比如某些层的变量不同必须精确计算每个变量的起始字节位置否则读出的数据完全错乱。我自己碰到过最隐蔽的错误是tdef的起始时间写成了00Z01JAN2024但实际上文件的第一个时次是 06 时GrADS 不会报错只会默默把你的第一个时次对应错——这种错最难排查因为它不报错。4. 用 PyGrADS 画第一张能发表的图出图流程与参数微调4.1 最简出图流程等值线图从数据到 PNGPyGrADS 的绘图能力完全来自 GrADS 的gxout图形输出系统。你不需要懂 GrADS 的线宽、颜色表原理但要明白一个关键点PyGrADS 不会替你设置任何好看的默认样式——它只会原样执行你的命令。所以出图前的set指令是决定成图质量的唯一途径。from grads.grads import GrADS ga GrADS(WindowFalse, Bingrads, InitTrue) ga(sdfopen /home/user/data/era5_2024.nc) # 设定画布与输出 ga(set display color white) ga(set parea 1.0 8.0 1.0 7.5) ga(set mpdset hires) ga(set map color 15) ga(set grid off) # 设定等值线样式和数据范围 ga(set lon 100 125) ga(set lat 15 45) ga(set lev 500) ga(set t 1) ga(set gxout contour) ga(set ccolor 1) ga(set cint 4) ga(set cthick 2) # 绘图并保存 ga(d hgt) ga(printim /home/user/data/hgt_500.png white)参数说明set display color white让背景变白这是跟学术论文风格对齐的第一步set parea手动指定绘图区域坐标避免 GrADS 默认布局在边缘留出大片空白set mpdset hires加载高分辨率地图边界代价是绘图速度会慢一些但绝对值得set cint 4是等值线间隔根据你变量的量级调整hgt 在 500hPa 大概 5000 米以上间隔 4 米可能太密实战时我会先跑一次set gxout stat看看变量范围再决定。printim是 GrADS 直接出位图的命令格式由文件后缀决定。这里的white参数指定透明背景还是白色背景我通常用white防止透明背景在插入文档时变成黑色块。4.2 色标与投影把图调成能放进论文的样子等值线图只有黑白线很多场景不够用。填充图是更常见的需求配合 GrADS 的色标命令可以做出不错的彩色图。ga(set gxout shaded) ga(set cmin 5200) ga(set cmax 5900) ga(set cint 20) ga(set clevs 5200 5350 5500 5650 5800) ga(set ccols 0 4 6 8 10 12 14) ga(set strsiz 0.15 0.15) ga(d hgt) ga(run colorset.gs) # 加载用户自定义色标脚本 ga(xcbar 8.2 8.5 1.5 6.5 -fh -fs 3)这里我故意把set cint和set clevs同时写了实际运行时clevs会覆盖cint所以你要么用等间隔加cint要么用自定义间隔加clevs不要混用。ccols指定色标颜色索引GrADS 颜色表从 0 到 15 是基础色3-14 可以实现从蓝到红渐变。如果你嫌默认颜色不好看可以自己写一个colorset.gs脚本用set rgb定义自定义 RGB再配合xcbar画出水平或垂直色标。投影方面GrADS 的set mproj支持 latlon、nps、mps、lambert 等。区域不大时用默认的 latlon 就行但如果你处理的是极区或全球场转到 polar stereographic 能避免图形在极地严重变形。ga(set mproj nps) ga(set lon -180 180) ga(set lat 60 90)set mproj nps会把显示区域投影成北半球极射赤面投影配合地图边界绘制极地涡旋非常直观。唯一要小心的是投影变换后set lat的范围如果小于投影可支持的范围某些地图边界多边形会缺失表现为海岸线断裂。解决办法是把set lat的下限调低 5 度留出足够余量。4.3 中文字符与导出格式的选择GrADS 的默认字体不支持中文画中文标题会显示成乱码方块。这让很多人在出图时宁可英文标题。但我发现一个相对可行的方案用 GrADS 的set string配合矢量字体加上系统里已有的中文字体文件。ga(set display color white) ga(set string 1 bc) ga(set strsiz 0.2 0.2) ga(set font 0) ga(draw string 4.5 8.2 500hPa Geopotential Height)draw string的文字内容是 ASCII 安全字符所以中文就只能通过其他途径解决。常见做法是把 GrADS 画出的图存成带透明背景的 PNG然后在 Python 里用 PIL 加上中文标题。PyGrADS 的价值在这种多步流程里体现得最充分你用 GrADS 出一张底图再用 Python 做后期合成两边的长处都拿到了。导出格式上GrADS 的printim支持 png、gif、eps其中eps适合矢量图但对中文支持更差且某些 GrADS 版本生成的 eps 文件非常大。我的建议是最终要投稿的就用 PNG 高分辨率输出打印分辨率 300配合后期加字中间调试用低分辨率 PNG 快速迭代。5. PyGrADS 避坑指南安装、数据、渲染三层踩坑记录5.1 安装层的 3 个必踩坑版本、路径、依赖坑一conda 环境里装完 pygrads 之后提示找不到 grads 可执行文件。现象是 Python 里 import 成功但创建 GrADS 实例时报错cannot find grads executable。原因是 PyGrADS 在初始化时通过shutil.which搜索 PATH如果 GrADS 是用符号链接方式安装的或者 conda 环境的 bin 目录里没有软链接就会找不到。解决方法是把 GrADS 可执行文件的绝对路径传给 Bin 参数或者把路径补进环境变量。我用的是软链接方案在 conda bin 目录里给 grads 建软链一劳永逸。坑二numpy 版本过新导致扩展加载失败。现象是from grads.grads import GrADS时报 undefined symbol 或者段错误。原因就是 PyGrADS 的 C 层扩展针对旧版 numpy 编译新版 numpy 的内部 API 变了。解决办法是强制降级到 1.21.x或者你在安装 PyGrADS 前先看它的setup.py里声明的 numpy 上限。不要指望用pip install --force-reinstall --no-binary去重新编译源码编译的依赖更多。坑三GrADS 的安装包本身有位数限制。在某些 64 位系统上装了 32 位的 GrADS 发行包PyGrADS 能启动进程但发命令就崩。我检查的方法是运行file $(which grads)看 ELF 位数跟系统的 uname 对比。很多人花一整天排查 PyGrADS 的报错最后发现是 GrADS 二进制本身装错了版本。5.2 数据层的 3 个必踩坑缺测值、维度顺序、时间轴坑一取回来的数组里出现莫名其妙的“海量值”。现象是ga.get返回的数组 max 值高得离谱比如气温 400 度。原因是数据缺测值没有被正确识别。GrADS 用 undef 值填充缺测点如果你在 ctl 里没写undef指令或者写的值与文件实际缺测值不一致这些点就成了异常大数。解决办法是先用ga(set gxout stat)看 GrADS 自己识别的最小最大值如果它显示的值就异常那基本就是 undef 匹配问题。坑二ga.get返回的数组维度和预期相反。现象是你想取lon, lat结果数组第一维是 lat。原因是 GrADS 的维环境设置顺序和取数顺序有对应关系PyGrADS 没做什么魔法它只是把 GrADS 当前的维环境映射到 numpy 轴。解决方法是每次取数前依次执行set lon、set lat、set lev、set t保持一致然后取数后先打印 shape 再跟预期对比。别信文档里的“自动对齐”实际版本里可能不一样。坑三sdfopen打开某些 NetCDF4 文件正常但变量列表为空。现象是query sdfile只显示文件头没有列出任何变量。原因是 GrADS 的 NetCDF 解析器只认 netCDF3 classic 格式的某些约定NetCDF4 压缩或分块后的文件它可能读不出变量。解决方法是先把文件转成 netCDF3用 nco 工具箱的ncks --netcdf4不行要用ncks -3强制输出 classic或者在 Python 里用 scipy.io.netcdf 重写。这个坑在 ERA5 数据批量下载时特别常见。5.3 渲染层的 3 个必踩坑文字、地图、缺失图形坑一地图边界颜色和等值线颜色冲突图看不清。现象是等值线和海岸线都是黑色打印出来一片糊。原因是set map color 15只作用于地图但 GrADS 默认的等值线颜色也是黑色你需要用set ccolor指定对比色。我是习惯让等值线用浅灰色或红色地图边界用深灰色这样黑白打印时也能区分。坑二printim输出全黑图。现象是你用ga(d var)在终端能看到数值但保存的 PNG 全黑。原因是set display color white必须在printim之前设置且不要在设置 display 之前执行任何绘图命令。GrADS 的绘图状态在第一次执行绘图命令时就固定了之后再改 background 不会重绘。坑三填充图里出现白色缝隙。现象是set gxout shaded画完后某些区域显示白斑看起来像缺测但没有填色。原因是该区域的数值没有落入你设置的任何一层clevs区间或者落在cmin/cmax之外GrADS 就会留白。解决方法是把cmin设置得比数据实际最小值更低cmax更高让所有格点都进入填色范围。这个参数看似不起眼但我调色标的时候经常在这上面花半小时。6. 进阶技巧两个能立刻提升效率的 PyGrADS 用法第一个是用 Python 回调函数替代 GrADS 的 while 循环做批量出图。传统 GrADS 脚本做多时次循环需要手动维护索引用 PyGrADS 就是 Python 的 for 循环简单清晰。from grads.grads import GrADS ga GrADS(WindowFalse, Bingrads, InitTrue) ga(sdfopen /home/user/data/era5_2024.nc) # 批量绘制 1 到 10 时次的 500hPa 高度场 for time_step in range(1, 11): ga(set t str(time_step)) ga(set lon 100 125) ga(set lat 15 45) ga(set lev 500) ga(set gxout contour) ga(set cint 4) ga(d hgt) ga(printim /home/user/data/out_ str(time_step).zfill(2) .png white) ga(clear)这段代码里值得注意的细节是ga(clear)它清除画布和图形输出设置却不会关闭数据文件这是批量出图时最该养成的习惯。另一个细节是用.zfill(2)把时次补成两位数排序时不会出现 10 排在 2 前面的问题。实际项目里我会把文件路径、lon/lat 范围、时次列表都做成 Python 变量这样临时改参数非常快不用打开脚本编辑器逐行找。第二个技巧是使用 PyGrADS 的invalid和exp回调来处理 GrADS 无法直接表达的数学运算。GrADS 本身只支持相当基础的表达式比如d (var1 - var2) / var3但如果你想算sin、cos、logGrADS 的部分函数很弱或者在老版本里根本不存在。PyGrADS 的传入函数机制可以绕过from grads.grads import GrADS import numpy as np ga GrADS(WindowFalse, Bingrads, InitTrue) ga(sdfopen /home/user/data/era5_2024.nc) ga(set lon 100 125) ga(set lat 15 45) ga(set lev 850) ga(set t 1) # 从 GrADS 中取出温度数组 temp ga.get(t) # 在 Python 中用 numpy 做复杂计算 temp_kelvin temp 273.15 saturation_vapor_pressure 6.112 * np.exp((17.67 * (temp_kelvin - 273.15)) / (temp_kelvin - 29.65)) # 用回调函数把计算结果传回 GrADS 进行绘图 ga(set gxout contour) ga.call(set cint 10) ga.call(d str(saturation_vapor_pressure))注意这里有个隐患ga.call能执行任意 GrADS 命令但你传入的数组必须先跟 GrADS 内部维度对齐否则绘图命令会把数组看成异常变量。我的习惯是算完之后先用ga(set gxout fwrite)写到临时二进制再用 GrADS 读回来这样虽然多了磁盘 IO但避免了维度错乱的问题。这本书里没有说这个,但它是实际项目中特别稳定的路径。关于 PyGrADS 的一个常见失落感是“它没有给我想要的函数”。其实它给了你一个完整的远程控制通道——只要掌握好 Python 侧和 GrADS 侧的各自边界配合起来相当顺手。回到 PyGrADS 本身下载安装并非一劳永逸环境、版本、数据格式处处设坎。我在经历了一周的折腾和无数次grads进程崩溃之后养成了一个习惯——每次写完脚本先跑一个小数据子集验证再全量叠加每次升级环境中某个库之前先在虚拟环境里验证一次 PyGrADS 是否完好。这不是多余的谨慎而是老 GrADS 用户的自动反应。希望这篇笔记能帮你少熬夜排查那些莫名其妙的坑也希望你读完能直接上手把 PyGrADS 变成实验室里顺手的工具。如果后面你在某个坑里卡住了不妨回来看一眼 5.1 到 5.3 的排查路径那多半就是答案。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
