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
Solution:
- Use CDO to convert the GRIB file 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
- https://github.com/corteva/rioxarray/issues/135
- https://docs.xarray.dev/en/stable/generated/xarray.DataArray.assign_coords.html
Requirements
Install the following in a Conda environment:
conda install --channel conda-forge python=3.10 cdo rioxarray geopandas dask cartopy matplotlib nc-time-axis pooch
Procedure
Download sample data
Import libraries:
import os
from datetime import date, datetime
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
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"
KNOWN_HASH = "3106063013265c1b1e9f535938aac7e391e2926b0df9ec15e2ed97e7fd0b565f"
pooch.retrieve(
url=URL,
known_hash=KNOWN_HASH,
fname="MERA_PRODYEAR_2015_06_11_105_2_0_FC3hr.grb",
path=DATA_DIR
)
Path to the example data file:
BASE_FILE_NAME = os.path.join(
DATA_DIR, "MERA_PRODYEAR_2015_06_11_105_2_0_FC3hr"
)
Read the original GRIB data
data = xr.open_dataset(
f"{BASE_FILE_NAME}.grb", decode_coords="all", chunks="auto"
)
View the GRIB data:
data

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 GRIB data:
plt.figure(figsize=(9, 7))
(data["t"][0][0] - 273.15).plot(
robust=True, levels=15, cmap="Spectral_r",
cbar_kwargs=dict(label="Temperature [°C]")
)
plt.tight_layout()
plt.show()
# use lon/lat
plt.figure(figsize=(9, 7))
(data["t"][0][0] - 273.15).plot(
robust=True, levels=15, cmap="Spectral_r", x="longitude", y="latitude",
cbar_kwargs=dict(label="Temperature [°C]")
)
plt.tight_layout()
plt.show()


Convert GRIB to NetCDF using CDO
os.system(f"cdo -f nc copy {BASE_FILE_NAME}.grb {BASE_FILE_NAME}.nc")
cdo copy: Processed 186250320 values from 1 variable over 720 timesteps [6.63s 328MB].
Read the NetCDF file:
data = xr.open_dataset(
f"{BASE_FILE_NAME}.nc", decode_coords="all", chunks="auto"
)
View the data:
data

View the coordinate reference system (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]]')
Save CRS info:
data_crs = data.rio.crs
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

Plot data:
plt.figure(figsize=(9, 7))
(data["var11"][0][0] - 273.15).plot(
robust=True, levels=15, cmap="Spectral_r",
cbar_kwargs=dict(label="Temperature [°C]")
)
plt.tight_layout()
plt.show()
# in lon/lat
plt.figure(figsize=(9, 7))
ax = plt.axes(projection=lambert_conformal)
(data["var11"][0][0] - 273.15).plot(
ax=ax, robust=True, levels=15, cmap="Spectral_r",
x="x", y="y", transform=lambert_conformal,
cbar_kwargs=dict(label="Temperature [°C]")
)
ax.gridlines(
draw_labels=dict(bottom="x", left="y"),
color="lightslategrey",
linewidth=.5,
x_inline=False,
y_inline=False
)
ax.coastlines(resolution="10m", color="darkslategrey", linewidth=.75)
plt.tight_layout()
plt.show()


Clip to boundary of Ireland
Ireland boundary (derived from NUTS 2021):
ie = gpd.read_file(
os.path.join("data", "boundary", "boundaries.gpkg"),
layer="NUTS_Ireland_ITM"
)
Clip:
data_ie = data.rio.clip(ie.buffer(2500).to_crs(data_crs))
Plot:
plt.figure(figsize=(9, 7))
ax = plt.axes(projection=lambert_conformal)
(data_ie["var11"][0][0] - 273.15).plot(
ax=ax, robust=True, levels=15, cmap="Spectral_r",
x="x", y="y", transform=lambert_conformal,
cbar_kwargs=dict(label="Temperature [°C]")
)
ax.gridlines(
draw_labels=dict(bottom="x", left="y"),
color="lightslategrey",
linewidth=.5,
x_inline=False,
y_inline=False
)
ax.coastlines(resolution="10m", color="darkslategrey", linewidth=.75)
plt.tight_layout()
plt.show()

Subset data at T=12:00 to create facet plots:
time_list = []
for time in data_ie["time"].values:
if "12:00" in str(time):
time_list.append(time)
data_sub = data_ie.sel(time=time_list)
fig = (data_sub["var11"] - 273.15).plot(
x="x", y="y", col="time", col_wrap=5, cmap="Spectral_r",
robust=True, cbar_kwargs=dict(aspect=40, label="Temperature [°C]"),
levels=15, transform=lambert_conformal,
subplot_kws=dict(projection=lambert_conformal)
)
for ax in fig.axes.flat:
ie.to_crs(lambert_conformal).boundary.plot(
ax=ax, color="darkslategrey", linewidth=.5
)
plt.show()

Extract time series for a point
Cork Airport met station coordinates:
LON = -8.48611
LAT = 51.84722
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
(557449.2236199714, 441074.0899205642)
Extract data for the nearest grid cell to the point:
data_ts = data.sel(dict(x=XLON, y=YLAT), method="nearest")
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["var11"] - 273.15))
plt.ylabel("Temperature [°C]")
plt.tight_layout()
plt.show()

Leave a comment
Your email address will not be published. Required fields are marked *.