Handling Met Éireann Reanalysis (MÉRA) GRIB data

https://www.met.ie/climate/available-data/mera

Summary

Issues with MÉRA GRIB 1 files when read using software packages such as the Xarray Python library:

• The data uses 0/360 longitude coordinates instead of -180/180
  • data spans both negative and positive longitudes
• Projection information is present in the GRIB file but is not parsed when the data is read
• Lon/lat are multidimensional coordinates (y, x)
  • data with two-dimensional coordinates cannot be spatially selected (e.g. extracting data for a certain point, or clipping with a geometry)
  • x and y correspond to the index of the cells and are not coordinates in Lambert Conformal Conic projection

My solution:

• Use CDO to convert the GRIB 1 files to netCDF first
  • projection info is parsed and can be read by Xarray
  • one-dimensional coordinates in Lambert Conformal Conic projection
  • data can now be indexed or selected both spatially and temporally using Xarray
    • some metadata are lost (e.g. variable name and attributes) but these can be reassigned manually
    • this method combines all three time steps in the example 3-h forecast data (data needs to be split prior to conversion to avoid this)

Relevant links:

• https://docs.xarray.dev/en/stable/examples/multidimensional-coords.html
• https://docs.xarray.dev/en/stable/examples/ERA5-GRIB-example.html
• https://scitools.org.uk/cartopy/docs/latest/reference/projections.html
• https://confluence.ecmwf.int/display/OIFS/How+to+convert+GRIB+to+netCDF

Requirements

Install the following in a Conda environment:

conda install --channel conda-forge python=3.10 cdo rioxarray cfgrib netcdf4 geopandas dask cartopy matplotlib nc-time-axis pooch

Procedure

Import libraries:

import os
from datetime import date, datetime, timezone
import cartopy.crs as ccrs
import geopandas as gpd
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import matplotlib.units as munits
import numpy as np
import pooch
import xarray as xr

Download sample data

Create a directory to store data:

DATA_DIR = os.path.join("data", "MERA", "sample")
os.makedirs(DATA_DIR, exist_ok=True)

Download sample GRIB data (2 m temperature; 3-h forecasts):

URL = "https://www.met.ie/downloads/MERA_PRODYEAR_2015_06_11_105_2_0_FC3hr.grb"
FILE_NAME = "MERA_PRODYEAR_2015_06_11_105_2_0_FC3hr"

if not os.path.isfile(os.path.join(DATA_DIR, FILE_NAME)):
    pooch.retrieve(
        url=URL,
        known_hash=None,
        fname=f"{FILE_NAME}.grb",
        path=DATA_DIR
    )

    with open(
        os.path.join(DATA_DIR, f"{FILE_NAME}.txt"), "w", encoding="utf-8"
    ) as outfile:
        outfile.write(
            f"Data downloaded on: {datetime.now(tz=timezone.utc)}\n"
            f"Download URL: {URL}"
        )

Path to the example data file:

BASE_FILE_NAME = os.path.join(DATA_DIR, FILE_NAME)

Read the original GRIB data

data = xr.open_dataset(
    f"{BASE_FILE_NAME}.grb",
    decode_coords="all", chunks="auto", engine="cfgrib"
)

View the GRIB data:

data

MÉRA GRIB data read using Xarray.

Check if there’s any coordinate reference system (CRS) information in the data (there should be no output):

data.rio.crs

Save variable attributes for later:

t_attrs = data["t"].attrs

Convert longitude format

Convert 0/360 deg to -180/180 deg longitudes, and reassign attributes:

long_attrs = data.longitude.attrs
data = data.assign_coords(longitude=(((data.longitude + 180) % 360) - 180))
data.longitude.attrs = long_attrs

Visualise the data:

plt.figure(figsize=(9, 7))
(data.isel(time=0, step=2)["t"] - 273.15).plot.contourf(
    robust=True, cmap="Spectral_r", levels=11,
    cbar_kwargs={"label": "Temperature [°C]"}
)
plt.tight_layout()
plt.xlabel(None)
plt.ylabel(None)
plt.title(f"time={data.isel(time=0, step=2)['t']['time'].values}")
plt.show()

# use lon/lat
plt.figure(figsize=(9, 7))
(data.isel(time=0, step=2)["t"] - 273.15).plot.contourf(
    robust=True, cmap="Spectral_r", x="longitude", y="latitude",
    cbar_kwargs={"label": "Temperature [°C]"}, levels=11
)
plt.xlabel(None)
plt.ylabel(None)
plt.title(f"time={data.isel(time=0, step=2)['t']['time'].values}")
plt.tight_layout()
plt.show()

Plots of MÉRA GRIB data read using Xarray.

Convert GRIB to netCDF using CDO

Keeping only the third forecast step of the FC3hr data:

os.system(
    f"cdo -s -f nc4c -copy -seltimestep,3/{len(data['time']) * 3}/3 "
    f"{BASE_FILE_NAME}.grb {BASE_FILE_NAME}.nc"
)

Read the netCDF file:

data = xr.open_dataset(
    f"{BASE_FILE_NAME}.nc", decode_coords="all", chunks="auto"
)

Reassign attributes and rename variables:

data["var11"].attrs = t_attrs
data = data.rename({"var11": "t"})

View the data:

data

MÉRA netCDF data (converted from GRIB using CDO) read using Xarray.

View the CRS:

data.rio.crs

CRS.from_wkt('PROJCS["undefined",GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6367470,0]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_1SP"],PARAMETER["latitude_of_origin",53.5],PARAMETER["central_meridian",5],PARAMETER["scale_factor",1],PARAMETER["false_easting",1481641.67696368],PARAMETER["false_northing",537326.063885016],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

Plot using Lambert Conformal Conic projection

Define Lambert Conformal Conic projection for plots and transformations using metadata:

lambert_conformal = ccrs.LambertConformal(
    false_easting=data["Lambert_Conformal"].attrs["false_easting"],
    false_northing=data["Lambert_Conformal"].attrs["false_northing"],
    standard_parallels=[data["Lambert_Conformal"].attrs["standard_parallel"]],
    central_longitude=(
        data["Lambert_Conformal"].attrs["longitude_of_central_meridian"]
    ),
    central_latitude=(
        data["Lambert_Conformal"].attrs["latitude_of_projection_origin"]
    )
)

lambert_conformal

Lambert Conformal Conic projection for the MÉRA data.

<cartopy.crs.LambertConformal object at 0x7fc44026b100>

Plot data:

plt.figure(figsize=(9, 7))
(data.isel(time=0, height=0)["t"] - 273.15).plot.contourf(
    robust=True, cmap="Spectral_r", levels=11,
    cbar_kwargs={"label": "Temperature [°C]"}
)
plt.xlabel(None)
plt.ylabel(None)
plt.title(f"time={data.isel(time=0, height=0)['t']['time'].values}")
plt.tight_layout()
plt.show()

# in lon/lat
plt.figure(figsize=(9, 7))
ax = plt.axes(projection=lambert_conformal)
(data.isel(time=0, height=0)["t"] - 273.15).plot.contourf(
    ax=ax, robust=True, cmap="Spectral_r", x="x", y="y", levels=11,
    transform=lambert_conformal, cbar_kwargs={"label": "Temperature [°C]"}
)
ax.gridlines(
    draw_labels={"bottom": "x", "left": "y"},
    color="lightslategrey",
    linewidth=.5,
    x_inline=False,
    y_inline=False
)
ax.coastlines(resolution="10m", color="darkslategrey", linewidth=.75)
plt.title(f"time={data.isel(time=0, height=0)['t']['time'].values}")
plt.tight_layout()
plt.show()

Plots of MÉRA netCDF data (converted from GRIB using CDO) read using Xarray.

Clip to boundary of Ireland

Ireland boundary (derived from NUTS 2021):

ie = gpd.read_file(
  os.path.join("data", "boundary", "boundaries.gpkg"),
  layer="NUTS_RG_01M_2021_2157_IE"
)

Clip:

data_ie = data.rio.clip(
    ie.buffer(1).to_crs(lambert_conformal), all_touched=True
)

View data:

data_ie

MÉRA netCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.

Find number of grid cells with data:

len(
    data_ie.isel(time=0, height=0)["t"].values.flatten()[
        np.isfinite(data_ie.isel(time=0, height=0)["t"].values.flatten())
    ]
)

14490

Plot:

plt.figure(figsize=(9, 7))
ax = plt.axes(projection=ccrs.EuroPP())
(data_ie.isel(time=0, height=0)["t"] - 273.15).plot.contourf(
    ax=ax, robust=True, cmap="Spectral_r", x="x", y="y", levels=8,
    transform=lambert_conformal, cbar_kwargs={"label": "Temperature [°C]"}
)
ax.gridlines(
    draw_labels={"bottom": "x", "left": "y"},
    color="lightslategrey",
    linewidth=.5,
    x_inline=False,
    y_inline=False
)
ax.coastlines(resolution="10m", color="darkslategrey", linewidth=.75)
plt.title(
    "MERA_FC3hr, " +
    f"time={str(data_ie.isel(time=0, height=0)['t']['time'].values)[:19]}"
)
plt.tight_layout()
plt.show()

Plot of MÉRA netCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.

Extract time series for a grid point

Moorepark, Fermoy met station coordinates:

LON, LAT = -8.26389, 52.16389

Transform coordinates from lon/lat to Lambert Conformal Conic projection:

XLON, YLAT = lambert_conformal.transform_point(
    x=LON, y=LAT, src_crs=ccrs.PlateCarree()
)

XLON, YLAT

(579021.4574730808, 472854.55040691514)

Extract data for the nearest grid cell to the point:

data_ts = data_ie.sel({"x": XLON, "y": YLAT}, method="nearest")

View data:

data_ts

MÉRA netCDF data (converted from GRIB using CDO) for the Moorepark, Fermoy met station grid point read using Xarray.

Plot:

converter = mdates.ConciseDateConverter()
munits.registry[np.datetime64] = converter
munits.registry[date] = converter
munits.registry[datetime] = converter

plt.figure(figsize=(12, 4))
plt.plot(data_ts["time"], (data_ts["t"] - 273.15))
plt.ylabel("Temperature [°C]")
plt.title(f"MERA_FC3hr, lon={LON}, lat={LAT}")
plt.tight_layout()
plt.show()

MÉRA time series for the Moorepark, Fermoy met station.