首页
学习
活动
专区
圈层
工具
发布
社区首页 >专栏 >代码实战 | 使用 Herbie 库获取 GFS 数据

代码实战 | 使用 Herbie 库获取 GFS 数据

作者头像
用户11172986
发布2026-04-24 19:30:45
发布2026-04-24 19:30:45
460
举报
文章被收录于专栏:气python风雨气python风雨

代码实战 | 使用 Herbie 库获取 GFS 数据

问题场景

在进行气象数据分析和可视化时,获取高质量的数值预报数据是第一步也是最重要的一步。GFS(Global Forecast System)是美国国家环境预报中心(NCEP)运行的全球数值天气预报模型,提供全球范围的高分辨率预报数据。

然而,直接从 NOAA 服务器下载 GFS 数据面临以下挑战:

  • 数据源众多:AWS、Google Cloud、Azure、NOMADS 等多个数据源
  • 文件格式复杂:GRIB2 格式需要专门的解码工具
  • 数据量大:完整预报文件通常数百 MB,下载耗时
  • 索引文件:需要处理 .idx 索引文件来定位特定变量

技术方案

Herbie 是一个专门为气象数据获取设计的 Python 库,它完美解决了上述问题:

  • 自动搜索多个数据源,选择最快的可用源
  • 支持 GRIB2 格式数据的自动解码
  • 可以下载完整文件或按变量子集下载,节省带宽
  • 与 xarray 无缝集成,直接返回可分析的数据结构
  • 支持历史数据和实时预报数据

本文将详细介绍如何使用 Herbie 库获取 GFS 数据的完整流程。


环境准备

代码语言:javascript
复制
# 安装必要的库
!pip install herbie-data xarray cfgrib matplotlib cartopy -q

1. 基础用法:获取 GFS 数据

首先导入 Herbie 库并创建一个 Herbie 对象来获取指定日期的 GFS 数据。

代码语言:javascript
复制
from herbie import Herbie
import matplotlib.pyplot as plt

# 创建Herbie对象获取2025年1月20日的GFS数据
# model='gfs' 指定使用GFS模型
# product='pgrb2.0p25' 指定使用0.25度分辨率的公共字段产品
H = Herbie("2025-01-20", model="gfs", product="pgrb2.0p25")

print(f"找到的数据源: {H.SOURCES}")
print(f"\n可用产品: {H.PRODUCTS}")
代码语言:javascript
复制
/opt/conda/lib/python3.11/site-packages/gribapi/__init__.py:23: UserWarning: ecCodes 2.42.0 or higher is recommended. You are running version 2.29.0
  warnings.warn(

✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
找到的数据源: {'aws': 'https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000', 'nomads': 'https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000', 'google': 'https://storage.googleapis.com/global-forecast-system/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000', 'azure': 'https://noaagfs.blob.core.windows.net/gfs/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000', 'ncar_rda': 'https://data.rda.ucar.edu/d084001/2025/20250120/gfs.0p25.2025012000.f000.grib2'}

可用产品: {'pgrb2.0p25': 'common fields, 0.25 degree resolution', 'pgrb2.0p50': 'common fields, 0.50 degree resolution', 'pgrb2.1p00': 'common fields, 1.00 degree resolution', 'pgrb2b.0p25': 'uncommon fields, 0.25 degree resolution', 'pgrb2b.0p50': 'uncommon fields, 0.50 degree resolution', 'pgrb2b.1p00': 'uncommon fields, 1.00 degree resolution', 'pgrb2full.0p50': 'combined grids of 0.50 resolution', 'sfluxgrb': 'surface flux fields, T1534 Semi-Lagrangian grid', 'goesimpgrb2.0p25': ', 0.50 degree resolution'}

2. 查看数据内容

使用 inventory() 方法可以查看 GRIB2 文件中包含的所有变量。

代码语言:javascript
复制
# 查看所有可用变量
inventory = H.inventory()
print(f"总共有 {len(inventory)} 个变量")
print("\n前10个变量:")
print(inventory.head(10))
代码语言:javascript
复制
Downloading inventory file from self.idx='https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000.idx'
总共有 696 个变量

前10个变量:
   grib_message  start_byte   end_byte            range reference_time  \
0             1           0  1037405.0        0-1037405     2025-01-20
1             2     1037406  1120810.0  1037406-1120810     2025-01-20
2             3     1120811  1441033.0  1120811-1441033     2025-01-20
3             4     1441034  1655963.0  1441034-1655963     2025-01-20
4             5     1655964  1748833.0  1655964-1748833     2025-01-20
5             6     1748834  1794026.0  1748834-1794026     2025-01-20
6             7     1794027  2544142.0  1794027-2544142     2025-01-20
7             8     2544143  3294995.0  2544143-3294995     2025-01-20
8             9     3294996  4118960.0  3294996-4118960     2025-01-20
9            10     4118961  4865850.0  4118961-4865850     2025-01-20

  valid_time variable              level forecast_time  \
0 2025-01-20    PRMSL     mean sea level           anl
1 2025-01-20     CLMR     1 hybrid level           anl
2 2025-01-20     ICMR     1 hybrid level           anl
3 2025-01-20     RWMR     1 hybrid level           anl
4 2025-01-20     SNMR     1 hybrid level           anl
5 2025-01-20     GRLE     1 hybrid level           anl
6 2025-01-20     REFD     1 hybrid level           anl
7 2025-01-20     REFD     2 hybrid level           anl
8 2025-01-20     REFC  entire atmosphere           anl
9 2025-01-20      VIS            surface           anl

                   search_this
0    :PRMSL:mean sea level:anl
1     :CLMR:1 hybrid level:anl
2     :ICMR:1 hybrid level:anl
3     :RWMR:1 hybrid level:anl
4     :SNMR:1 hybrid level:anl
5     :GRLE:1 hybrid level:anl
6     :REFD:1 hybrid level:anl
7     :REFD:2 hybrid level:anl
8  :REFC:entire atmosphere:anl
9             :VIS:surface:anl

3. 按变量筛选数据

可以通过变量名筛选查看特定变量。

代码语言:javascript
复制
# 查看所有温度相关的变量
temp_vars = H.inventory("TMP")
print("温度相关变量:")
print(temp_vars)
代码语言:javascript
复制
温度相关变量:
     grib_message  start_byte     end_byte                range  \
15             16     7898420    8481263.0      7898420-8481263
25             26    15349681   16179307.0    15349681-16179307
35             36    23521109   24408730.0    23521109-24408730
45             46    31648088   32531711.0    31648088-32531711
55             56    39676875   40525934.0    39676875-40525934
65             66    47702988   48530386.0    47702988-48530386
75             76    55335665   56156267.0    55335665-56156267
85             86    62811677   63621500.0    62811677-63621500
95             96    70585870   71385192.0    70585870-71385192
105           106    77988476   78805651.0    77988476-78805651
115           116    85547631   86339817.0    85547631-86339817
125           126    93115800   93909741.0    93115800-93909741
135           136   100841452  101611595.0  100841452-101611595
145           146   108534989  109288588.0  108534989-109288588
155           156   116359924  117110627.0  116359924-117110627
165           166   124295367  125043401.0  124295367-125043401
175           176   132197134  132926676.0  132197134-132926676
185           186   140332556  141089177.0  140332556-141089177
195           196   149202321  149944562.0  149202321-149944562
211           212   161392352  162151414.0  161392352-162151414
221           222   170417515  171183673.0  170417515-171183673
237           238   181208357  181960078.0  181208357-181960078
253           254   191385689  192146110.0  191385689-192146110
269           270   201321723  202069264.0  201321723-202069264
285           286   211573654  212337832.0  211573654-212337832
301           302   222473645  223209024.0  222473645-223209024
317           318   233114570  233848115.0  233114570-233848115
333           334   243573396  244300695.0  243573396-244300695
349           350   253967193  254692060.0  253967193-254692060
365           366   264830999  265559433.0  264830999-265559433
381           382   275838132  276570737.0  275838132-276570737
397           398   286547234  287286504.0  286547234-287286504
413           414   297466771  298213588.0  297466771-298213588
429           430   308433498  309206502.0  308433498-309206502
445           446   319593674  320401662.0  319593674-320401662
461           462   331084948  331925983.0  331084948-331925983
477           478   343073010  343927424.0  343073010-343927424
493           494   355166245  356021972.0  355166245-356021972
509           510   367181428  368031908.0  367181428-368031908
526           527   379061027  379907405.0  379061027-379907405
541           542   389311162  390162811.0  389311162-390162811
562           563   403038770  403607852.0  403038770-403607852
579           580   410223309  411089041.0  410223309-411089041
583           584   413583689  414139320.0  413583689-414139320
618           619   435271822  436227251.0  435271822-436227251
627           628   444552556  445795573.0  444552556-445795573
636           637   453482050  454327197.0  453482050-454327197
641           642   458262516  459108662.0  458262516-459108662
644           645   461041540  461929340.0  461041540-461929340
647           648   463926179  464785750.0  463926179-464785750
650           651   466784040  467642193.0  466784040-467642193
657           658   473188548  474022851.0  473188548-474022851
670           671   483523818  484387626.0  483523818-484387626
683           684   491197285  491375166.0  491197285-491375166
686           687   492649622  493281520.0  492649622-493281520
692           693   497304595  497952385.0  497304595-497952385

    reference_time valid_time variable                          level  \
15      2025-01-20 2025-01-20      TMP                        0.01 mb
25      2025-01-20 2025-01-20      TMP                        0.02 mb
35      2025-01-20 2025-01-20      TMP                        0.04 mb
45      2025-01-20 2025-01-20      TMP                        0.07 mb
55      2025-01-20 2025-01-20      TMP                         0.1 mb
65      2025-01-20 2025-01-20      TMP                         0.2 mb
75      2025-01-20 2025-01-20      TMP                         0.4 mb
85      2025-01-20 2025-01-20      TMP                         0.7 mb
95      2025-01-20 2025-01-20      TMP                           1 mb
105     2025-01-20 2025-01-20      TMP                           2 mb
115     2025-01-20 2025-01-20      TMP                           3 mb
125     2025-01-20 2025-01-20      TMP                           5 mb
135     2025-01-20 2025-01-20      TMP                           7 mb
145     2025-01-20 2025-01-20      TMP                          10 mb
155     2025-01-20 2025-01-20      TMP                          15 mb
165     2025-01-20 2025-01-20      TMP                          20 mb
175     2025-01-20 2025-01-20      TMP                          30 mb
185     2025-01-20 2025-01-20      TMP                          40 mb
195     2025-01-20 2025-01-20      TMP                          50 mb
211     2025-01-20 2025-01-20      TMP                          70 mb
221     2025-01-20 2025-01-20      TMP                         100 mb
237     2025-01-20 2025-01-20      TMP                         150 mb
253     2025-01-20 2025-01-20      TMP                         200 mb
269     2025-01-20 2025-01-20      TMP                         250 mb
285     2025-01-20 2025-01-20      TMP                         300 mb
301     2025-01-20 2025-01-20      TMP                         350 mb
317     2025-01-20 2025-01-20      TMP                         400 mb
333     2025-01-20 2025-01-20      TMP                         450 mb
349     2025-01-20 2025-01-20      TMP                         500 mb
365     2025-01-20 2025-01-20      TMP                         550 mb
381     2025-01-20 2025-01-20      TMP                         600 mb
397     2025-01-20 2025-01-20      TMP                         650 mb
413     2025-01-20 2025-01-20      TMP                         700 mb
429     2025-01-20 2025-01-20      TMP                         750 mb
445     2025-01-20 2025-01-20      TMP                         800 mb
461     2025-01-20 2025-01-20      TMP                         850 mb
477     2025-01-20 2025-01-20      TMP                         900 mb
493     2025-01-20 2025-01-20      TMP                         925 mb
509     2025-01-20 2025-01-20      TMP                         950 mb
526     2025-01-20 2025-01-20      TMP                         975 mb
541     2025-01-20 2025-01-20      TMP                        1000 mb
562     2025-01-20 2025-01-20      TMP                        surface
579     2025-01-20 2025-01-20      TMP               2 m above ground
583     2025-01-20 2025-01-20    APTMP               2 m above ground
618     2025-01-20 2025-01-20      TMP                     tropopause
627     2025-01-20 2025-01-20      TMP                       max wind
636     2025-01-20 2025-01-20      TMP              80 m above ground
641     2025-01-20 2025-01-20      TMP             100 m above ground
644     2025-01-20 2025-01-20      TMP    1829 m above mean sea level
647     2025-01-20 2025-01-20      TMP    2743 m above mean sea level
650     2025-01-20 2025-01-20      TMP    3658 m above mean sea level
657     2025-01-20 2025-01-20      TMP           30-0 mb above ground
670     2025-01-20 2025-01-20      TMP              0.995 sigma level
683     2025-01-20 2025-01-20   ICETMP                        surface
686     2025-01-20 2025-01-20      TMP   PV=2e-06 (Km^2/kg/s) surface
692     2025-01-20 2025-01-20      TMP  PV=-2e-06 (Km^2/kg/s) surface

    forecast_time                             search_this
15            anl                        :TMP:0.01 mb:anl
25            anl                        :TMP:0.02 mb:anl
35            anl                        :TMP:0.04 mb:anl
45            anl                        :TMP:0.07 mb:anl
55            anl                         :TMP:0.1 mb:anl
65            anl                         :TMP:0.2 mb:anl
75            anl                         :TMP:0.4 mb:anl
85            anl                         :TMP:0.7 mb:anl
95            anl                           :TMP:1 mb:anl
105           anl                           :TMP:2 mb:anl
115           anl                           :TMP:3 mb:anl
125           anl                           :TMP:5 mb:anl
135           anl                           :TMP:7 mb:anl
145           anl                          :TMP:10 mb:anl
155           anl                          :TMP:15 mb:anl
165           anl                          :TMP:20 mb:anl
175           anl                          :TMP:30 mb:anl
185           anl                          :TMP:40 mb:anl
195           anl                          :TMP:50 mb:anl
211           anl                          :TMP:70 mb:anl
221           anl                         :TMP:100 mb:anl
237           anl                         :TMP:150 mb:anl
253           anl                         :TMP:200 mb:anl
269           anl                         :TMP:250 mb:anl
285           anl                         :TMP:300 mb:anl
301           anl                         :TMP:350 mb:anl
317           anl                         :TMP:400 mb:anl
333           anl                         :TMP:450 mb:anl
349           anl                         :TMP:500 mb:anl
365           anl                         :TMP:550 mb:anl
381           anl                         :TMP:600 mb:anl
397           anl                         :TMP:650 mb:anl
413           anl                         :TMP:700 mb:anl
429           anl                         :TMP:750 mb:anl
445           anl                         :TMP:800 mb:anl
461           anl                         :TMP:850 mb:anl
477           anl                         :TMP:900 mb:anl
493           anl                         :TMP:925 mb:anl
509           anl                         :TMP:950 mb:anl
526           anl                         :TMP:975 mb:anl
541           anl                        :TMP:1000 mb:anl
562           anl                        :TMP:surface:anl
579           anl               :TMP:2 m above ground:anl
583           anl             :APTMP:2 m above ground:anl
618           anl                     :TMP:tropopause:anl
627           anl                       :TMP:max wind:anl
636           anl              :TMP:80 m above ground:anl
641           anl             :TMP:100 m above ground:anl
644           anl    :TMP:1829 m above mean sea level:anl
647           anl    :TMP:2743 m above mean sea level:anl
650           anl    :TMP:3658 m above mean sea level:anl
657           anl           :TMP:30-0 mb above ground:anl
670           anl              :TMP:0.995 sigma level:anl
683           anl                     :ICETMP:surface:anl
686           anl   :TMP:PV=2e-06 (Km^2/kg/s) surface:anl
692           anl  :TMP:PV=-2e-06 (Km^2/kg/s) surface:anl

4. 下载并读取数据

使用 xarray() 方法下载并读取特定变量到 xarray Dataset 中。

代码语言:javascript
复制
# 下载并读取2米高度的温度数据
# ":TMP:2 m above" 是GRIB2变量的搜索字符串
ds = H.xarray(":TMP:2 m above")

print("数据集信息:")
print(ds)
print(f"\n数据形状: {ds.t2m.shape}")
print(f"温度范围: {ds.t2m.min().values:.2f} - {ds.t2m.max().values:.2f} K")
print(f"温度范围(摄氏度): {ds.t2m.min().values - 273.15:.2f} - {ds.t2m.max().values - 273.15:.2f} °C")
代码语言:javascript
复制
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'cfradial1' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'datamet' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'furuno' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'gamic' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'hpl' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'iris' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'metek' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'nexradlevel2' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'odim' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'radolan' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)
/opt/conda/lib/python3.11/site-packages/xarray/backends/plugins.py:109: RuntimeWarning: Engine 'rainbow' loading failed:
No module named 'xarray.core.merge'
  external_backend_entrypoints = backends_dict_from_pkg(entrypoints_unique)

数据集信息:
<xarray.Dataset> Size: 4MB
Dimensions:              (latitude: 721, longitude: 1440)
Coordinates:
  * latitude             (latitude) float64 6kB 90.0 89.75 89.5 ... -89.75 -90.0
  * longitude            (longitude) float64 12kB 0.0 0.25 0.5 ... 359.5 359.8
    time                 datetime64[ns] 8B 2025-01-20
    step                 timedelta64[ns] 8B 00:00:00
    heightAboveGround    float64 8B 2.0
    valid_time           datetime64[ns] 8B 2025-01-20
    gribfile_projection  object 8B None
Data variables:
    t2m                  (latitude, longitude) float32 4MB 242.4 242.4 ... 253.5
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    model:                   gfs
    product:                 pgrb2.0p25
    description:             NOAA Global Forecast System (GFS)
    remote_grib:             https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20...
    local_grib:              /home/mw/data/gfs/20250120/subset_f7ef22b0__gfs....
    search:                  :TMP:2 m above

数据形状: (721, 1440)
温度范围: 221.08 - 310.62 K
温度范围(摄氏度): -52.07 - 37.47 °C

5. 可视化数据

使用 matplotlib 和 cartopy 绘制温度分布图。

代码语言:javascript
复制
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

# 创建地图投影
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())

# 添加地图特征
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.OCEAN, color='lightblue', alpha=0.5)
ax.add_feature(cfeature.LAND, color='lightgray', alpha=0.5)

# 绘制温度数据(转换为摄氏度)
temp_c = ds.t2m - 273.15
im = ax.pcolormesh(
    ds.longitude,
    ds.latitude,
    temp_c,
    transform=ccrs.PlateCarree(),
    cmap='RdYlBu_r',
    vmin=-40,
    vmax=40,
    shading='auto'
)

# 添加颜色条
cbar = plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05, aspect=50)
cbar.set_label('温度 (°C)', fontsize=12)

# 设置标题
ax.set_title(f"GFS 2米温度预报\n{ds.t2m.time.values}", fontsize=14, fontweight='bold')
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

plt.tight_layout()
plt.show()
代码语言:javascript
复制
<Figure size 1200x800 with 2 Axes>

GFS 2米温度预报

6. 获取多个预报时效数据

GFS 提供不同预报时效的数据,可以指定 fxx 参数获取特定时效的预报。

代码语言:javascript
复制
# 获取6小时预报
H_f6 = Herbie("2025-01-20 00:00", model="gfs", product="pgrb2.0p25", fxx=6)
ds_f6 = H_f6.xarray(":TMP:2 m above")

# 获取12小时预报
H_f12 = Herbie("2025-01-20 00:00", model="gfs", product="pgrb2.0p25", fxx=12)
ds_f12 = H_f12.xarray(":TMP:2 m above")

# 获取24小时预报
H_f24 = Herbie("2025-01-20 00:00", model="gfs", product="pgrb2.0p25", fxx=24)
ds_f24 = H_f24.xarray(":TMP:2 m above")

print(f"6小时预报时间: {ds_f6.t2m.time.values}")
print(f"12小时预报时间: {ds_f12.t2m.time.values}")
print(f"24小时预报时间: {ds_f24.t2m.time.values}")
代码语言:javascript
复制
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F06[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
Downloading inventory file from self.idx='https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f006.idx'
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F12[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
Downloading inventory file from self.idx='https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f012.idx'
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F24[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
Downloading inventory file from self.idx='https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f024.idx'
6小时预报时间: 2025-01-20T00:00:00.000000000
12小时预报时间: 2025-01-20T00:00:00.000000000
24小时预报时间: 2025-01-20T00:00:00.000000000

7. 对比不同时效的预报

对比分析场预报和 6 小时预报的温度差异。

代码语言:javascript
复制
# 计算温度变化
temp_change = ds_f6.t2m - ds.t2m

# 绘制温度变化图
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())

ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)

im = ax.pcolormesh(
    ds_f6.longitude,
    ds_f6.latitude,
    temp_change,
    transform=ccrs.PlateCarree(),
    cmap='RdBu_r',
    vmin=-10,
    vmax=10,
    shading='auto'
)

cbar = plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05, aspect=50)
cbar.set_label('温度变化 (K)', fontsize=12)

ax.set_title(f"GFS 6小时温度变化预报\n{ds.t2m.time.values} → {ds_f6.t2m.time.values}", fontsize=14, fontweight='bold')
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

plt.tight_layout()
plt.show()
代码语言:javascript
复制
<Figure size 1200x800 with 2 Axes>

GFS 6小时温度变化预报

8. 获取多个变量

可以同时获取多个气象变量进行分析。

代码语言:javascript
复制
# 获取2米温度和相对湿度
ds_temp = H.xarray(":TMP:2 m above")
ds_rh = H.xarray(":RH:2 m above")

# 创建子图
fig, axes = plt.subplots(1, 2, figsize=(16, 6), subplot_kw={'projection': ccrs.PlateCarree()})

# 温度图
axes[0].add_feature(cfeature.COASTLINE)
axes[0].add_feature(cfeature.BORDERS)
im1 = axes[0].pcolormesh(
    ds_temp.longitude,
    ds_temp.latitude,
    ds_temp.t2m - 273.15,
    transform=ccrs.PlateCarree(),
    cmap='RdYlBu_r',
    vmin=-40,
    vmax=40,
    shading='auto'
)
cbar1 = plt.colorbar(im1, ax=axes[0], orientation='horizontal', pad=0.05)
cbar1.set_label('温度 (°C)')
axes[0].set_title('2米温度', fontweight='bold')
axes[0].set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

# 相对湿度图
axes[1].add_feature(cfeature.COASTLINE)
axes[1].add_feature(cfeature.BORDERS)
im2 = axes[1].pcolormesh(
    ds_rh.longitude,
    ds_rh.latitude,
    ds_rh.r2,
    transform=ccrs.PlateCarree(),
    cmap='YlGnBu',
    vmin=0,
    vmax=100,
    shading='auto'
)
cbar2 = plt.colorbar(im2, ax=axes[1], orientation='horizontal', pad=0.05)
cbar2.set_label('相对湿度 (%)')
axes[1].set_title('2米相对湿度', fontweight='bold')
axes[1].set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

plt.suptitle(f"GFS气象要素预报\n{ds_temp.t2m.time.values}", fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
代码语言:javascript
复制
<Figure size 1600x600 with 4 Axes>

GFS气象要素预报

9. 保存处理后的数据

可以将处理后的 xarray Dataset 保存为 NetCDF 格式,便于后续分析。

代码语言:javascript
复制
# 创建输出目录
output_dir = "./gfs_data"
os.makedirs(output_dir, exist_ok=True)

# 保存为NetCDF文件
output_file = f"{output_dir}/gfs_20250120_t2m.nc"
ds.to_netcdf(output_file)

print(f"数据已保存到: {output_file}")
print(f"文件大小: {os.path.getsize(output_file) / (1024 * 1024):.2f} MB")

# 验证保存的数据
ds_loaded = H.xarray(":TMP:2 m above")
print(f"\n重新加载的数据形状: {ds_loaded.t2m.shape}")
print(f"数据一致性检查: {np.allclose(ds.t2m.values, ds_loaded.t2m.values)}")
代码语言:javascript
复制
数据已保存到: ./gfs_data/gfs_20250120_t2m.nc
文件大小: 4.00 MB

重新加载的数据形状: (721, 1440)
数据一致性检查: True

10. 指定数据源

如果需要使用特定的数据源,可以通过 overwrite 参数指定。

代码语言:javascript
复制
# 从AWS数据源获取
H_aws = Herbie(
    "2025-01-20",
    model="gfs",
    product="pgrb2.0p25",
    overwrite=True,
    save_dir="./gfs_aws"
)

print(f"AWS数据源: {H_aws.SOURCES.get('aws', 'Not available')}")

# 从Google Cloud获取
H_google = Herbie(
    "2025-01-20",
    model="gfs",
    product="pgrb2.0p25",
    overwrite=True,
    save_dir="./gfs_google"
)

print(f"Google Cloud数据源: {H_google.SOURCES.get('google', 'Not available')}")
代码语言:javascript
复制
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
AWS数据源: https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ aws[0m
Google Cloud数据源: https://storage.googleapis.com/global-forecast-system/gfs.20250120/00/atmos/gfs.t00z.pgrb2.0p25.f000

11. 获取历史数据

Herbie 支持获取历史 GFS 数据,对于 2021 年之前的数据会自动从 NCEI 或 RDA 存档获取。

代码语言:javascript
复制
# 获取2020年的历史数据
H_history = Herbie("2020-06-01", model="gfs")

print(f"历史数据源: {H_history.SOURCES}")
print(f"可用产品: {H_history.PRODUCTS}")

# 查看历史数据的内容
inventory_history = H_history.inventory("GRD:80 m above")
print("\n80米高度风场变量:")
print(inventory_history)
代码语言:javascript
复制
✅ Found ┊ model=gfs ┊ [3mproduct=0.5-degree[0m ┊ [38;2;41;130;13m2020-Jun-01 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ ncei_analysis[0m ┊ [38;2;255;153;0m[3mIDX @ ncei_analysis[0m
历史数据源: {'ncei_analysis': 'https://www.ncei.noaa.gov/data/global-forecast-system/access/grid-004-0.5-degree/analysis/202006/20200601/gfs_4_20200601_0000_000.grb2', 'ncei_forecast': 'https://www.ncei.noaa.gov/data/global-forecast-system/access/grid-004-0.5-degree/forecast/202006/20200601/gfs_4_20200601_0000_000.grb2', 'ncar_rda': 'https://data.rda.ucar.edu/d084001/2020/20200601/gfs.0p25.2020060100.f000.grib2', 'ncei_historical_analysis': 'https://www.ncei.noaa.gov/data/global-forecast-system/access/historical/analysis/202006/20200601/gfsanl_4_20200601_0000_000.grb2'}
可用产品: {'0.5-degree': '0.5 degree grid', '1.0-degree': '1.0 degree grid'}
Downloading inventory file from self.idx='https://www.ncei.noaa.gov/data/global-forecast-system/access/grid-004-0.5-degree/analysis/202006/20200601/gfs_4_20200601_0000_000.grb2.inv'

80米高度风场变量:
     grib_message  start_byte    end_byte              range reference_time  \
465           466    80208994  80496003.0  80208994-80496003     2020-06-01
466           467    80496004  80780823.0  80496004-80780823     2020-06-01

    valid_time variable              level forecast_time  \
465 2020-06-01     UGRD  80 m above ground           anl
466 2020-06-01     VGRD  80 m above ground           anl

                     search_this
465  :UGRD:80 m above ground:anl
466  :VGRD:80 m above ground:anl

12. 批量获取数据

演示如何批量获取多个时间点的数据。

代码语言:javascript
复制
from datetime import datetime, timedelta

# 生成时间序列
start_date = datetime(2025, 1, 20)
end_date = datetime(2025, 1, 22)
dates = [start_date + timedelta(days=i) for i in range(3)]

# 批量获取数据
datasets = []
for date in dates:
    print(f"正在获取 {date.strftime('%Y-%m-%d')} 的数据...")
    H_batch = Herbie(date.strftime('%Y-%m-%d'), model="gfs", product="pgrb2.0p25")
    ds_batch = H_batch.xarray(":TMP:2 m above")
    datasets.append(ds_batch)

print(f"\n成功获取 {len(datasets)} 个时间点的数据")
for i, ds in enumerate(datasets):
    print(f"  {i+1}. {ds.t2m.time.values}")
代码语言:javascript
复制
正在获取 2025-01-20 的数据...
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-20 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ local[0m
正在获取 2025-01-21 的数据...
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-21 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ local[0m
正在获取 2025-01-22 的数据...
✅ Found ┊ model=gfs ┊ [3mproduct=pgrb2.0p25[0m ┊ [38;2;41;130;13m2025-Jan-22 00:00 UTC[92m F00[0m ┊ [38;2;255;153;0m[3mGRIB2 @ aws[0m ┊ [38;2;255;153;0m[3mIDX @ local[0m

成功获取 3 个时间点的数据
  1. 2025-01-20T00:00:00.000000000
  2. 2025-01-21T00:00:00.000000000
  3. 2025-01-22T00:00:00.000000000

13. 区域数据提取

演示如何提取特定区域的数据进行分析。

代码语言:javascript
复制
# 提取中国区域的数据
lon_min, lon_max = 70, 140
lat_min, lat_max = 10, 55

# 使用xarray的sel方法选择区域
ds_china = ds.sel(
    longitude=slice(lon_min, lon_max),
    latitude=slice(lat_max, lat_min)  # 注意latitude顺序
)

print(f"中国区域数据形状: {ds_china.t2m.shape}")
print(f"中国区域平均温度: {(ds_china.t2m.mean() - 273.15).values:.2f} °C")

# 绘制中国区域温度图
fig = plt.figure(figsize=(10, 8))
ax = plt.axes(projection=ccrs.PlateCarree())

# ax.add_feature(cfeature.COASTLINE)
# ax.add_feature(cfeature.BORDERS)

from meteva.base.tool.plot_tools import add_china_map_2basemap
add_china_map_2basemap(ax, name="river", edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "河流"
add_china_map_2basemap(ax, name="nation",edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "国界"
add_china_map_2basemap(ax, name="province", edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "省界"
add_china_map_2basemap(ax, name="world", edgecolor='k', lw=0.5, encoding='gbk', grid0=None)  # "省界"
ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

im = ax.pcolormesh(
    ds_china.longitude,
    ds_china.latitude,
    ds_china.t2m - 273.15,
    transform=ccrs.PlateCarree(),
    cmap='RdYlBu_r',
    vmin=-20,
    vmax=40,
    shading='auto'
)

cbar = plt.colorbar(im, ax=ax, orientation='horizontal', pad=0.05, aspect=50)
cbar.set_label('温度 (°C)', fontsize=12)
ax.set_title(f"中国区域GFS 2米温度预报\n{ds_china.t2m.time.values}", fontsize=14, fontweight='bold')

plt.tight_layout()
plt.show()
代码语言:javascript
复制
中国区域数据形状: (181, 281)
中国区域平均温度: 3.45 °C

<Figure size 1000x800 with 2 Axes>

中国区域GFS 2米温度预报

小结

通过本文的学习,你应该掌握了:

  1. 安装和配置 Herbie 库
  2. 获取 GFS 实时和历史数据
  3. 查看和筛选数据变量
  4. 将数据读取为 xarray Dataset
  5. 数据可视化和地图绘制
  6. 多时效预报数据获取和对比
  7. 多变量数据获取
  8. 数据保存为 NetCDF 格式
  9. 指定数据源
  10. 批量数据获取
  11. 区域数据提取

Herbie 库大大简化了气象数据的获取流程,是气象数据分析和可视化工作的利器。

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2026-03-16,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气python风雨 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 代码实战 | 使用 Herbie 库获取 GFS 数据
    • 问题场景
    • 技术方案
    • 环境准备
    • 1. 基础用法:获取 GFS 数据
    • 2. 查看数据内容
    • 3. 按变量筛选数据
    • 4. 下载并读取数据
    • 5. 可视化数据
    • 6. 获取多个预报时效数据
    • 7. 对比不同时效的预报
    • 8. 获取多个变量
    • 9. 保存处理后的数据
    • 10. 指定数据源
    • 11. 获取历史数据
    • 12. 批量获取数据
    • 13. 区域数据提取
    • 小结
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档