Priors

The prior mean emission from a gridded inventory (EDGAR, GFEI, EPA), averaged over the analysis T-cells, and the v1 per-tile background for column observations. Tier 0. load_inventory needs pip install 'plumax[io]'.

priors

Priors for the inversions (Tier 0 prerequisites).

• :func:load_inventory — a gridded inventory's area-mean flux on the analysis T-cells, as the prior mean emission [kg m⁻² s⁻¹].
• :func:per_tile_background — v1 background for column observations, the median within square tiles of the analysis frame.

per_tile_background

per_tile_background(columns: Float[ndarray, ' obs'], receptors: Float[ndarray, 'obs 2'], *, tile_m: float = 10000.0) -> Float[np.ndarray, ' obs']

Median of the observations in each square tile, given to every receptor.

b_n = median{ y_m : tile(𝐱_m) = tile(𝐱_n) }, tile(𝐱) = (⌊x / ℓ⌋, ⌊y / ℓ⌋).

The median is not moved by a plume covering less than half a tile.

Parameters:

Name	Type	Description	Default
columns	Float[ndarray, ' obs']	Observed columns y_n, any units.	required
receptors	Float[ndarray, 'obs 2']	Receptor positions 𝐱_n = (x, y) in frame metres.	required
tile_m	float	Tile edge ℓ [m].	10000.0

Returns:

Type	Description
Float[ndarray, ' obs']	Background b_n per receptor, in the units of columns.

Raises:

Type	Description
GridMismatchError	If receptors does not match columns.
PlumaxError	If tile_m is not positive.

Examples:

Two tiles, one with a plume pixel the median ignores:

>>> import numpy as np
>>> receptors = np.array(
...     [[1e3, 0.0], [2e3, 0.0], [3e3, 0.0], [12e3, 0.0], [15e3, 0.0]]
... )
>>> columns = np.array([1880.0, 1880.0, 1960.0, 1900.0, 1900.0])
>>> per_tile_background(columns, receptors).tolist()
[1880.0, 1880.0, 1880.0, 1900.0, 1900.0]

One tile spanning the scene gives the scene median:

>>> per_tile_background(columns, receptors, tile_m=1e6).tolist()
[1900.0, 1900.0, 1900.0, 1900.0, 1900.0]

Source code in src/plumax/priors/_src/background.py

def per_tile_background(
    columns: Float[np.ndarray, " obs"],
    receptors: Float[np.ndarray, "obs 2"],
    *,
    tile_m: float = 10_000.0,
) -> Float[np.ndarray, " obs"]:
    """Median of the observations in each square tile, given to every receptor.

    b_n = median{ y_m : tile(𝐱_m) = tile(𝐱_n) },  tile(𝐱) = (⌊x / ℓ⌋, ⌊y / ℓ⌋).

    The median is not moved by a plume covering less than half a tile.

    Args:
        columns: Observed columns y_n, any units.
        receptors: Receptor positions 𝐱_n = (x, y) in frame metres.
        tile_m: Tile edge ℓ [m].

    Returns:
        Background b_n per receptor, in the units of ``columns``.

    Raises:
        GridMismatchError: If ``receptors`` does not match ``columns``.
        PlumaxError: If ``tile_m`` is not positive.

    Examples:
        Two tiles, one with a plume pixel the median ignores:

        >>> import numpy as np
        >>> receptors = np.array(
        ...     [[1e3, 0.0], [2e3, 0.0], [3e3, 0.0], [12e3, 0.0], [15e3, 0.0]]
        ... )
        >>> columns = np.array([1880.0, 1880.0, 1960.0, 1900.0, 1900.0])
        >>> per_tile_background(columns, receptors).tolist()
        [1880.0, 1880.0, 1880.0, 1900.0, 1900.0]

        One tile spanning the scene gives the scene median:

        >>> per_tile_background(columns, receptors, tile_m=1e6).tolist()
        [1900.0, 1900.0, 1900.0, 1900.0, 1900.0]
    """
    values = np.asarray(columns, dtype=float)
    positions = np.asarray(receptors, dtype=float)
    if positions.shape != (values.size, 2):
        raise GridMismatchError(
            f"per_tile_background: receptors must be ({values.size}, 2); "
            f"got {positions.shape}"
        )
    if not tile_m > 0.0:
        raise PlumaxError("per_tile_background: `tile_m` must be positive")
    # tile(𝐱) = (⌊x / ℓ⌋, ⌊y / ℓ⌋)
    tiles = np.floor(positions / tile_m).astype(np.int64)
    tile_x, tile_y = einx.id("obs c -> c obs", tiles)
    # one integer per tile, so a 1-D unique groups receptors by tile
    tile_key = tile_x * (2**32) + tile_y
    _, tile_label = np.unique(tile_key, return_inverse=True)
    background = np.empty_like(values)
    for label in np.unique(tile_label):
        in_tile = tile_label == label
        # b = median of the tile's observations
        background[in_tile] = np.median(values[in_tile])
    return background

load_inventory

load_inventory(path: str | Path, *, frame: LocalFrame, grid: PlumeGrid3D, variable: str | None = None, units: str | None = None, time_index: int = 0, supersample: int = 8) -> Float[Array, 'y x']

Area-mean inventory flux over each analysis T-cell, as a prior mean.

F̄c = (1/|c|) ∫_c F dA ≈ (1/s²) Σ F(𝐱_k),}^{s²

with 𝐱_k an s × s lattice of points evenly spaced in frame metres over cell c. The sampling error falls as 1/s at inventory-cell edges; the frame is conformal and locally near-equidistant, so equal frame areas are equal ground areas.

Algorithm

1. Read the inventory with :func:xrreader.open_ch4_inventory, which converts to kg m⁻² s⁻¹ with ascending lat and lon in [-180, 180).
2. Lay an s × s lattice over every analysis cell and map it to lat/lon.
3. Unwrap the raster's longitudes to start after their widest gap, look up the inventory cell holding each point (longitudes modulo 360° from the first edge) and average over each analysis cell's lattice.

Parameters:

Name	Type	Description	Default
path	str | Path	Inventory NetCDF (one sector or a total).	required
frame	LocalFrame	Frame the analysis grid is laid out in.	required
grid	PlumeGrid3D	Analysis grid, in frame metres.	required
variable	str | None	Flux variable name; inferred by the reader when None.	None
units	str | None	Override of the file's units attribute.	None
time_index	int	Time step to use when the file has a time axis.	0
supersample	int	Sub-points per cell edge, s.	8

Returns:

Type	Description
Float[Array, 'y x']	F̄ [kg m⁻² s⁻¹] on the horizontal T-cells, in the grid's dtype.

Raises:

Type	Description
MissingExtraError	If the io extra (xrreader) is not installed.
DomainError	If a lattice point falls outside the inventory raster, or the raster's longitudes are unevenly spaced.

Source code in src/plumax/priors/_src/inventory.py

def load_inventory(
    path: str | Path,
    *,
    frame: LocalFrame,
    grid: PlumeGrid3D,
    variable: str | None = None,
    units: str | None = None,
    time_index: int = 0,
    supersample: int = 8,
) -> Float[Array, "y x"]:
    """Area-mean inventory flux over each analysis T-cell, as a prior mean.

    F̄_c = (1/|c|) ∫_c F dA  ≈  (1/s²) Σ_{k=1}^{s²} F(𝐱_k),

    with ``𝐱_k`` an ``s × s`` lattice of points evenly spaced in frame metres
    over cell ``c``. The sampling error falls as 1/s at inventory-cell edges;
    the frame is conformal and locally near-equidistant, so equal frame areas
    are equal ground areas.

    Algorithm:
        1. Read the inventory with :func:`xrreader.open_ch4_inventory`, which
           converts to kg m⁻² s⁻¹ with ascending lat and lon in [-180, 180).
        2. Lay an ``s × s`` lattice over every analysis cell and map it to
           lat/lon.
        3. Unwrap the raster's longitudes to start after their widest gap,
           look up the inventory cell holding each point (longitudes modulo
           360° from the first edge) and average over each analysis cell's
           lattice.

    Args:
        path: Inventory NetCDF (one sector or a total).
        frame: Frame the analysis grid is laid out in.
        grid: Analysis grid, in ``frame`` metres.
        variable: Flux variable name; inferred by the reader when ``None``.
        units: Override of the file's units attribute.
        time_index: Time step to use when the file has a ``time`` axis.
        supersample: Sub-points per cell edge, s.

    Returns:
        F̄ [kg m⁻² s⁻¹] on the horizontal T-cells, in the grid's dtype.

    Raises:
        MissingExtraError: If the ``io`` extra (xrreader) is not installed.
        DomainError: If a lattice point falls outside the inventory raster,
            or the raster's longitudes are unevenly spaced.
    """
    try:
        import xrreader
    except ImportError as err:
        raise MissingExtraError(
            "load_inventory reads through xrreader; install it with "
            "`pip install 'plumax[io]'`"
        ) from err

    # --- read ---
    dataset = xrreader.open_ch4_inventory(path, variable=variable, units=units)
    flux = dataset["ch4_emission_flux"]
    if "time" in flux.dims:
        flux = flux.isel(time=time_index)
    flux = flux.transpose("lat", "lon")
    lat_edges = _edges(flux["lat"].values)
    lon_centres, flux_values = _unwrap_longitudes(flux["lon"].values, flux.values)
    lon_edges = _edges(lon_centres)

    # --- sub-cell lattice ---
    x = np.asarray(grid.x, dtype=float)
    y = np.asarray(grid.y, dtype=float)
    # u_k = (k + 1/2) / s − 1/2, the sub-cell centres in cell units
    unit_offsets = (np.arange(supersample) + 0.5) / supersample - 0.5
    # x_{i,k} = x_i + u_k Δx
    sub_x = einx.add("x, k -> (x k)", x, unit_offsets * grid.dx)
    # y_{j,k} = y_j + u_k Δy
    sub_y = einx.add("y, k -> (y k)", y, unit_offsets * grid.dy)
    lattice_x, lattice_y = np.meshgrid(sub_x, sub_y)
    lat, lon = frame.to_latlon(lattice_x, lattice_y)

    # --- look up and average ---
    # λ ← λ₀ + ((λ − λ₀) mod 360°): sample longitudes on the raster's
    # unwrapped axis; a global raster then covers every point
    lon = lon_edges[0] + np.mod(lon - lon_edges[0], 360.0)
    row = np.searchsorted(lat_edges, lat, side="right") - 1
    col = np.searchsorted(lon_edges, lon, side="right") - 1
    outside_rows = (row < 0) | (row >= lat_edges.size - 1)
    outside_cols = (col < 0) | (col >= lon_edges.size - 1)
    if np.any(outside_rows | outside_cols):
        raise DomainError(
            "load_inventory: the analysis grid extends beyond the inventory raster"
        )
    # F(𝐱_k)
    samples = einx.get_at("[lat lon], ys xs, ys xs -> ys xs", flux_values, row, col)
    # F̄_c = mean over the s × s lattice of cell c
    cell_mean = einx.mean(
        "(y [sy]) (x [sx]) -> y x", samples, sy=supersample, sx=supersample
    )
    return jnp.asarray(cell_mean, dtype=grid.x.dtype)