GRIB: getting latitudes, longitudes and values (reduced Gaussian grid)

In this example we will work with GRIB data on a reduced Gaussian GRIB.

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_rgg_O32.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	reduced_gg
1	ecmf	t	isobaricInhPa	850	20180801	1200	0	an	0	reduced_gg
2	ecmf	t	isobaricInhPa	700	20180801	1200	0	an	0	reduced_gg
3	ecmf	t	isobaricInhPa	500	20180801	1200	0	an	0	reduced_gg
4	ecmf	t	isobaricInhPa	400	20180801	1200	0	an	0	reduced_gg
5	ecmf	t	isobaricInhPa	300	20180801	1200	0	an	0	reduced_gg

In the data each field is defined on an O32 octahedral reduced Gaussian grid. These grids cannot be represented as a 2D array so the field shape is always 1 dimensional:

ds[0].shape

(5248,)

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.

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

(numpy.ndarray, (3, 5248))

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 (5248,)).

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

(5248,)

87.86379883923263

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

(5248,)

0.0

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

(5248,)

271.55796813964844

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 (5248,)).

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

(numpy.ndarray, (8, 5248))

Latitudes can be accessed as llv[0].

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

(5248,)

87.86379883923263

Longitudes can be accessed as llv[1].

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

(5248,)

0.0

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

print(llv[2:].shape)

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

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

(6, 5248)
271.55796813964844
[271.55796814 277.64489746 272.19042969 255.98562622 243.70146179
 228.54637146]

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=(5248,)
lon shape=(5248,)

The field values can be accessed with to_numpy().

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

(5248,)

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)

(5248,)
(5248,)
(5248,)

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=(5248,)
lon shape=(5248,)

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, 5248)

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

v[:,0]

array([271.55796814, 277.64489746, 272.19042969, 255.98562622,
       243.70146179, 228.54637146])

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]

array([271.55798, 277.6449 , 272.19043, 255.98563, 243.70146, 228.54637],
      dtype=float32)

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

llv[:,0]

array([ 87.8638 ,   0.     , 271.55798, 277.6449 , 272.19043, 255.98563,
       243.70146, 228.54637], dtype=float32)