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_in = ekd.from_source("sample", "tuv_pl_rgg_O32.grib").to_fieldlist()
ds = ds_in.sel({"parameter.variable": "t"})
ds.ls()

	parameter.variable	time.valid_datetime	time.base_datetime	time.step	vertical.level	vertical.level_type	ensemble.member	geography.grid_type
0	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	1000	pressure	0	reduced_gg
1	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	850	pressure	0	reduced_gg
2	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	700	pressure	0	reduced_gg
3	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	500	pressure	0	reduced_gg
4	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	400	pressure	0	reduced_gg
5	t	2018-08-01 12:00:00	2018-08-01 12:00:00	0 days	300	pressure	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() methods:

• earthkit.data.core.field.Field.data()

• earthkit.data.core.fieldlist.FieldList.data()

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

np.float64(87.86379883923263)

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

(5248,)

np.float64(0.0)

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

(5248,)

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

np.float64(87.86379883923263)

Longitudes can be accessed as llv[1].

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

(5248,)

np.float64(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 latlon() and to_numpy() methods

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

Fields

latlons() returns a dict:

lat, lon = ds[0].geography.latlons()
lat.shape, lon.shape

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

lat, lon = ds[0].geography.latlons(flatten=True)
v = ds[0].to_numpy(flatten=True)
lat.shape, lon.shape, v.shape

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

FieldLists

latlons() 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:

lat, lon = ds.geography.latlons()
lat.shape, lon.shape

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