Skip to content

Domain

domain

Geographic domains and regional cubes of gridded fields.

Cube dataclass

One hour of gridded fields on a regional grid, levels top to bottom.

Source code in packages/barograph-met/src/baro/met/domain.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
@dataclass
class Cube:
    """One hour of gridded fields on a regional grid, levels top to bottom."""

    source: str
    hour: np.datetime64
    lat: NDArray[np.float64]  # descending
    lon: NDArray[np.float64]  # unwrapped, ascending
    vertical: VerticalCoordinate
    fields: dict[str, NDArray[np.float32]]  # (level, lat, lon)
    ps: NDArray[np.float32] | None = None  # surface pressure (Pa), (lat, lon)
    run: np.datetime64 | None = None  # forecast base time (None for analyses)

    def __post_init__(self) -> None:
        shape = (len(self.vertical), len(self.lat), len(self.lon))
        for name, arr in self.fields.items():
            if arr.shape != shape:
                raise ValueError(f"field {name} has shape {arr.shape}, expected {shape}")
        if self.ps is not None and self.ps.shape != shape[1:]:
            raise ValueError("surface pressure grid does not match the fields")
        if len(self.lat) < _MIN_GRID_POINTS or len(self.lon) < _MIN_GRID_POINTS:
            raise ValueError("a cube needs at least 2x2 grid points")

    @property
    def nbytes(self) -> int:
        """Memory held by the fields and surface pressure."""
        total = sum(a.nbytes for a in self.fields.values())
        return total + (0 if self.ps is None else self.ps.nbytes)
nbytes property

Memory held by the fields and surface pressure.

Domain dataclass

Latitude/longitude box; longitudes may be expressed beyond 180 (no wrap inside).

Source code in packages/barograph-met/src/baro/met/domain.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
@dataclass(frozen=True)
class Domain:
    """Latitude/longitude box; longitudes may be expressed beyond 180 (no wrap inside)."""

    lat_min: float
    lat_max: float
    lon_min: float
    lon_max: float

    def __post_init__(self) -> None:
        if not (-LAT_LIMIT <= self.lat_min <= self.lat_max <= LAT_LIMIT):
            raise ValueError(f"invalid latitude range {self.lat_min}..{self.lat_max}")
        if not (0.0 <= self.lon_max - self.lon_min <= FULL_TURN):
            raise ValueError(f"invalid longitude range {self.lon_min}..{self.lon_max}")

    @classmethod
    def from_points(cls, lat: ArrayLike, lon: ArrayLike, margin: float = 0.5) -> Domain:
        """Smallest box around the points (shortest longitude arc), plus ``margin`` degrees."""
        la = np.asarray(lat, dtype=np.float64)
        lo = np.asarray(lon, dtype=np.float64)
        if la.size == 0:
            raise ValueError("no points")
        best: tuple[float, float] | None = None
        for base in (-180.0, 0.0):  # [-180, 180) and [0, 360) framings
            framed = base + np.mod(lo - base, 360.0)
            span = (float(framed.min()), float(framed.max()))
            if best is None or span[1] - span[0] < best[1] - best[0]:
                best = span
        assert best is not None
        lon_min, lon_max = best[0] - margin, best[1] + margin
        if lon_max - lon_min > FULL_TURN:
            lon_min, lon_max = -180.0, 180.0
        return cls(
            max(float(la.min()) - margin, -90.0),
            min(float(la.max()) + margin, 90.0),
            lon_min,
            lon_max,
        )

    def to_json(self) -> dict[str, float]:
        """Bounds as a JSON-friendly dict."""
        return {k: getattr(self, k) for k in ("lat_min", "lat_max", "lon_min", "lon_max")}

    @classmethod
    def whole_globe(cls) -> Domain:
        """The domain covering the entire globe."""
        return cls(-90.0, 90.0, 0.0, 360.0)
from_points(lat, lon, margin=0.5) classmethod

Smallest box around the points (shortest longitude arc), plus margin degrees.

Source code in packages/barograph-met/src/baro/met/domain.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
@classmethod
def from_points(cls, lat: ArrayLike, lon: ArrayLike, margin: float = 0.5) -> Domain:
    """Smallest box around the points (shortest longitude arc), plus ``margin`` degrees."""
    la = np.asarray(lat, dtype=np.float64)
    lo = np.asarray(lon, dtype=np.float64)
    if la.size == 0:
        raise ValueError("no points")
    best: tuple[float, float] | None = None
    for base in (-180.0, 0.0):  # [-180, 180) and [0, 360) framings
        framed = base + np.mod(lo - base, 360.0)
        span = (float(framed.min()), float(framed.max()))
        if best is None or span[1] - span[0] < best[1] - best[0]:
            best = span
    assert best is not None
    lon_min, lon_max = best[0] - margin, best[1] + margin
    if lon_max - lon_min > FULL_TURN:
        lon_min, lon_max = -180.0, 180.0
    return cls(
        max(float(la.min()) - margin, -90.0),
        min(float(la.max()) + margin, 90.0),
        lon_min,
        lon_max,
    )
to_json()

Bounds as a JSON-friendly dict.

Source code in packages/barograph-met/src/baro/met/domain.py 
59
60
61
def to_json(self) -> dict[str, float]:
    """Bounds as a JSON-friendly dict."""
    return {k: getattr(self, k) for k in ("lat_min", "lat_max", "lon_min", "lon_max")}
whole_globe() classmethod

The domain covering the entire globe.

Source code in packages/barograph-met/src/baro/met/domain.py 
63
64
65
66
@classmethod
def whole_globe(cls) -> Domain:
    """The domain covering the entire globe."""
    return cls(-90.0, 90.0, 0.0, 360.0)

GridIndex dataclass

Where a domain sits in a regular grid (latitude descending, longitude ascending).

A grid spanning 360 degrees of longitude is periodic (columns wrap around); a regional grid is clipped: the part of the domain outside it is dropped and points there are later flagged OUTSIDE_DOMAIN.

Source code in packages/barograph-met/src/baro/met/domain.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
@dataclass(frozen=True)
class GridIndex:
    """Where a domain sits in a regular grid (latitude descending, longitude ascending).

    A grid spanning 360 degrees of longitude is periodic (columns wrap around);
    a regional grid is clipped: the part of the domain outside it is dropped and
    points there are later flagged ``OUTSIDE_DOMAIN``.
    """

    rows: slice
    cols: NDArray[np.int64]  # indices into the global longitude axis (may wrap)
    lat: NDArray[np.float64]
    lon: NDArray[np.float64]  # unwrapped, strictly increasing

    @classmethod
    def locate(
        cls, domain: Domain, grid_lat: NDArray[np.float64], grid_lon: NDArray[np.float64]
    ) -> GridIndex:
        """Rows, columns and coordinates of ``domain`` in a regular grid."""
        dlat = float(grid_lat[1] - grid_lat[0])
        dlon = float(grid_lon[1] - grid_lon[0])
        if dlat >= 0 or dlon <= 0:
            raise ValueError("expected descending latitudes and ascending longitudes")
        nlat, nlon = len(grid_lat), len(grid_lon)
        i0 = int(np.floor((domain.lat_max - grid_lat[0]) / dlat + 1e-9))
        i1 = int(np.ceil((domain.lat_min - grid_lat[0]) / dlat - 1e-9))
        # A wholly outside regional domain still needs a valid two-row sliver:
        # sampling flags its points rather than indexing beyond the archive.
        i0 = min(max(i0, 0), nlat - 2)
        i1 = min(max(i1, i0 + 1), nlat - 1)
        periodic = abs(nlon * dlon - FULL_TURN) < _PERIODIC_TOLERANCE
        lon0 = float(grid_lon[0])
        # periodic: keep the caller's framing (columns wrap); regional: the framing
        # of the domain that overlaps the grid most
        start = domain.lon_min if periodic else _regional_start(domain, lon0, dlon, nlon)
        stop = start + (domain.lon_max - domain.lon_min)
        j0 = int(np.floor((start - lon0) / dlon + 1e-9))
        j1 = int(np.ceil((stop - lon0) / dlon - 1e-9))
        if periodic:
            j1 = min(max(j1, j0 + 1), j0 + nlon)
        else:
            j0, j1 = max(j0, 0), min(j1, nlon - 1)
            if j1 - j0 < 1:
                # Inside: keep the location's neighbour. Outside: use the
                # nearest edge sliver so sampling can flag the point.
                j0 = min(max(j0, 0), nlon - 2)
                j1 = j0 + 1
        j = np.arange(j0, j1 + 1, dtype=np.int64)
        return cls(
            rows=slice(i0, i1 + 1),
            cols=np.mod(j, nlon),
            lat=np.asarray(grid_lat[i0 : i1 + 1], dtype=np.float64),
            lon=(float(grid_lon[0]) + j * dlon).astype(np.float64),
        )
locate(domain, grid_lat, grid_lon) classmethod

Rows, columns and coordinates of domain in a regular grid.

Source code in packages/barograph-met/src/baro/met/domain.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
@classmethod
def locate(
    cls, domain: Domain, grid_lat: NDArray[np.float64], grid_lon: NDArray[np.float64]
) -> GridIndex:
    """Rows, columns and coordinates of ``domain`` in a regular grid."""
    dlat = float(grid_lat[1] - grid_lat[0])
    dlon = float(grid_lon[1] - grid_lon[0])
    if dlat >= 0 or dlon <= 0:
        raise ValueError("expected descending latitudes and ascending longitudes")
    nlat, nlon = len(grid_lat), len(grid_lon)
    i0 = int(np.floor((domain.lat_max - grid_lat[0]) / dlat + 1e-9))
    i1 = int(np.ceil((domain.lat_min - grid_lat[0]) / dlat - 1e-9))
    # A wholly outside regional domain still needs a valid two-row sliver:
    # sampling flags its points rather than indexing beyond the archive.
    i0 = min(max(i0, 0), nlat - 2)
    i1 = min(max(i1, i0 + 1), nlat - 1)
    periodic = abs(nlon * dlon - FULL_TURN) < _PERIODIC_TOLERANCE
    lon0 = float(grid_lon[0])
    # periodic: keep the caller's framing (columns wrap); regional: the framing
    # of the domain that overlaps the grid most
    start = domain.lon_min if periodic else _regional_start(domain, lon0, dlon, nlon)
    stop = start + (domain.lon_max - domain.lon_min)
    j0 = int(np.floor((start - lon0) / dlon + 1e-9))
    j1 = int(np.ceil((stop - lon0) / dlon - 1e-9))
    if periodic:
        j1 = min(max(j1, j0 + 1), j0 + nlon)
    else:
        j0, j1 = max(j0, 0), min(j1, nlon - 1)
        if j1 - j0 < 1:
            # Inside: keep the location's neighbour. Outside: use the
            # nearest edge sliver so sampling can flag the point.
            j0 = min(max(j0, 0), nlon - 2)
            j1 = j0 + 1
    j = np.arange(j0, j1 + 1, dtype=np.int64)
    return cls(
        rows=slice(i0, i1 + 1),
        cols=np.mod(j, nlon),
        lat=np.asarray(grid_lat[i0 : i1 + 1], dtype=np.float64),
        lon=(float(grid_lon[0]) + j * dlon).astype(np.float64),
    )