Nithiya Streethran

Nithiya Streethran

Environmental research and random stuff

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:

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
MÉRA GRIB data read using Xarray.
MÉRA GRIB data read using Xarray.

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()
Plot of MÉRA GRIB data read using Xarray. Plot of MÉRA GRIB data with lon/lat coordinates read using Xarray.
Plots of MÉRA GRIB data read using Xarray.

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
MÉRA NetCDF data (converted from GRIB using CDO) read using Xarray.
MÉRA NetCDF data (converted from GRIB using CDO) read using Xarray.

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
Lambert Conformal Conic projection for the MÉRA data.
Lambert Conformal Conic projection for the MÉRA data.

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()
Plot of MÉRA NetCDF data (converted from GRIB using CDO) read using Xarray. Plot of MÉRA NetCDF data (converted from GRIB using CDO) with lon/lat coordinates read using Xarray.
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_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()
Plot of MÉRA NetCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.
Plot of MÉRA NetCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.

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()
Facet plots of MÉRA NetCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.
Facet plots of MÉRA NetCDF data (converted from GRIB using CDO) clipped to the boundary of Ireland read using Xarray.

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()
MÉRA time series for the Cork Airport met station.
MÉRA time series for the Cork Airport met station.

Leave a comment

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

Markdown is supported.