3
\$\begingroup\$

If you would like to download some grib data for yourself, the script below will pull girb data from the google api over a date range at an hourly interval.

Extract some hrrr grib2 data

from pathlib import Path
from typing import Iterator
from datetime import datetime
from shutil import copyfileobj
import warnings

from requests import Session, HTTPError
import pandas as pd


def google_api_hrrr_grib2_data(start: datetime, end: datetime) -> Iterator[str]:
    """
    url generator function for googleapis high-resolution-rapid-refresh dataset
    """
    base_url = "https://storage.googleapis.com/high-resolution-rapid-refresh/"
    date_range = pd.date_range(start, end, freq="h")
    yield from base_url + date_range.strftime("hrrr.%Y%m%d/conus/hrrr.t%Hz.wrfnatf00.grib2")

def download_hrrr(start: datetime, end: datetime, path: Path):
    """
    iterate over urls and save files to a Path directory
    """
    # request context manager
    with Session() as session:
        # iteratate over the generator function
        for url in google_api_hrrr_grib2_data(start=start, end=end):
            # add the filename to the path object
            save_to = path / ".".join(url.replace("hrrr.", "").split("/")[-3:])

            try:
                # make a http get request to the url
                res = session.get(url, stream=True)
                # on non 200 status code raise HTTPError
                res.raise_for_status()
                # save the file to the directory
                with save_to.open("wb") as fileout:
                    copyfileobj(res.raw, fileout)
                print("grib2 file saved at ", save_to)

            except HTTPError:
                warnings.warn(f"Warning: failed to download {url}")



data_hrrr = Path("data/hrrr/")
if __name__ == "__main__":
    download_hrrr(start=datetime(2022, 6, 15), end=datetime(2022, 6, 16), path=data_hrrr)

Transform the data into python objects

This environment requires xarray, cfgrib, and ideally dask.

dask is not required but speeds up the process.

if not using dask the chunks argument to xr.open_dataset will need to be set to None.

I wanted to try and make use of pythons decorator function. I have used it from several other libraries but never really implemented one to this extent myself.

The standard dimensions for every parameter are latitude, longitude, & valid_time

griblib/common.py

from typing import Callable, Literal, overload, Iterator
import xarray as xr

class GribBase:
    def __init__(self, files: list[str]) -> None:
        self._file_list = files

    def __repr__(self):
        return "{0}.propertys({1})".format(
            self.__class__.__name__,
            ", ".join(attr for attr in self.__dir__() if not attr.startswith("_")),
        )

    def iterfiles(self)->Iterator[str]:
        yield from self._file_list

def filter_by_level(level: str):
    """decorator"""

    lat_lon_vt = {"latitude", "longitude", "valid_time"}

    def generator(grib:GribBase, **kwargs):
        for file in grib.iterfiles():
            with xr.open_dataset(
                file,
                engine="cfgrib",
                **kwargs,
            ) as ds:
                if lat_lon_vt.issubset(ds.coords):
                    yield ds.drop_vars(coord for coord in ds.coords if coord not in lat_lon_vt)

    def func_wrapper(func: Callable[["GribBase"], dict[str, str]]):
        """the func() is the returned value from the function"""

        @overload
        def key_filter(
            grib: "GribBase",
            name: str = ...,
            stepType: Literal["max", "instant"] = ...,
            shortName: str = ...,
            standard_name: str = ...,
            **kwargs:str,
        ) -> xr.Dataset:
            ...

        def key_filter(grib: "GribBase", **kwargs: str) -> xr.Dataset:
            """the wrapped func"""

            default_return = func(grib)

            if not default_return:
                filter_by_keys = {"typeOfLevel": level} | kwargs
            else:
                filter_by_keys = {"typeOfLevel": level} | default_return | kwargs

            objs = generator(
                grib,
                filter_by_keys=filter_by_keys,
                chunks={},
            )
            return xr.concat(objs, dim="valid_time")

        return key_filter

    return func_wrapper

I went ahead an omitted the abc Docstring class

with the decorator function I can specify the level argument that I want to filter grib data and set an optional default return value. This is only necessary for certain variables.

griblib/hrrr.py

import xarray as xr
from typing import Literal
from griblib.common import GribBase, filter_by_level
## from griblib._abc import DocStrings



## class ByLevel(GribBase, DocStrings):
class ByLevel(GribBase):
    @filter_by_level("atmosphere")
    def atmosphere(self):
        """returns all attributes at with the shared `atmosphere` level"""

    @filter_by_level("isothermal")
    def isothermal(self):
        """returns all attributes at with the shared `isothermal` level"""
        return {"stepType": "instant"}

    @filter_by_level("hybrid")
    def hybrid(self):
        """returns all attributes at with the shared `hybrid` level"""

    @filter_by_level("heightAboveGroundLayer")
    def height_above_ground_layer(self):
        """returns all attributes at with the shared `heightAboveGroundLayer` level"""


class HRRR(ByLevel):
    def geopotential_height(self, level="isothermal") -> xr.Dataset:
        match level:
            case "isothermal":
                return super().isothermal(name="Geopotential Height", stepType="instant")
            case "hybrid":
                return super().hybrid(name="Geopotential Height")

    def pressure(self) -> xr.Dataset:
        return super().hybrid(name="Pressure")

    def unknown(self) -> xr.Dataset:
        return super().hybrid(name="unknown")

    def mixing_ratio(self, kind: Literal["rain", "snow", "cloud"] = "rain") -> xr.Dataset:
        return super().hybrid(name=f"{kind.title()} mixing ratio")

    def graupel(self) -> xr.Dataset:
        return super().hybrid(name="Graupel (snow pellets)")

    def particulate_matter(self, kind: Literal["fine", "coarse"] = "fine") -> xr.Dataset:
        return super().hybrid(name=f"Particulate matter ({kind})")

    def fraction_of_cloud_cover(self) -> xr.Dataset:
        return super().hybrid(name="Fraction of cloud cover")

    def temperature(self) -> xr.Dataset:
        return super().hybrid(name="Temperature")

    def specific_humidity(self) -> xr.Dataset:
        return super().hybrid(name="Specific humidity")

    def u_component_of_wind(self) -> xr.Dataset:
        return super().hybrid(name="U component of wind")

    def v_component_of_wind(self):
        return super().hybrid(name="V component of wind")

    def vertical_velocity(self) -> xr.Dataset:
        return super().hybrid(name="Vertical velocity")

    def turbulent_kinetic_energy(self) -> xr.Dataset:
        return super().hybrid(name="Turbulent kinetic energy")

    def derived_radar_reflectivity(self, step_type: Literal["max", "instant"] = "instant") -> xr.Dataset:
        return super().isothermal(name="Derived radar reflectivity", stepType=step_type)

main.py

from glob import glob

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

from griblib.hrrr import HRRR


def open_gribs(files: list[str]) -> HRRR:
    return HRRR(files if isinstance(files, list) else [files])


def scatter_plot(ds: xr.Dataset, parameter: str, **kwargs):
    columns = {
        v.attrs["GRIB_shortName"]: v.attrs["long_name"].lower().replace(" ", "_") for v in ds.data_vars.values()
    }
    for vt, df in ds.to_dataframe().rename(columns=columns).groupby("valid_time"):
        plt.rcParams["axes.grid"] = True
        ax = plt.subplot(2, 1, 1, projection=ccrs.PlateCarree())
        plt.figure(vt.value)

        ax.coastlines("50m")

        ax.set_extent(
            [
                *df["longitude"].agg(["min", "max"]),
                *df["latitude"].agg(["min", "max"]),
            ],
            ccrs.PlateCarree(),
        )

        df.plot(
            kind="scatter",
            y="latitude",
            x="longitude",
            c=parameter,
            ax=ax,
            alpha=np.where(df[parameter] > 0, 0.2, 0),
            cmap=plt.get_cmap("jet"),
            s=0.25,
            figsize=(16, 12),
            **kwargs,
        )


if __name__ == "__main__":
    hrrr = open_gribs(glob("data/hrrr/*.grib2")[:2])
    hrrr.derived_radar_reflectivity().pipe(scatter_plot, parameter="derived_radar_reflectivity")

radar-image-one radar-image-two

\$\endgroup\$
0

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.