GEFS

This notebook demonstrates using the Global Ensemble Forecast System. It is available on AWS for 2017 to present.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

from herbie import Herbie
from paint.standard2 import cm_tmp, cm_wind, cm_wave_height, cm_pcp
from toolbox import EasyMap, pc

The GEFS model product could be any of the following: - 'atmos.5' - Half degree atmos PRIMARY fields (pgrb2ap5); ~83 most common variables. - 'atmos.5b' - Half degree atmos SECONDARY fields (pgrb2bp5); ~500 least common variables - 'atmos.25' - Quarter degree atmos PRIMARY fields (pgrb2sp25); ~35 most common variables - 'wave' - Global wave products. - 'chem.5' - Chemistry fields on 0.5 degree grid - 'chem.25' - Chemistry fields on 0.25 degree grid

You also must specify a member. For the atmos output, this should be something like 0 or "c00" for control member, 1-30 or "p01"-"p30" for members 1-30, or 'avg' or 'mean' for the ensemble mean, 'spr' for ensemble spread.

H = Herbie(
    "2023-01-04 12:00",
    model="gefs",
    product="atmos.5",
    member="mean",
)
H

✅ Found ┊ model=gefs ┊ product=atmos.5 ┊ 2023-Jan-04 12:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws

▌▌Herbie GEFS model atmos.5 product initialized 2023-Jan-04 12:00 UTC F00 ┊ source=aws

# Show all the available products (from the model template file)
H.PRODUCTS

{'atmos.5': 'Half degree atmos PRIMARY fields (pgrb2ap5); ~83 most common variables.',
 'atmos.5b': 'Half degree atmos SECONDARY fields (pgrb2bp5); ~500 least common variables',
 'atmos.25': 'Quarter degree atmos PRIMARY fields (pgrb2sp25); ~35 most common variables',
 'wave': 'Global wave products.',
 'chem.5': 'Chemistry fields on 0.5 degree grid',
 'chem.25': 'Chemistry fields on 0.25 degree grid'}

df = H.inventory()
df

	grib_message	start_byte	end_byte	range	reference_time	valid_time	variable	level	forecast_time	?	search_this
0	1	0	224944	0-224944	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	10 mb	anl	ens mean	:HGT:10 mb:anl:ens mean
1	2	224944	353986	224944-353986	2023-01-04 12:00:00	2023-01-04 12:00:00	TMP	10 mb	anl	ens mean	:TMP:10 mb:anl:ens mean
2	3	353986	387271	353986-387271	2023-01-04 12:00:00	2023-01-04 12:00:00	RH	10 mb	anl	ens mean	:RH:10 mb:anl:ens mean
3	4	387271	641627	387271-641627	2023-01-04 12:00:00	2023-01-04 12:00:00	UGRD	10 mb	anl	ens mean	:UGRD:10 mb:anl:ens mean
4	5	641627	881486	641627-881486	2023-01-04 12:00:00	2023-01-04 12:00:00	VGRD	10 mb	anl	ens mean	:VGRD:10 mb:anl:ens mean
...	...	...	...	...	...	...	...	...	...	...	...
66	67	13100932	13361479	13100932-13361479	2023-01-04 12:00:00	2023-01-04 12:00:00	VGRD	10 m above ground	anl	ens mean	:VGRD:10 m above ground:anl:ens mean
67	68	13361479	13520800	13361479-13520800	2023-01-04 12:00:00	2023-01-04 12:00:00	PWAT	entire atmosphere (considered as a single layer)	anl	ens mean	:PWAT:entire atmosphere (considered as a singl...
68	69	13520800	13665977	13520800-13665977	2023-01-04 12:00:00	2023-01-04 12:00:00	CAPE	180-0 mb above ground	anl	ens mean	:CAPE:180-0 mb above ground:anl:ens mean
69	70	13665977	13762089	13665977-13762089	2023-01-04 12:00:00	2023-01-04 12:00:00	CIN	180-0 mb above ground	anl	ens mean	:CIN:180-0 mb above ground:anl:ens mean
70	71	13762089		13762089-	2023-01-04 12:00:00	2023-01-04 12:00:00	PRMSL	mean sea level	anl	ens mean	:PRMSL:mean sea level:anl:ens mean

71 rows × 11 columns

df.variable.unique()

array(['HGT', 'TMP', 'RH', 'UGRD', 'VGRD', 'VVEL', 'PRES', 'TSOIL',
       'SOILW', 'WEASD', 'SNOD', 'ICETK', 'PWAT', 'CAPE', 'CIN', 'PRMSL'],
      dtype=object)

H.inventory("HGT")

	grib_message	start_byte	end_byte	range	reference_time	valid_time	variable	level	forecast_time	?	search_this
0	1	0	224944	0-224944	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	10 mb	anl	ens mean	:HGT:10 mb:anl:ens mean
5	6	881486	1113634	881486-1113634	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	50 mb	anl	ens mean	:HGT:50 mb:anl:ens mean
10	11	1829912	2058354	1829912-2058354	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	100 mb	anl	ens mean	:HGT:100 mb:anl:ens mean
15	16	2853920	3078454	2853920-3078454	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	200 mb	anl	ens mean	:HGT:200 mb:anl:ens mean
20	21	3750171	3972303	3750171-3972303	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	250 mb	anl	ens mean	:HGT:250 mb:anl:ens mean
25	26	4665076	4885223	4665076-4885223	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	300 mb	anl	ens mean	:HGT:300 mb:anl:ens mean
30	31	5575673	5819613	5575673-5819613	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	500 mb	anl	ens mean	:HGT:500 mb:anl:ens mean
35	36	6505417	6760483	6505417-6760483	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	700 mb	anl	ens mean	:HGT:700 mb:anl:ens mean
40	41	7642597	7908462	7642597-7908462	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	850 mb	anl	ens mean	:HGT:850 mb:anl:ens mean
46	47	9144964	9419015	9144964-9419015	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	925 mb	anl	ens mean	:HGT:925 mb:anl:ens mean
55	56	11237567	11523719	11237567-11523719	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	1000 mb	anl	ens mean	:HGT:1000 mb:anl:ens mean
57	58	11777787	11930300	11777787-11930300	2023-01-04 12:00:00	2023-01-04 12:00:00	HGT	surface	anl	ens mean	:HGT:surface:anl:ens mean

ds = H.xarray("PWAT")
ds

<xarray.Dataset>
Dimensions:                (latitude: 361, longitude: 720)
Coordinates:
    time                   datetime64[ns] 2023-01-04T12:00:00
    step                   timedelta64[ns] 00:00:00
    atmosphereSingleLayer  float64 0.0
  * latitude               (latitude) float64 90.0 89.5 89.0 ... -89.5 -90.0
  * longitude              (longitude) float64 0.0 0.5 1.0 ... 358.5 359.0 359.5
    valid_time             datetime64[ns] 2023-01-04T12:00:00
Data variables:
    pwat                   (latitude, longitude) float32 4.3 4.3 4.3 ... 1.0 1.0
    gribfile_projection    object None
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          2
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    model:                   gefs
    product:                 atmos.5
    description:             Global Ensemble Forecast System (GEFS)
    remote_grib:             https://noaa-gefs-pds.s3.amazonaws.com/gefs.2023...
    local_grib:              /home/blaylock/data/gefs/20230104/subset_f4ef1a3...
    search:            PWAT

xarray.Dataset

• Dimensions:
  • latitude: 361
  • longitude: 720
• Coordinates: (6)
  • time
    ()
    datetime64[ns]
    2023-01-04T12:00:00

    long_name :

    initial time of forecast

    standard_name :

    forecast_reference_time

    array('2023-01-04T12:00:00.000000000', dtype='datetime64[ns]')

  • step
    ()
    timedelta64[ns]
    00:00:00

    long_name :

    time since forecast_reference_time

    standard_name :

    forecast_period

    array(0, dtype='timedelta64[ns]')

  • atmosphereSingleLayer
    ()
    float64
    0.0

    long_name :

    original GRIB coordinate for key: level(atmosphereSingleLayer)

    units :

    1

    array(0.)

  • latitude
    (latitude)
    float64
    90.0 89.5 89.0 ... -89.5 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([ 90. ,  89.5,  89. , ..., -89. , -89.5, -90. ])

  • longitude
    (longitude)
    float64
    0.0 0.5 1.0 ... 358.5 359.0 359.5

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5])

  • valid_time
    ()
    datetime64[ns]
    2023-01-04T12:00:00

    standard_name :

    time

    long_name :

    time

    array('2023-01-04T12:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (2)
  • pwat
    (latitude, longitude)
    float32
    4.3 4.3 4.3 4.3 ... 1.0 1.0 1.0 1.0

    GRIB_paramId :

    3054

    GRIB_dataType :

    pf

    GRIB_numberOfPoints :

    259920

    GRIB_typeOfLevel :

    atmosphereSingleLayer

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_NV :

    0

    GRIB_Nx :

    720

    GRIB_Ny :

    361

    GRIB_cfName :

    unknown

    GRIB_cfVarName :

    pwat

    GRIB_gridDefinitionDescription :

    Latitude/longitude. Also called equidistant cylindrical, or Plate Carree

    GRIB_iDirectionIncrementInDegrees :

    0.5

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.5

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    90.0

    GRIB_latitudeOfLastGridPointInDegrees :

    -90.0

    GRIB_longitudeOfFirstGridPointInDegrees :

    0.0

    GRIB_longitudeOfLastGridPointInDegrees :

    359.5

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    Precipitable water

    GRIB_parameterName :

    Precipitable water

    GRIB_parameterUnits :

    kg m-2

    GRIB_shortName :

    pwat

    GRIB_stepRange :

    0

    GRIB_totalNumber :

    30

    GRIB_units :

    kg m**-2

    long_name :

    Precipitable water

    units :

    kg m**-2

    standard_name :

    unknown

    grid_mapping :

    gribfile_projection

    array([[4.3, 4.3, 4.3, ..., 4.3, 4.3, 4.3],
           [4.7, 4.7, 4.7, ..., 4.7, 4.7, 4.7],
           [4.8, 4.8, 4.8, ..., 4.8, 4.8, 4.8],
           ...,
           [1.2, 1.2, 1.2, ..., 1.2, 1.2, 1.2],
           [1.1, 1.1, 1.1, ..., 1.1, 1.1, 1.1],
           [1. , 1. , 1. , ..., 1. , 1. , 1. ]], dtype=float32)

  • gribfile_projection
    ()
    object
    None

    crs_wkt :

    GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]

    semi_major_axis :

    6371229.0

    semi_minor_axis :

    6371229.0

    inverse_flattening :

    0.0

    reference_ellipsoid_name :

    unknown

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    unknown

    grid_mapping_name :

    latitude_longitude

    long_name :

    GEFS model grid projection

    array(None, dtype=object)

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 90.0,  89.5,  89.0,  88.5,  88.0,  87.5,  87.0,  86.5,  86.0,
                   85.5,
                  ...
                  -85.5, -86.0, -86.5, -87.0, -87.5, -88.0, -88.5, -89.0, -89.5,
                  -90.0],
                 dtype='float64', name='latitude', length=361))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([  0.0,   0.5,   1.0,   1.5,   2.0,   2.5,   3.0,   3.5,   4.0,
                    4.5,
                  ...
                  355.0, 355.5, 356.0, 356.5, 357.0, 357.5, 358.0, 358.5, 359.0,
                  359.5],
                 dtype='float64', name='longitude', length=720))

• Attributes: (12)

  GRIB_edition :

  2

  GRIB_centre :

  kwbc

  GRIB_centreDescription :

  US National Weather Service - NCEP

  GRIB_subCentre :

  2

  Conventions :

  CF-1.7

  institution :

  US National Weather Service - NCEP

  model :

  gefs

  product :

  atmos.5

  description :

  Global Ensemble Forecast System (GEFS)

  remote_grib :

  https://noaa-gefs-pds.s3.amazonaws.com/gefs.20230104/12/atmos/pgrb2ap5/geavg.t12z.pgrb2a.0p50.f000

  local_grib :

  /home/blaylock/data/gefs/20230104/subset_f4ef1a3b__geavg.t12z.pgrb2a.0p50.f000

  search :

  PWAT

ds2 = H.xarray("HGT:500 mb")
ds2

<xarray.Dataset>
Dimensions:              (latitude: 361, longitude: 720)
Coordinates:
    time                 datetime64[ns] 2023-01-04T12:00:00
    step                 timedelta64[ns] 00:00:00
    isobaricInhPa        float64 500.0
  * latitude             (latitude) float64 90.0 89.5 89.0 ... -89.0 -89.5 -90.0
  * longitude            (longitude) float64 0.0 0.5 1.0 ... 358.5 359.0 359.5
    valid_time           datetime64[ns] 2023-01-04T12:00:00
Data variables:
    gh                   (latitude, longitude) float32 5.248e+03 ... 5.023e+03
    gribfile_projection  object None
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          2
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    model:                   gefs
    product:                 atmos.5
    description:             Global Ensemble Forecast System (GEFS)
    remote_grib:             https://noaa-gefs-pds.s3.amazonaws.com/gefs.2023...
    local_grib:              /home/blaylock/data/gefs/20230104/subset_f4efbd3...
    search:            HGT:500 mb

xarray.Dataset

• Dimensions:
  • latitude: 361
  • longitude: 720
• Coordinates: (6)
  • time
    ()
    datetime64[ns]
    2023-01-04T12:00:00

    long_name :

    initial time of forecast

    standard_name :

    forecast_reference_time

    array('2023-01-04T12:00:00.000000000', dtype='datetime64[ns]')

  • step
    ()
    timedelta64[ns]
    00:00:00

    long_name :

    time since forecast_reference_time

    standard_name :

    forecast_period

    array(0, dtype='timedelta64[ns]')

  • isobaricInhPa
    ()
    float64
    500.0

    long_name :

    pressure

    units :

    hPa

    positive :

    down

    stored_direction :

    decreasing

    standard_name :

    air_pressure

    array(500.)

  • latitude
    (latitude)
    float64
    90.0 89.5 89.0 ... -89.5 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([ 90. ,  89.5,  89. , ..., -89. , -89.5, -90. ])

  • longitude
    (longitude)
    float64
    0.0 0.5 1.0 ... 358.5 359.0 359.5

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5])

  • valid_time
    ()
    datetime64[ns]
    2023-01-04T12:00:00

    standard_name :

    time

    long_name :

    time

    array('2023-01-04T12:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (2)
  • gh
    (latitude, longitude)
    float32
    5.248e+03 5.248e+03 ... 5.023e+03

    GRIB_paramId :

    156

    GRIB_dataType :

    pf

    GRIB_numberOfPoints :

    259920

    GRIB_typeOfLevel :

    isobaricInhPa

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_NV :

    0

    GRIB_Nx :

    720

    GRIB_Ny :

    361

    GRIB_cfName :

    geopotential_height

    GRIB_cfVarName :

    gh

    GRIB_gridDefinitionDescription :

    Latitude/longitude. Also called equidistant cylindrical, or Plate Carree

    GRIB_iDirectionIncrementInDegrees :

    0.5

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.5

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    90.0

    GRIB_latitudeOfLastGridPointInDegrees :

    -90.0

    GRIB_longitudeOfFirstGridPointInDegrees :

    0.0

    GRIB_longitudeOfLastGridPointInDegrees :

    359.5

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    Geopotential height

    GRIB_parameterName :

    Geopotential height

    GRIB_parameterUnits :

    gpm

    GRIB_shortName :

    gh

    GRIB_stepRange :

    0

    GRIB_totalNumber :

    30

    GRIB_units :

    gpm

    long_name :

    Geopotential height

    units :

    gpm

    standard_name :

    geopotential_height

    grid_mapping :

    gribfile_projection

    array([[5248.3467, 5248.3467, 5248.3467, ..., 5248.3467, 5248.3467,
            5248.3467],
           [5248.8066, 5248.8467, 5248.8867, ..., 5248.687 , 5248.727 ,
            5248.767 ],
           [5250.187 , 5250.247 , 5250.3066, ..., 5249.967 , 5250.027 ,
            5250.107 ],
           ...,
           [5028.507 , 5028.727 , 5028.967 , ..., 5027.8267, 5028.067 ,
            5028.287 ],
           [5025.707 , 5025.8267, 5025.927 , ..., 5025.3867, 5025.487 ,
            5025.607 ],
           [5022.787 , 5022.787 , 5022.787 , ..., 5022.787 , 5022.787 ,
            5022.787 ]], dtype=float32)

  • gribfile_projection
    ()
    object
    None

    crs_wkt :

    GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]

    semi_major_axis :

    6371229.0

    semi_minor_axis :

    6371229.0

    inverse_flattening :

    0.0

    reference_ellipsoid_name :

    unknown

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    unknown

    grid_mapping_name :

    latitude_longitude

    long_name :

    GEFS model grid projection

    array(None, dtype=object)

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 90.0,  89.5,  89.0,  88.5,  88.0,  87.5,  87.0,  86.5,  86.0,
                   85.5,
                  ...
                  -85.5, -86.0, -86.5, -87.0, -87.5, -88.0, -88.5, -89.0, -89.5,
                  -90.0],
                 dtype='float64', name='latitude', length=361))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([  0.0,   0.5,   1.0,   1.5,   2.0,   2.5,   3.0,   3.5,   4.0,
                    4.5,
                  ...
                  355.0, 355.5, 356.0, 356.5, 357.0, 357.5, 358.0, 358.5, 359.0,
                  359.5],
                 dtype='float64', name='longitude', length=720))

• Attributes: (12)

  GRIB_edition :

  2

  GRIB_centre :

  kwbc

  GRIB_centreDescription :

  US National Weather Service - NCEP

  GRIB_subCentre :

  2

  Conventions :

  CF-1.7

  institution :

  US National Weather Service - NCEP

  model :

  gefs

  product :

  atmos.5

  description :

  Global Ensemble Forecast System (GEFS)

  remote_grib :

  https://noaa-gefs-pds.s3.amazonaws.com/gefs.20230104/12/atmos/pgrb2ap5/geavg.t12z.pgrb2a.0p50.f000

  local_grib :

  /home/blaylock/data/gefs/20230104/subset_f4efbd36__geavg.t12z.pgrb2a.0p50.f000

  search :

  HGT:500 mb

ax = EasyMap("50m", crs=ds.herbie.crs, figsize=[10, 10]).STATES().BORDERS().ax
p = ax.pcolormesh(
    ds.longitude, ds.latitude, ds.pwat, transform=pc, cmap="Greens", vmin=0
)
plt.colorbar(
    p,
    ax=ax,
    orientation="horizontal",
    pad=0.01,
    shrink=0.8,
    label=f"{ds.pwat.GRIB_name}({ds.pwat.GRIB_units})",
)

# Add 500 hPa Geopotential height
ax.contour(
    ds2.longitude,
    ds2.latitude,
    ds2.gh,
    colors="k",
    linewidths=0.5,
    levels=range(0, 10000, 120),
)

ax.set_title(
    f"{ds.model.upper()} ensemble mean\nValid: {ds.valid_time.dt.strftime('%H:%M UTC %d %b %Y').item()}",
    loc="left",
)
ax.set_title(ds.pwat.GRIB_name, loc="right")

Text(1.0, 1.0, 'Precipitable water')

# Get a specific member
H5 = Herbie(
    "2022-01-01",
    model="gefs",
    product="atmos.5",
    member=5,
)
H5.xarray("TMP:2 m")

## If a user wants all the members in
## a single Dataframe, I'll let the user
## concat it themselves.

✅ Found ┊ model=gefs ┊ product=atmos.5 ┊ 2022-Jan-01 00:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws

<xarray.Dataset>
Dimensions:              (latitude: 361, longitude: 720)
Coordinates:
    number               int64 5
    time                 datetime64[ns] 2022-01-01
    step                 timedelta64[ns] 00:00:00
    heightAboveGround    float64 2.0
  * latitude             (latitude) float64 90.0 89.5 89.0 ... -89.0 -89.5 -90.0
  * longitude            (longitude) float64 0.0 0.5 1.0 ... 358.5 359.0 359.5
    valid_time           datetime64[ns] 2022-01-01
Data variables:
    t2m                  (latitude, longitude) float32 244.5 244.5 ... 252.3
    gribfile_projection  object None
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          2
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    model:                   gefs
    product:                 atmos.5
    description:             Global Ensemble Forecast System (GEFS)
    remote_grib:             https://noaa-gefs-pds.s3.amazonaws.com/gefs.2022...
    local_grib:              /home/blaylock/data/gefs/20220101/subset_04ef8ac...
    search:            TMP:2 m

xarray.Dataset

• Dimensions:
  • latitude: 361
  • longitude: 720
• Coordinates: (7)
  • number
    ()
    int64
    5

    long_name :

    ensemble member numerical id

    units :

    1

    standard_name :

    realization

    array(5)

  • time
    ()
    datetime64[ns]
    2022-01-01

    long_name :

    initial time of forecast

    standard_name :

    forecast_reference_time

    array('2022-01-01T00:00:00.000000000', dtype='datetime64[ns]')

  • step
    ()
    timedelta64[ns]
    00:00:00

    long_name :

    time since forecast_reference_time

    standard_name :

    forecast_period

    array(0, dtype='timedelta64[ns]')

  • heightAboveGround
    ()
    float64
    2.0

    long_name :

    height above the surface

    units :

    m

    positive :

    up

    standard_name :

    height

    array(2.)

  • latitude
    (latitude)
    float64
    90.0 89.5 89.0 ... -89.5 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([ 90. ,  89.5,  89. , ..., -89. , -89.5, -90. ])

  • longitude
    (longitude)
    float64
    0.0 0.5 1.0 ... 358.5 359.0 359.5

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5])

  • valid_time
    ()
    datetime64[ns]
    2022-01-01

    standard_name :

    time

    long_name :

    time

    array('2022-01-01T00:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (2)
  • t2m
    (latitude, longitude)
    float32
    244.5 244.5 244.5 ... 252.3 252.3

    GRIB_paramId :

    167

    GRIB_dataType :

    pf

    GRIB_numberOfPoints :

    259920

    GRIB_typeOfLevel :

    heightAboveGround

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_NV :

    0

    GRIB_Nx :

    720

    GRIB_Ny :

    361

    GRIB_cfName :

    air_temperature

    GRIB_cfVarName :

    t2m

    GRIB_gridDefinitionDescription :

    Latitude/longitude. Also called equidistant cylindrical, or Plate Carree

    GRIB_iDirectionIncrementInDegrees :

    0.5

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.5

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    90.0

    GRIB_latitudeOfLastGridPointInDegrees :

    -90.0

    GRIB_longitudeOfFirstGridPointInDegrees :

    0.0

    GRIB_longitudeOfLastGridPointInDegrees :

    359.5

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    2 metre temperature

    GRIB_parameterName :

    Temperature

    GRIB_parameterUnits :

    K

    GRIB_shortName :

    2t

    GRIB_stepRange :

    0

    GRIB_totalNumber :

    30

    GRIB_units :

    K

    long_name :

    2 metre temperature

    units :

    K

    standard_name :

    air_temperature

    grid_mapping :

    gribfile_projection

    array([[244.54944, 244.54944, 244.54944, ..., 244.54944, 244.54944,
            244.54944],
           [244.62944, 244.62944, 244.62944, ..., 244.61943, 244.61943,
            244.62944],
           [245.09943, 245.11943, 245.12944, ..., 245.04944, 245.06943,
            245.07944],
           ...,
           [251.91943, 251.89943, 251.87944, ..., 251.97943, 251.95943,
            251.93944],
           [252.22943, 252.21944, 252.20943, ..., 252.25943, 252.24944,
            252.23944],
           [252.25943, 252.25943, 252.25943, ..., 252.25943, 252.25943,
            252.25943]], dtype=float32)

  • gribfile_projection
    ()
    object
    None

    crs_wkt :

    GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]

    semi_major_axis :

    6371229.0

    semi_minor_axis :

    6371229.0

    inverse_flattening :

    0.0

    reference_ellipsoid_name :

    unknown

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    unknown

    grid_mapping_name :

    latitude_longitude

    long_name :

    GEFS model grid projection

    array(None, dtype=object)

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 90.0,  89.5,  89.0,  88.5,  88.0,  87.5,  87.0,  86.5,  86.0,
                   85.5,
                  ...
                  -85.5, -86.0, -86.5, -87.0, -87.5, -88.0, -88.5, -89.0, -89.5,
                  -90.0],
                 dtype='float64', name='latitude', length=361))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([  0.0,   0.5,   1.0,   1.5,   2.0,   2.5,   3.0,   3.5,   4.0,
                    4.5,
                  ...
                  355.0, 355.5, 356.0, 356.5, 357.0, 357.5, 358.0, 358.5, 359.0,
                  359.5],
                 dtype='float64', name='longitude', length=720))

• Attributes: (12)

  GRIB_edition :

  2

  GRIB_centre :

  kwbc

  GRIB_centreDescription :

  US National Weather Service - NCEP

  GRIB_subCentre :

  2

  Conventions :

  CF-1.7

  institution :

  US National Weather Service - NCEP

  model :

  gefs

  product :

  atmos.5

  description :

  Global Ensemble Forecast System (GEFS)

  remote_grib :

  https://noaa-gefs-pds.s3.amazonaws.com/gefs.20220101/00/atmos/pgrb2ap5/gep05.t00z.pgrb2a.0p50.f000

  local_grib :

  /home/blaylock/data/gefs/20220101/subset_04ef8acd__gep05.t00z.pgrb2a.0p50.f000

  search :

  TMP:2 m

GEFS-Reanalysis Data

This section demonstrates using the Global Ensemble Forecast System (GEFS) reanalysis: 2000-2019 dataset.

The GEFS version 12 reanalysis is available on Amazon Web Services and can be retrieved with Herbie.

Note: The GEFS directory structure is different than other models that Herbie can access, and that makes Herbie’s access to these files little awkward. Instead of grouping the GRIB fields by forecast hour where there are many different variables for the same lead time in the same file, the GEFS files are grouped into the same variable per file with each GRIB message being a different lead time. This changes the way a user would use Herbie to access GEFS data–a user will need to supply a “variable_level” argument to access a full file. For subsetting by specific grib messages, you will use the “search” argument to key in on the message of interest. You will still need to give a value for “fxx” to tell Herbie which directory to look for.

Yeah, it’s a little different paradigm for Herbie, but we can work with it. It may be nice to write a herbie.tool to help make calls to Herbie a little more simple (like a custom bulk_download script)

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

from herbie import Herbie
from paint.standard2 import cm_tmp, cm_wind, cm_wave_height
from toolbox import EasyMap, pc

Retrieve a full file

Download a full GRIB2 file to your local system.

Remember to specify the following: - fxx is the lead time lead time, a number between 0 and 384. If you are getting the full file, this only tells Herbie what folder to look in (Days:0-10 or Days:10-16). - member is the ensemble member. 0 is the control and a value between 1, 2, 3, or 4 is a perturbation member. - variable_level is the name of the file to obtain.

H = Herbie(
    "2017-03-14", model="gefs_reforecast", fxx=12, member=0, variable_level="tmp_2m"
)
H.download(verbose=True)

✅ Found ┊ model=gefs_reforecast ┊ product=GEFSv12/reforecast ┊ 2017-Mar-14 00:00 UTC F12 ┊ GRIB2 @ aws ┊ IDX @ aws
👨🏻‍🏭 Created directory: [/home/blaylock/data/gefs_reforecast/20170314]
✅ Success! Downloaded GEFS_REFORECAST from aws                 
        src: /home/blaylock/data/gefs_reforecast/20170314/tmp_2m_2017031400_c00.grib2
        dst: /home/blaylock/data/gefs_reforecast/20170314/tmp_2m_2017031400_c00.grib2

PosixPath('/home/blaylock/data/gefs_reforecast/20170314/tmp_2m_2017031400_c00.grib2')

Download/Retrieve a subset

Subsetting uses the search argument to parse out information from the GRIB2’s index file. Since the variables of each message in a file are all the same, we need to set the search to key in on the lead time we are interested in.

Look at the index file to see how to key in on a specific GRIB message.

# Look at the search_this column of the index DataFrame
H.inventory().search_this

grib_message
1.0     :TMP:2 m above ground:3 hour fcst:ENS=low-res ctl
2.0     :TMP:2 m above ground:6 hour fcst:ENS=low-res ctl
3.0     :TMP:2 m above ground:9 hour fcst:ENS=low-res ctl
4.0     :TMP:2 m above ground:12 hour fcst:ENS=low-res...
5.0     :TMP:2 m above ground:15 hour fcst:ENS=low-res...
                              ...
76.0    :TMP:2 m above ground:228 hour fcst:ENS=low-re...
77.0    :TMP:2 m above ground:231 hour fcst:ENS=low-re...
78.0    :TMP:2 m above ground:234 hour fcst:ENS=low-re...
79.0    :TMP:2 m above ground:237 hour fcst:ENS=low-re...
80.0    :TMP:2 m above ground:240 hour fcst:ENS=low-re...
Name: search_this, Length: 80, dtype: object

# Get the 15-h forecast
ds = H.xarray(":15 hour fcst:")

ds

<xarray.Dataset>
Dimensions:              (latitude: 721, longitude: 1440)
Coordinates:
    number               int32 0
    time                 datetime64[ns] 2017-03-14
    step                 timedelta64[ns] 15:00:00
    heightAboveGround    float64 2.0
  * latitude             (latitude) float64 90.0 89.75 89.5 ... -89.75 -90.0
  * longitude            (longitude) float64 0.0 0.25 0.5 ... 359.2 359.5 359.8
    valid_time           datetime64[ns] 2017-03-14T15:00:00
Data variables:
    t2m                  (latitude, longitude) float32 ...
    gribfile_projection  object None
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          2
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    model:                   gefs
    product:                 GEFSv12/reforecast
    description:             Global Ensemble Forecast System (GEFS)
    remote_grib:             C:\Users\blayl_depgywe\data\gefs\20170314\tmp_2m...
    local_grib:              C:\Users\blayl_depgywe\data\gefs\20170314\tmp_2m...
    search:            :15 hour fcst:

xarray.Dataset

• Dimensions:
  • latitude: 721
  • longitude: 1440
• Coordinates: (7)
  • number
    ()
    int32
    0

    long_name :

    ensemble member numerical id

    units :

    1

    standard_name :

    realization

    array(0)

  • time
    ()
    datetime64[ns]
    2017-03-14

    long_name :

    initial time of forecast

    standard_name :

    forecast_reference_time

    array('2017-03-14T00:00:00.000000000', dtype='datetime64[ns]')

  • step
    ()
    timedelta64[ns]
    15:00:00

    long_name :

    time since forecast_reference_time

    standard_name :

    forecast_period

    array(54000000000000, dtype='timedelta64[ns]')

  • heightAboveGround
    ()
    float64
    2.0

    long_name :

    height above the surface

    units :

    m

    positive :

    up

    standard_name :

    height

    array(2.)

  • latitude
    (latitude)
    float64
    90.0 89.75 89.5 ... -89.75 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ])

  • longitude
    (longitude)
    float64
    0.0 0.25 0.5 ... 359.2 359.5 359.8

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([0.0000e+00, 2.5000e-01, 5.0000e-01, ..., 3.5925e+02, 3.5950e+02,
           3.5975e+02])

  • valid_time
    ()
    datetime64[ns]
    ...

    standard_name :

    time

    long_name :

    time

    array('2017-03-14T15:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (2)
  • t2m
    (latitude, longitude)
    float32
    ...

    GRIB_paramId :

    167

    GRIB_dataType :

    cf

    GRIB_numberOfPoints :

    1038240

    GRIB_typeOfLevel :

    heightAboveGround

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_NV :

    0

    GRIB_Nx :

    1440

    GRIB_Ny :

    721

    GRIB_cfName :

    air_temperature

    GRIB_cfVarName :

    t2m

    GRIB_gridDefinitionDescription :

    Latitude/longitude. Also called equidistant cylindrical, or Plate Carree

    GRIB_iDirectionIncrementInDegrees :

    0.25

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.25

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    90.0

    GRIB_latitudeOfLastGridPointInDegrees :

    -90.0

    GRIB_longitudeOfFirstGridPointInDegrees :

    0.0

    GRIB_longitudeOfLastGridPointInDegrees :

    359.75

    GRIB_missingValue :

    9999

    GRIB_name :

    2 metre temperature

    GRIB_parameterName :

    Temperature

    GRIB_parameterUnits :

    K

    GRIB_shortName :

    2t

    GRIB_stepRange :

    15

    GRIB_totalNumber :

    10

    GRIB_units :

    K

    long_name :

    2 metre temperature

    units :

    K

    standard_name :

    air_temperature

    grid_mapping :

    gribfile_projection

    [1038240 values with dtype=float32]

  • gribfile_projection
    ()
    object
    None

    crs_wkt :

    GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]

    semi_major_axis :

    6371229.0

    semi_minor_axis :

    6371229.0

    inverse_flattening :

    0.0

    reference_ellipsoid_name :

    unknown

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    unknown

    grid_mapping_name :

    latitude_longitude

    long_name :

    GEFS model grid projection

    array(None, dtype=object)

• Attributes: (12)

  GRIB_edition :

  2

  GRIB_centre :

  kwbc

  GRIB_centreDescription :

  US National Weather Service - NCEP

  GRIB_subCentre :

  2

  Conventions :

  CF-1.7

  institution :

  US National Weather Service - NCEP

  model :

  gefs

  product :

  GEFSv12/reforecast

  description :

  Global Ensemble Forecast System (GEFS)

  remote_grib :

  C:\Users\blayl_depgywe\data\gefs\20170314\tmp_2m_2017031400_c00.grib2

  local_grib :

  C:\Users\blayl_depgywe\data\gefs\20170314\tmp_2m_2017031400_c00.grib2.subset_ac3478d69a3c81fa62e60f5c3696165a4e5e6ac4

  search :

  :15 hour fcst:

ax = EasyMap("50m", crs=ds.herbie.crs, figsize=[10, 10]).STATES().BORDERS().ax
p = ax.pcolormesh(
    ds.longitude, ds.latitude, ds.t2m, transform=pc, **cm_tmp(units="K").cmap_kwargs
)
plt.colorbar(
    p, ax=ax, orientation="horizontal", pad=0.05, **cm_tmp(units="K").cbar_kwargs
)

ax.set_title(
    f"{ds.model.upper()}: {H.product_description}\nValid: {ds.valid_time.dt.strftime('%H:%M UTC %d %b %Y').item()}",
    loc="left",
)
ax.set_title(ds.t2m.GRIB_name, loc="right")

C:\Users\blayl_depgywe\miniconda3\envs\herbie\lib\site-packages\metpy\xarray.py:353: UserWarning: More than one time coordinate present for variable "gribfile_projection".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'

Text(1.0, 1.0, '2 metre temperature')

What are valid values for variable_level?

You need to look at the file structure within the GEFS bucket to know what is avaialble. We can use s3fs to tell us what we can use for our variable_level argument.

import s3fs
import pandas as pd

# List files in the GEFS bucket for a day
fs = s3fs.S3FileSystem(anon=True)
files = fs.ls(
    path="noaa-gefs-retrospective/GEFSv12/reforecast/2015/2015010100/c00/Days:1-10"
)

# var_lev prefix
var_lev = [i.split("/")[-1].split("_") for i in files if i.endswith(".grib2")]

variable_levels_df = pd.DataFrame(var_lev, columns=["variable", "level", "a", "b", "c"])
variable_levels_df

	variable	level	a	b	c
0	acpcp	sfc	2015010100	c00.grib2	None
1	apcp	sfc	2015010100	c00.grib2	None
2	cape	sfc	2015010100	c00.grib2	None
3	cin	sfc	2015010100	c00.grib2	None
4	dlwrf	sfc	2015010100	c00.grib2	None
5	dswrf	sfc	2015010100	c00.grib2	None
6	gflux	sfc	2015010100	c00.grib2	None
7	gust	sfc	2015010100	c00.grib2	None
8	hgt	ceiling	2015010100	c00.grib2	None
9	hgt	hybr	2015010100	c00.grib2	None
10	hgt	pres	2015010100	c00.grib2	None
11	hgt	pres	abv700mb	2015010100	c00.grib2
12	hgt	sfc	2015010100	c00.grib2	None
13	hlcy	hgt	2015010100	c00.grib2	None
14	lhtfl	sfc	2015010100	c00.grib2	None
15	ncpcp	sfc	2015010100	c00.grib2	None
16	pbl	hgt	2015010100	c00.grib2	None
17	pres	hybr	2015010100	c00.grib2	None
18	pres	msl	2015010100	c00.grib2	None
19	pres	pvor	2015010100	c00.grib2	None
20	pres	sfc	2015010100	c00.grib2	None
21	pvort	isen	2015010100	c00.grib2	None
22	pwat	eatm	2015010100	c00.grib2	None
23	rh	hybr	2015010100	c00.grib2	None
24	sfcr	sfc	2015010100	c00.grib2	None
25	shtfl	sfc	2015010100	c00.grib2	None
26	soilw	bgrnd	2015010100	c00.grib2	None
27	spfh	2m	2015010100	c00.grib2	None
28	spfh	pres	2015010100	c00.grib2	None
29	spfh	pres	abv700mb	2015010100	c00.grib2
30	tcdc	eatm	2015010100	c00.grib2	None
31	tmax	2m	2015010100	c00.grib2	None
32	tmin	2m	2015010100	c00.grib2	None
33	tmp	2m	2015010100	c00.grib2	None
34	tmp	hybr	2015010100	c00.grib2	None
35	tmp	pres	2015010100	c00.grib2	None
36	tmp	pres	abv700mb	2015010100	c00.grib2
37	tmp	pvor	2015010100	c00.grib2	None
38	tmp	sfc	2015010100	c00.grib2	None
39	tozne	eatm	2015010100	c00.grib2	None
40	tsoil	bgrnd	2015010100	c00.grib2	None
41	uflx	sfc	2015010100	c00.grib2	None
42	ugrd	hgt	2015010100	c00.grib2	None
43	ugrd	hybr	2015010100	c00.grib2	None
44	ugrd	pres	2015010100	c00.grib2	None
45	ugrd	pres	abv700mb	2015010100	c00.grib2
46	ugrd	pvor	2015010100	c00.grib2	None
47	ulwrf	sfc	2015010100	c00.grib2	None
48	ulwrf	tatm	2015010100	c00.grib2	None
49	uswrf	sfc	2015010100	c00.grib2	None
50	vflx	sfc	2015010100	c00.grib2	None
51	vgrd	hgt	2015010100	c00.grib2	None
52	vgrd	hybr	2015010100	c00.grib2	None
53	vgrd	pres	2015010100	c00.grib2	None
54	vgrd	pres	abv700mb	2015010100	c00.grib2
55	vgrd	pvor	2015010100	c00.grib2	None
56	vvel	pres	2015010100	c00.grib2	None
57	vvel	pres	abv700mb	2015010100	c00.grib2
58	watr	sfc	2015010100	c00.grib2	None
59	weasd	sfc	2015010100	c00.grib2	None

# These are the available variables
variable_levels_df.variable.unique()

array(['acpcp', 'apcp', 'cape', 'cin', 'dlwrf', 'dswrf', 'gflux', 'gust',
       'hgt', 'hlcy', 'lhtfl', 'ncpcp', 'pbl', 'pres', 'pvort', 'pwat',
       'rh', 'sfcr', 'shtfl', 'soilw', 'spfh', 'tcdc', 'tmax', 'tmin',
       'tmp', 'tozne', 'tsoil', 'uflx', 'ugrd', 'ulwrf', 'uswrf', 'vflx',
       'vgrd', 'vvel', 'watr', 'weasd'], dtype=object)

# These are the available levels
variable_levels_df.level.unique()

array(['sfc', 'ceiling', 'hybr', 'pres', 'hgt', 'msl', 'pvor', 'isen',
       'eatm', 'bgrnd', '2m', 'tatm'], dtype=object)