GRIB: getting latitudes, longitudes and values (regular LL grid)

In this example we will work with GRIB data on a regular latitude-longitude grid.

import earthkit.data as ekd

We will only use the temperature fields, so we read the file and extract them with sel().

ds = ekd.from_source("sample", "tuv_pl.grib").sel(param="t")
ds.ls()

	centre	shortName	typeOfLevel	level	dataDate	dataTime	stepRange	dataType	number	gridType
0	ecmf	t	isobaricInhPa	1000	20180801	1200	0	an	0	regular_ll
1	ecmf	t	isobaricInhPa	850	20180801	1200	0	an	0	regular_ll
2	ecmf	t	isobaricInhPa	700	20180801	1200	0	an	0	regular_ll
3	ecmf	t	isobaricInhPa	500	20180801	1200	0	an	0	regular_ll
4	ecmf	t	isobaricInhPa	400	20180801	1200	0	an	0	regular_ll
5	ecmf	t	isobaricInhPa	300	20180801	1200	0	an	0	regular_ll

In the data each field is defined on a regular latitude-longitude grid. These grids can be represented as a 2D array so the field shape is always 2 dimensional:

ds[0].shape

(7, 12)

Using the data() method

The simplest way to access the latitudes, longitudes and values is using the data() method. It works both on fields and fieldlists.

Fields

data() returns the latitude, longitude and values arrays from a field as an ndarray. By default the field’s shape is kept.

llv = ds[0].data()
type(llv), llv.shape

(numpy.ndarray, (3, 7, 12))

Here llv[0] contains the latitudes, llv[1] the longitudes, while llv[2] the values. Each of these arrays has the same shape that of the field (7, 12)).

print(llv[0].shape)
# first latitude
llv[0, 0, 0]

(7, 12)

90.0

print(llv[1].shape)
# first longitude
llv[1, 0, 0]

(7, 12)

0.0

print(llv[2].shape)
# first value
llv[2, 0, 0]

(7, 12)

272.5641784667969

When using the flatten keyword 1D arrays are returned:

llv = ds[0].data(flatten=True)
llv.shape

(3, 84)

FieldLists

data() only works on a fieldlist if all the fields have the same grid. The first two elements of the resulting ndarray are the latitude and longitude arrays (shared between fields), while the rest of the elements are the value arrays per field. Since we have 6 fields in our data the size of the first axis of the resulting ndarray is 2+6=8. Each of these arrays has the same shape that of the field (7, 12)).

llv = ds.data()
type(llv), llv.shape

(numpy.ndarray, (8, 7, 12))

Latitudes can be accessed as llv[0].

print(llv[0].shape)
# first latitude
llv[0, 0, 0]

(7, 12)

90.0

Longitudes can be accessed as llv[1].

print(llv[1].shape)
# first longitude
llv[1, 0, 0]

(7, 12)

0.0

The field values can be accessed as llv[2:].

print(llv[2:].shape)

# first value in first field
print(llv[2, 0, 0])

# first value from all the fields
print(llv[2:, 0, 0])

(6, 7, 12)
272.5641784667969
[272.56417847 272.53916931 271.26531982 255.84306335 244.00323486
 226.65315247]

When using the flatten keyword 1D arrays are returned for latitude and longitude, and the values array will be flattened per field:

llv = ds.data(flatten=True)
llv.shape

(8, 84)

Using the to_latlon() and to_numpy() methods

We can also get the latitudes, longitudes and values by using the to_latlon() and to_to_numpy() methods.

Fields

to_latlon() returns a dict:

ll = ds[0].to_latlon()
for k, v in ll.items():
    print(f"{k} shape={v.shape}")

lat shape=(7, 12)
lon shape=(7, 12)

The field values can be accessed with to_numpy().

v = ds[0].to_numpy()
v.shape

(7, 12)

By default both methods keep the field’s shape, but we can use the flatten keyword to get 1D arrays:

ll = ds[0].to_latlon(flatten=True)
v = ds[0].to_numpy(flatten=True)
print(ll["lat"].shape)
print(ll["lon"].shape)
print(v.shape)

(84,)
(84,)
(84,)

FieldLists

to_latlon() only works on a fieldlist if all the fields have the same grid. In this case it returns the same dict that we would get for any of the fields:

ll = ds.to_latlon()
for k, v in ll.items():
    print(f"{k} shape={v.shape}")

lat shape=(7, 12)
lon shape=(7, 12)

to_numpy() only works on a fieldlist if all the fields has the same grid. It returns the array of the field value arrays.

v = ds.to_numpy()
v.shape

(6, 7, 12)

E.g. the first value from each field can be extracted as:

v[:,0,0]

array([272.56417847, 272.53916931, 271.26531982, 255.84306335,
       244.00323486, 226.65315247])

Specifying the array type

For all the methods above we can set the array type with the dtype keyword both for fields and fieldlists:

import numpy as np
v = ds.to_numpy(dtype=np.float32)
v[:,0,0]

array([272.56418, 272.53918, 271.26532, 255.84306, 244.00323, 226.65315],
      dtype=float32)

llv = ds.data(dtype=np.float32)

llv[:,0,0]

array([ 90.     ,   0.     , 272.56418, 272.53918, 271.26532, 255.84306,
       244.00323, 226.65315], dtype=float32)