Field arithmetic¶
earthkit-data supports standard arithmetic operators directly on Field and FieldList objects. Operations produce new objects and never modify their inputs.
The table below lists the supported operators.
Operator |
Example |
Description |
|---|---|---|
|
|
Addition |
|
|
Subtraction |
|
|
Multiplication |
|
|
True division |
|
|
Floor division |
|
|
Modulo |
|
|
Power |
|
|
Unary negation |
|
|
Comparison (returns boolean array field) |
[1]:
import numpy as np
import earthkit.data as ekd
Load data¶
Load a multi-field GRIB file so we have several fields to work with.
[2]:
ds = ekd.from_source("sample", "tuv_pl.grib").to_fieldlist()
len(ds)
[2]:
18
[3]:
# Extract individual fields
t_fields = ds.sel({"parameter.variable": "t"}) # temperature at all levels
u_fields = ds.sel({"parameter.variable": "u"}) # u-component of wind
v_fields = ds.sel({"parameter.variable": "v"}) # v-component of wind
t = t_fields[0] # single temperature field
u = u_fields[0] # single u-wind field
v = v_fields[0] # single v-wind field
print("Shape:", t.geography.shape())
print("Level:", t.vertical.level())
print("Param:", t.parameter.variable())
Shape: (7, 12)
Level: 1000
Param: t
Scalar arithmetic¶
Arithmetic with a scalar operates element-wise on every grid point. The result is a new Field whose metadata is inherited from the original.
[4]:
# Convert temperature from Kelvin to Celsius
t_celsius = t - 273.15
print(f"T(K) mean: {t.values.mean():.2f}")
print(f"T(°C) mean: {t_celsius.values.mean():.2f}")
T(K) mean: 279.71
T(°C) mean: 6.56
[5]:
# Scaling and offset
t_scaled = t * 0.5 + 100
print(f"scaled mean: {t_scaled.values.mean():.2f}")
# Field as a right hand side operand
t_scaled = 0.5 * t + 100
print(f"scaled mean: {t_scaled.values.mean():.2f}")
# Power
t_sq = t**2
print(f"squared mean: {t_sq.values.mean():.2f}")
# Unary negation
t_neg = -t
np.testing.assert_allclose(t_neg.values, -t.values)
print(f"negative mean: {t_neg.values.mean():.2f}")
scaled mean: 239.85
scaled mean: 239.85
squared mean: 78623.10
negative mean: -279.71
Field–field arithmetic¶
When both operands are Field objects with the same grid, operations are applied point-wise. A common example is computing wind speed from its components.
[6]:
# Wind speed: sqrt(u^2 + v^2)
wind_speed = (u**2 + v**2) ** 0.5
print(f"u mean: {u.values.mean():.2f} m/s")
print(f"v mean: {v.values.mean():.2f} m/s")
print(f"wind speed mean: {wind_speed.values.mean():.2f} m/s")
u mean: -0.38 m/s
v mean: -0.07 m/s
wind speed mean: 7.45 m/s
[7]:
# Difference between two fields (e.g. temperature anomaly relative to first level)
t_850 = t_fields.sel({"vertical.level": 850})[0]
t_500 = t_fields.sel({"vertical.level": 500})[0]
thickness = t_850 - t_500
print(f"T@850 mean: {t_850.values.mean():.2f} K")
print(f"T@500 mean: {t_500.values.mean():.2f} K")
print(f"diff mean: {thickness.values.mean():.2f} K")
T@850 mean: 272.73 K
T@500 mean: 252.89 K
diff mean: 19.84 K
FieldList arithmetic¶
Arithmetic between two FieldList objects applies the operation element-wise across corresponding fields (i.e. field 0 op field 0, field 1 op field 1, …). The two fieldlists must have the same length.
[8]:
# Convert all temperature fields from K to °C at once
t_celsius_all = t_fields - 273.15
print(f"Original length: {len(t_fields)}, result length: {len(t_celsius_all)}")
for i, (orig, conv) in enumerate(zip(t_fields, t_celsius_all)):
print(f" level {orig.vertical.level():4d}: {orig.values.mean():.2f} K -> {conv.values.mean():.2f} °C")
Original length: 6, result length: 6
level 1000: 279.71 K -> 6.56 °C
level 850: 272.73 K -> -0.42 °C
level 700: 265.17 K -> -7.98 °C
level 500: 252.89 K -> -20.26 °C
level 400: 241.53 K -> -31.62 °C
level 300: 227.84 K -> -45.31 °C
[9]:
# Element-wise wind speed across all levels
wind_speed_all = (u_fields**2 + v_fields**2) ** 0.5
print("Wind speed at each level (mean m/s):")
for f in wind_speed_all:
print(f" level {f.vertical.level():4d}: {f.values.mean():.2f}")
Wind speed at each level (mean m/s):
level 1000: 7.45
level 850: 9.14
level 700: 9.22
level 500: 12.39
level 400: 15.83
level 300: 19.13
Comparison operators¶
Comparison operators return a new Field containing boolean (0/1) values, which is useful as a mask.
[10]:
# Find grid points above freezing
above_freezing = t > 273.15
fraction = above_freezing.values.mean()
print(f"Fraction of grid points above 0°C: {fraction:.2%}")
# Use as a mask: zero out sub-freezing temperatures
t_warm = t * (t > 273.15)
print(f"Masked field min: {t_warm.values.min():.2f} K (should be 0 where frozen)")
Fraction of grid points above 0°C: 58.33%
Masked field min: 0.00 K (should be 0 where frozen)
Using numpy ufuncs¶
numpy universal functions work transparently on Field objects via earthkit.data.utils.compute.apply_ufunc().
You can also extract values with field.values, apply any numpy operation, and reconstruct
a field with field.set(values=...).
[11]:
# Using values directly with numpy, then setting back
log_t = t.set(values=np.log(t.values))
print(f"log(T) mean: {log_t.values.mean():.4f}")
# Verify the metadata is preserved
print(f"param: {log_t.parameter.variable()}, level: {log_t.vertical.level()}")
log(T) mean: 5.6312
param: t, level: 1000
apply_ufunc¶
apply_ufunc is earthkit-data’s mechanism for applying arbitrary callables — including numpy ufuncs and custom functions — directly to Field and FieldList objects without manually extracting and re-inserting values.
from earthkit.data.utils.compute import apply_ufunc
result = apply_ufunc(func, *args)
funcreceives plain numpy arrays and must return a numpy array.argscan beField,FieldList, scalar, or numpy array objects.earthkit-data iterates over matching fields, calls
funcon their value arrays, and wraps each result back into aField, copying metadata from the firstFieldListargument.If one
FieldListhas length 1, it is broadcast across all fields of the other.
[12]:
from earthkit.data.utils.compute import apply_ufunc
Unary functions¶
Pass a single Field or FieldList to apply any unary callable across every field.
[13]:
# np.sin applied to all temperature fields in one call
t_sin = apply_ufunc(np.sin, t_fields)
print(f"sin(T) fields: {len(t_sin)}")
for orig, result in zip(t_fields, t_sin):
print(f" level {orig.vertical.level():4d}: sin(T) mean = {result.values.mean():.6f}")
sin(T) fields: 6
level 1000: sin(T) mean = 0.398520
level 850: sin(T) mean = -0.029065
level 700: sin(T) mean = 0.161524
level 500: sin(T) mean = -0.067683
level 400: sin(T) mean = -0.096225
level 300: sin(T) mean = -0.189231
[14]:
# Custom unary function — clip temperatures to [250, 310] K
def clip_temp(arr):
return np.clip(arr, 250.0, 310.0)
t_clipped = apply_ufunc(clip_temp, t_fields)
for orig, clipped in zip(t_fields, t_clipped):
print(
f" level {orig.vertical.level():4d}: "
f"original [{orig.values.min():.1f}, {orig.values.max():.1f}] K -> "
f"clipped [{clipped.values.min():.1f}, {clipped.values.max():.1f}] K"
)
level 1000: original [240.6, 320.6] K -> clipped [250.0, 310.0] K
level 850: original [232.5, 304.5] K -> clipped [250.0, 304.5] K
level 700: original [223.3, 287.3] K -> clipped [250.0, 287.3] K
level 500: original [227.8, 271.8] K -> clipped [250.0, 271.8] K
level 400: original [216.0, 264.0] K -> clipped [250.0, 264.0] K
level 300: original [206.7, 250.7] K -> clipped [250.0, 250.7] K
Binary functions¶
Pass two Field / FieldList arguments for element-wise binary operations. np.hypot(u, v) computes wind speed more precisely than (u**2 + v**2)**0.5 because it avoids intermediate overflow.
[15]:
# Wind speed via np.hypot — avoids overflow for large values
wind_hypot = apply_ufunc(np.hypot, u_fields, v_fields)
wind_manual = (u_fields**2 + v_fields**2) ** 0.5
print("Wind speed comparison (hypot vs manual):")
for h, m in zip(wind_hypot, wind_manual):
print(f" level {h.vertical.level():4d}: hypot={h.values.mean():.4f} manual={m.values.mean():.4f}")
Wind speed comparison (hypot vs manual):
level 1000: hypot=7.4502 manual=7.4502
level 850: hypot=9.1436 manual=9.1436
level 700: hypot=9.2167 manual=9.2167
level 500: hypot=12.3901 manual=12.3901
level 400: hypot=15.8330 manual=15.8330
level 300: hypot=19.1302 manual=19.1302
[16]:
# Custom binary function — wind direction (meteorological convention, degrees from N)
def wind_direction(u_arr, v_arr):
# atan2 gives mathematical angle; convert to met convention (0=N, clockwise)
return (270.0 - np.degrees(np.arctan2(v_arr, u_arr))) % 360.0
wdir = apply_ufunc(wind_direction, u_fields, v_fields)
print("Wind direction at each level (mean degrees from N):")
for f in wdir:
print(f" level {f.vertical.level():4d}: {f.values.mean():.1f}°")
Wind direction at each level (mean degrees from N):
level 1000: 162.6°
level 850: 166.3°
level 700: 184.6°
level 500: 241.0°
level 400: 189.3°
level 300: 250.3°
Broadcasting a single field across a FieldList¶
When one argument has length 1 it is broadcast against every field in the other argument. Here we subtract the 1000 hPa temperature from every level to produce a temperature anomaly profile.
[17]:
t_1000 = t_fields.sel({"vertical.level": 1000}) # single-field FieldList (length 1)
print(f"Reference field count: {len(t_1000)}, all-level count: {len(t_fields)}")
# Broadcast: subtract t_1000 from every level
t_anomaly = apply_ufunc(lambda a, b: a - b, t_fields, t_1000)
print("\nTemperature anomaly relative to 1000 hPa (mean K):")
for orig, anom in zip(t_fields, t_anomaly):
print(f" level {orig.vertical.level():4d}: {anom.values.mean():+.2f} K")
Reference field count: 1, all-level count: 6
Temperature anomaly relative to 1000 hPa (mean K):
level 1000: +0.00 K
level 850: -6.98 K
level 700: -14.54 K
level 500: -26.82 K
level 400: -38.18 K
level 300: -51.86 K
[ ]: