🪄 Merge Datasets#

Often users want to get multiple files across several files (across time or forecast hours). Here is an example of how to get those data and merge them into a single xarray Dataset.

[2]:

from herbie import Herbie, FastHerbie
import xarray as xr
import pandas as pd

from itertools import chain

[3]:

# Create multiple Herbie objects for a range of dates
# Some data we want is in the RAP pressure grid file while other data is
# in the RAP native grid file.

dates = pd.date_range("2024-01-01", periods=3, freq="1H")

FH_prs = FastHerbie(dates, model="rap", product="awp130pgrb", fxx=[0])
FH_nat = FastHerbie(dates, model="rap", product="awp130bgrb", fxx=[0])

[28]:

FH_prs.file_exists, FH_nat.file_exists

[28]:

([▌▌Herbie RAP model awp130pgrb product initialized 2024-Jan-01 00:00 UTC F00 ┊ source=local,
  ▌▌Herbie RAP model awp130pgrb product initialized 2024-Jan-01 01:00 UTC F00 ┊ source=aws,
  ▌▌Herbie RAP model awp130pgrb product initialized 2024-Jan-01 02:00 UTC F00 ┊ source=aws],
 [▌▌Herbie RAP model awp130bgrb product initialized 2024-Jan-01 00:00 UTC F00 ┊ source=aws,
  ▌▌Herbie RAP model awp130bgrb product initialized 2024-Jan-01 01:00 UTC F00 ┊ source=aws,
  ▌▌Herbie RAP model awp130bgrb product initialized 2024-Jan-01 02:00 UTC F00 ┊ source=aws])

[4]:

# Get xarray data for pressure level data
ds_prs = [
    H.xarray("(?:TMP:2 m|GRD:10 m|DPT:2 m|GUST|TMP:1000 mb|TMP:500 mb)")
    for H in FH_prs.file_exists
]

# flatten the list of lists into just a list of Datasets
ds_prs = list(chain(*ds_prs))
len(ds_prs)

[4]:

[5]:

# Get xarray data for native level data
ds_nat = [H.xarray("(?:SOIL|VGTYP|TOSIL)") for H in FH_nat.file_exists]

# flatten the list of lists into just a list of Datasets
ds_nat = list(chain(*ds_nat))
len(ds_nat)

[5]:

[15]:

def merge_datasets(ds_list):
    """Merge list of Datasets together.

    Since cfgrib doesn't merge data in different "hypercubes", we will
    do the merge ourselves.

    Parameters
    ----------
    ds_list : list
        A list of xarray.Datasets, usually from the list of datasets
        returned by cfgrib when data is on multiple levels.
    """
    these = []
    for ds in ds_list:
        ds = ds.drop_vars("gribfile_projection")
        expand_dims = []
        for i in [
            "heightAboveGround",
            "time",
            "step",
            "isobaricInhPa",
            "depthBelowLandLayer",
        ]:
            if i in ds and i not in ds.dims:
                expand_dims.append(i)
        these.append(ds.expand_dims(expand_dims))
    return xr.merge(these, compat="override")

[16]:

# Merge prs and nat data into one dataframe
# (NOTE: I'm not 100% convinced I did the merge correct. Why are there
# so many NaN values?)
ds = merge_datasets(ds_prs + ds_nat)
ds

[21]:

# Get data at a single point
point = ds.sel(x=100, y=100).squeeze()
point

# (Note: to get data nearest a specific lat/lon point, you would have to
# do a nearest neighbor search on the latitude/longitude coordinates
# get the index, and then query the specific x/y location. Maybe a Ball
# Tree algorithm could be used for that.)

[29]:

point.u10