在进行气象数据分析和可视化时,获取高质量的数值预报数据是第一步也是最重要的一步。GFS(Global Forecast System)是美国国家环境预报中心(NCEP)运行的全球数值天气预报模型,提供全球范围的高分辨率预报数据。
然而,直接从 NOAA 服务器下载 GFS 数据面临以下挑战:
Herbie 是一个专门为气象数据获取设计的 Python 库,它完美解决了上述问题:
本文将详细介绍如何使用 Herbie 库获取 GFS 数据的完整流程。
# 安装必要的库
!pip install herbie-data xarray cfgrib matplotlib cartopy -q
首先导入 Herbie 库并创建一个 Herbie 对象来获取指定日期的 GFS 数据。
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}")
/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'}
使用 inventory() 方法可以查看 GRIB2 文件中包含的所有变量。
# 查看所有可用变量
inventory = H.inventory()
print(f"总共有 {len(inventory)} 个变量")
print("\n前10个变量:")
print(inventory.head(10))
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
可以通过变量名筛选查看特定变量。
# 查看所有温度相关的变量
temp_vars = H.inventory("TMP")
print("温度相关变量:")
print(temp_vars)
温度相关变量:
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
使用 xarray() 方法下载并读取特定变量到 xarray Dataset 中。
# 下载并读取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")
/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
使用 matplotlib 和 cartopy 绘制温度分布图。
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()
<Figure size 1200x800 with 2 Axes>
GFS 2米温度预报
GFS 提供不同预报时效的数据,可以指定 fxx 参数获取特定时效的预报。
# 获取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}")
✅ 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
对比分析场预报和 6 小时预报的温度差异。
# 计算温度变化
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()
<Figure size 1200x800 with 2 Axes>
GFS 6小时温度变化预报
可以同时获取多个气象变量进行分析。
# 获取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()
<Figure size 1600x600 with 4 Axes>
GFS气象要素预报
可以将处理后的 xarray Dataset 保存为 NetCDF 格式,便于后续分析。
# 创建输出目录
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)}")
数据已保存到: ./gfs_data/gfs_20250120_t2m.nc
文件大小: 4.00 MB
重新加载的数据形状: (721, 1440)
数据一致性检查: True
如果需要使用特定的数据源,可以通过 overwrite 参数指定。
# 从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')}")
✅ 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
Herbie 支持获取历史 GFS 数据,对于 2021 年之前的数据会自动从 NCEI 或 RDA 存档获取。
# 获取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)
✅ 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
演示如何批量获取多个时间点的数据。
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}")
正在获取 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
演示如何提取特定区域的数据进行分析。
# 提取中国区域的数据
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()
中国区域数据形状: (181, 281)
中国区域平均温度: 3.45 °C
<Figure size 1000x800 with 2 Axes>
中国区域GFS 2米温度预报
通过本文的学习,你应该掌握了:
Herbie 库大大简化了气象数据的获取流程,是气象数据分析和可视化工作的利器。