API Reference

SatelliteGridding

A Julia package for gridding satellite Level-2 data onto regular rectangular grids with footprint-aware oversampling.

Satellite instruments have irregularly shaped footprints (quadrilaterals defined by 4 corner coordinates) that can overlap multiple output grid cells. Rather than computing exact geometric intersections, this package subdivides each footprint into n×n sub-pixels and uses "point-in-box" binning. The mean per grid cell is computed incrementally using Welford's online algorithm for numerical stability.

Main functions

• grid_l2: Grid Level-2 data with footprint oversampling
• grid_center: Grid data using center coordinates only (e.g., MODIS)
• load_config: Load a TOML/JSON data source configuration

Example

using SatelliteGridding

config = load_config("tropomi_sif.toml")
grid_spec = GridSpec(lat_min=-90f0, lat_max=90f0, lon_min=-180f0, lon_max=180f0,
                     dlat=1f0, dlon=1f0)
time_spec = TimeSpec(DateTime("2019-01-01"), DateTime("2019-12-31"),
                     Dates.Day(8), 1.0f0)
grid_l2(config, grid_spec, time_spec; outfile="output.nc")

grid(config, grid_spec, time_spec; method=SubpixelGridding(), kwargs...)
grid(config, grid_spec, time_spec, method::AbstractGriddingMethod; kwargs...)

Dispatch-based library entry point for gridding. The CLI still exposes explicit l2 and center commands, while library callers can choose a gridding method with a first-class type such as SubpixelGridding, CircularFootprintGridding, or CenterPointGridding.

grid_l2(config, grid_spec, time_spec; kwargs...)

Grid Level-2 satellite data with footprint-aware oversampling.

Reads input NetCDF files as specified by config, applies quality filters, subdivides each satellite footprint into sub-pixels for proper overlap handling, and accumulates gridded averages. Output is written to a NetCDF4 file. Quadrilateral footprints are used by default; pass footprint_method=CircularFootprintGridding(...) for circular products such as GOSAT. Circular footprints can be read from lat_bnd/lon_bnd or from center lat/lon plus radius.

Arguments

• config::DataSourceConfig: Data source configuration (from TOML/JSON)
• grid_spec::GridSpec: Output grid geometry
• time_spec::TimeSpec: Temporal binning parameters

Keyword Arguments

• n_oversample::Union{Nothing,Int}=nothing: Sub-pixel factor. nothing = auto-compute from footprint/grid ratio.
• footprint_method::AbstractGriddingMethod=SubpixelGridding(n_oversample): Footprint geometry method.
• compute_std::Bool=false: Also compute per-cell standard deviation
• outfile::String="gridded_output.nc": Output file path
• backend: KernelAbstractions backend (default: nothing for sequential Welford path). Use CPU(), CUDABackend(), MetalBackend(), or another compatible KA backend.
• keep_going::Bool=false: Continue after per-file errors. By default, the first file error stops the run.

grid_center(config, grid_spec, time_spec; kwargs...)

Grid data using center coordinates only (no footprint bounds), suitable for MODIS-style data where each pixel maps to exactly one grid cell.

Keyword Arguments

• geo_table::Union{Nothing,String}=nothing: Path to a legacy monolithic geolocation lookup table (NetCDF)
• geo_cache::Union{Nothing,String}=nothing: Directory for generated per-tile MODIS geolocation cache
• geo_provider::Union{Symbol,String}=:auto: :auto, :variables, :lut, or :modis
• veg_indices::Bool=false: Compute vegetation indices (EVI, NDVI, NIRv, NDWI)
• compute_std::Bool=false: Compute per-cell standard deviation
• outfile::String="gridded_output.nc": Output file path
• keep_going::Bool=false: Continue after per-file errors. By default, the first file error stops the run.

load_config(config_path::String) -> DataSourceConfig

Load a configuration file (TOML or JSON) and return a DataSourceConfig.

The file must contain a [grid] section plus filePattern and folder. Footprint-aware L2 gridding usually uses [basic] entries for lat_bnd and lon_bnd. Circular footprints can alternatively use center lat/lon plus a radius variable or scalar [circle] radius. Center-coordinate gridding can use [basic] lat and lon, a geolocation lookup table, or generated MODIS sinusoidal geolocation. Filters are defined in an optional [filter] section using intuitive string expressions:

Example TOML

filePattern = "S5P_PAL__L2B_SIF____YYYYMMDD*.nc"
folder = "/path/to/data/"

[basic]
lat = "PRODUCT/latitude"
lon = "PRODUCT/longitude"
lat_bnd = "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds"
lon_bnd = "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds"

[grid]
sif_743 = "PRODUCT/SIF_743"
sif_735 = "PRODUCT/SIF_735"

[filter]
"PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_zenith_angle" = "< 80"
"PRODUCT/qa_value" = "> 0.5"
"PRODUCT/methane" = "1600 < x < 2200"
"PRODUCT/mode" = "== 1"

The legacy JSON format with separate filter_gt/filter_lt/filter_eq sections is still supported for backward compatibility.

supported_backend_names() -> Tuple

Return the backend names accepted by resolve_backend and the CLI. sequential uses the single-threaded Welford path; the others use the KernelAbstractions accumulation path.

backend_help_lines() -> Vector{String}

Short descriptions for CLI help and documentation.

resolve_backend(name) -> Union{Nothing, KernelAbstractions.Backend}

Resolve a user-facing backend name to a KernelAbstractions backend object. The optional GPU packages are weak dependencies and are loaded only when the corresponding backend is requested.

cuda_backend()

Return a CUDA KernelAbstractions backend after CUDA.jl has been loaded.

metal_backend()

Return a Metal KernelAbstractions backend after Metal.jl has been loaded.

Back-compat alias for RectangularGridSpec. All call sites using GridSpec(...) and GridSpec{T} keep working unchanged.

TimeSpec

Specification of the temporal gridding parameters.

Fields

• start_date, stop_date: Date range for processing
• time_step: Temporal bin size (Dates.Day or Dates.Month)
• oversample_temporal: Multiplier for the actual averaging window (>1 gives moving-average-like behavior)

DataSourceConfig

Configuration loaded from a TOML/JSON file that defines a satellite data source.

Fields

• basic: Maps internal keys ("lat", "lon", "lat_bnd", "lon_bnd", "radius") to variable paths in the NetCDF files
• grid_vars: Ordered mapping of output variable names to input variable paths (all will be gridded)
• filters: Quality filter rules parsed from [filter] section
• file_pattern: Glob pattern with YYYY/MM/DD/DOY placeholders for finding input files
• folder: Root folder for input data (may also contain YYYY/MM/DD placeholders)
• options: Optional processing settings from [options], [center], or [modis]

FilterRule

A single filter condition for a NetCDF variable.

Fields

• variable: Path to the NetCDF variable (may include group separators /)
• op: Operation — :lt, :gt, :eq, :between
• lo: Lower bound (threshold for :gt/:lt/:eq, low end for :between)
• hi: Upper bound (only used for :between)

Config syntax (TOML [filter] section)

[filter]
"solar_zenith_angle" = "< 80"         # less than
"qa_value" = "> 0.5"                  # greater than
"mode" = "== 1"                       # equality
"methane" = "1600 < x < 2200"         # range (between)

AbstractGriddingMethod

Dispatch marker for choosing how observations are mapped onto the output grid. Use SubpixelGridding for footprint-aware oversampling, CircularFootprintGridding for circular or near-circular footprints, CenterPointGridding for fast center-coordinate binning, and ExactIntersectionGridding as the future exact geometry hook.

SubpixelGridding(; n_oversample=nothing)

Footprint-aware gridding by sampling an n × n grid of subpoints inside each quadrilateral footprint. nothing keeps the existing automatic oversampling heuristic.

CircularFootprintGridding(; n_oversample=nothing)

Footprint-aware gridding for circular or near-circular footprints described by a center coordinate and four bounding coordinates. The implementation samples an n × n grid over the footprint bounding box and keeps points inside the inferred circle or ellipse. nothing keeps the existing automatic oversampling heuristic.

CenterPointGridding()

Fast center-coordinate gridding. Each observation or pixel contributes to the single output cell containing its center coordinate.

ExactIntersectionGridding()

Reserved method for exact footprint/grid-cell intersection gridding. This is not implemented yet, but the type keeps the public API open for an exact geometry backend without changing callers.

sort_corners_ccw(lat_corners, lon_corners) -> (sorted_lat, sorted_lon)

Sort 4 footprint corners into counter-clockwise order by computing the angle from the centroid to each corner. This ensures edges 1→2 and 4→3 are always opposite edges, which is required by compute_subpixels!.

Different satellite products store corners in different orders (e.g., TROPOMI uses CCW order while OCO-2 uses a different convention). This function normalizes them.

sort_corners_ccw!(lat_mat, lon_mat)

Sort corners in-place for all soundings in N×4 matrices.

sort_corners_ccw_ka!(backend, lat_mat, lon_mat)

KA kernel that sorts footprint corners into CCW order for all soundings in parallel. Uses a 5-comparator sorting network (optimal for 4 elements, no allocations).

• lat_mat, lon_mat: N×4 matrices of corner coordinates

compute_subpixels!(points, vert_lat, vert_lon, n, lats_0, lons_0, lats_1, lons_1)

Fill the n×n×2 array points with sub-pixel positions inside the quadrilateral defined by 4 corner coordinates (vert_lat[1:4], vert_lon[1:4]).

The algorithm:

1. Subdivide edge 1→2 into n baseline points
2. Subdivide edge 4→3 into n baseline points
3. For each pair of corresponding baseline points, subdivide the connecting line into n interior points

This produces n² sub-pixels that fill the footprint quadrilateral. Temporary buffers lats_0, lons_0, lats_1, lons_1 (each length n) must be provided.

floor_indices!(ix, iy, points)

Convert the floating-point sub-pixel positions in points to integer grid cell indices via floor. Results stored in ix and iy (each length n²).

compute_n_oversample(footprint_extent, cell_size; n_min=2, n_max=20)

Estimate an appropriate oversampling factor n based on the ratio of footprint size to grid cell size. Aims for ~3 sub-pixels per grid cell width.

Returns an integer in [n_min, n_max].

compute_footprint_indices_ka!(backend, ix_out, iy_out, skip_flag,
                              lat_corners, lon_corners, n)

KA kernel that computes n×n sub-pixel grid indices for each footprint. Also sets skip_flag[fp] = true if the footprint's longitude extent >= n (too wide to oversample meaningfully).

For footprints where all 4 corners fall in the same cell, the kernel stores the single cell index in position 1 and marks all other positions with ix=-1 (sentinel). The accumulation kernel handles these fast-path footprints by checking skip_flag.

• ix_out: (n_fp, n²) Int32 — lat grid cell indices
• iy_out: (n_fp, n²) Int32 — lon grid cell indices
• skip_flag: (n_fp,) Int32 — 0=oversample, 1=fast-path (single cell), 2=skip (too wide)
• lat_corners, lon_corners: (n_fp, 4) — fractional grid indices of corners
• n: oversampling factor

compute_circular_footprint_indices_ka!(backend, ix_out, iy_out, inside_count,
                                       skip_flag, center_lat, center_lon,
                                       lat_corners, lon_corners, n)

KA kernel that computes sampled grid-cell indices for circular or near-circular footprints. The four coordinates define the footprint bounding box or edge points around the center coordinate. Points outside the inferred circle/ellipse are marked with -1; inside_count[fp] stores the number of valid samples used for normalizing the scatter weight.

polygon_area(points)

Compute planar polygon area with the shoelace formula. Points are (x, y) tuples in grid-index coordinates.

exact_footprint_weights(lat_corners, lon_corners; n_lon, n_lat)

Compute exact planar overlap weights between one quadrilateral footprint and regular grid cells in fractional grid-index coordinates. The returned matrix has shape (n_lon, n_lat) and sums to one when the footprint is fully inside the domain.

circle_footprint_weights(center_lat, center_lon, lat_corners, lon_corners;
                         n_lon, n_lat, n_oversample)

Compute sampled per-cell weights for one circular or near-circular footprint in fractional grid-index coordinates. The four coordinates define the footprint bounding box or edge points. In grid-index space the sampled shape may be an ellipse when the output grid has unequal latitude and longitude resolutions.

circle_footprint_weights(center_lat, center_lon, radius_lat, radius_lon;
                         n_lon, n_lat, n_oversample)

Convenience overload for a circular or elliptical footprint specified directly by center coordinate and radius in fractional grid-index units.

footprint_weights(method, lat_corners, lon_corners; n_lon, n_lat)

Return per-cell footprint weights for a single footprint using the requested geometry method. This is a lightweight geometry-level API used by tests and by future exact-intersection gridders.

accumulate_footprint!(grid_data, grid_std, grid_weights, compute_std,
                      lat_idx, lon_idx, values, n_pixels, n_vars, n_oversample,
                      grid_spec, points_buf, ix_buf, iy_buf,
                      lats0_buf, lons0_buf, lats1_buf, lons1_buf)

Accumulate a batch of satellite observations into the output grid using footprint-aware oversampling and Welford's online averaging algorithm.

For each pixel:

• Fast path: If all 4 corner coordinates fall into the same grid cell, the full observation value is assigned to that cell (weight=1).
• Oversample path: If the footprint spans multiple cells, it is subdivided into n_oversample × n_oversample sub-pixels. Each sub-pixel contributes weight 1/n² to whichever grid cell it falls in.

The Welford algorithm maintains a running mean and sum-of-squared-deviations (M2), enabling numerically stable computation of mean and variance in a single pass.

Arguments

• grid_data: Output array (n_lon, n_lat, n_vars) — running mean
• grid_std: Output array (n_lon, n_lat, n_vars) — running M2 (sum of squared deviations)
• grid_weights: Output array (n_lon, n_lat) — accumulated weights
• compute_std: Whether to track variance
• lat_idx, lon_idx: Pre-computed fractional grid indices for corners, size (n_pixels, 4)
• values: Input data, size (n_pixels, n_vars)
• n_pixels: Number of observations in this batch
• n_vars: Number of variables being gridded
• n_oversample: Sub-pixel subdivision factor (n)
• Remaining arguments are pre-allocated work buffers

accumulate_circular_footprint!(grid_data, grid_std, grid_weights, compute_std,
                               center_lat_idx, center_lon_idx,
                               lat_idx, lon_idx, values,
                               n_pixels, n_vars, n_oversample)

Accumulate observations with circular or near-circular footprints. The center coordinate gives the footprint center and the four coordinates define its bounding box or edge points. Sampling is done in fractional grid-index space.

accumulate_batch!(backend, grid_sum, grid_weights,
                  lat_corners, lon_corners, values, n;
                  compute_std=false, grid_std=nothing, grid_mean=nothing)

Accumulate observations into the grid using the KA kernel pipeline. This function works identically on CPU and GPU backends.

The pipeline:

1. Sort corners into CCW order
2. Compute n×n subpixel grid indices per footprint
3. Scatter-accumulate weighted sums into grid cells (using atomics)

After processing all files, call finalize_mean! to convert sums to means.

For compute_std=true, a two-pass approach is used:

1. First call accumulate_batch! without std to accumulate sums
2. Call finalize_mean! to get the mean
3. Call accumulate_std_pass! with the computed mean to accumulate variance

Arguments

• backend: KernelAbstractions backend (CPU(), CUDABackend(), MetalBackend(), etc.)
• grid_sum: (nlon, nlat, n_vars) — weighted sum accumulator
• grid_weights: (nlon, nlat) — weight accumulator
• lat_corners, lon_corners: (n_fp, 4) — fractional grid indices of corners
• values: (nfp, nvars) — observation values
• n: oversampling factor

accumulate_circular_batch!(backend, grid_sum, grid_weights,
                           center_lat, center_lon,
                           lat_corners, lon_corners, values, n, n_vars)

Accumulate circular or near-circular footprints using the KA pipeline. This is the sum-based backend counterpart to accumulate_circular_footprint! and works with KA CPU, CUDA, Metal, and other compatible backends.

accumulate_center!(grid_data, grid_weights, lat, lon, values, grid_spec)

Accumulate observations using center coordinates only (no footprint bounds). Used for MODIS-style gridding where each pixel maps to exactly one grid cell. Simple sum accumulation — divide by grid_weights afterward to get the mean.

scatter_accumulate!(grid_sum, grid_weights, ix, iy, skip_flag, values, n, n_vars)

Sequential scatter-accumulate of weighted values into grid cells. Uses precomputed grid indices from compute_footprint_indices_ka!.

For each footprint:

• skip_flag == 1 (fast path): add weight=1, value to single cell
• skip_flag == 0 (oversample): add weight=1/n² to each subpixel's cell
• skip_flag == 2 (skip): do nothing

This function is sequential (no atomics needed) and works with any array type. For GPU backends, use the @kernel version with @atomic instead.

• grid_sum: (nlon, nlat, n_vars) — weighted sum accumulator
• grid_weights: (nlon, nlat) — weight accumulator
• ix, iy: (n_fp, n²) — precomputed grid indices
• skip_flag: (n_fp,) — 0=oversample, 1=fast-path, 2=skip
• values: (nfp, nvars) — observation values

scatter_accumulate_ka!(backend, grid_sum, grid_weights,
                       ix, iy, skip_flag, values, n, n_vars)

KA kernel version of scatter-accumulate using @atomic for thread safety. Required for GPU backends where multiple threads write to the same grid cells. For CPU backends, use scatter_accumulate! (sequential) instead.

scatter_accumulate_circular!(grid_sum, grid_weights, ix, iy, inside_count,
                             skip_flag, values, n, n_vars)

Sequential scatter-accumulate for circular footprints. Oversampled footprints use per-footprint weights 1 / inside_count[fp] because samples outside the circle/ellipse are masked out.

scatter_accumulate_circular_ka!(backend, grid_sum, grid_weights,
                                ix, iy, inside_count, skip_flag,
                                values, n, n_vars)

KA kernel version of circular-footprint scatter-accumulate using atomics for GPU backends.

finalize_mean!(grid_data, grid_weights)

Convert weighted-sum accumulator to mean: mean = sum / weight. Cells with near-zero weight are left as zero.

find_files(pattern::String, folder::String, date::DateTime) -> Vector{String}

Find input files matching the given pattern for a specific date.

Resolves YYYY, YY, MM, DD, and DOY placeholders in both pattern and folder, then globs for matching files. Returns only non-empty files.

read_nc_variable(dataset, path::String; bounds=false)

Read a variable from a NetCDF dataset, handling nested groups via / separators.

If bounds=true, reshape the result so that footprint vertices are in the last dimension, returning an N×4 matrix (N soundings, 4 corners).

read_nc_attribute(dataset, path::String, attr::String)

Read an attribute from a possibly group-nested variable.

read_variable_from_file(path, variable; bounds=false)

Read a variable from a NetCDF-like file. Falls back to the system libnetcdf for flat HDF4/HDF-EOS variables when the Julia NetCDF build cannot open them.

read_array_from_file(path, variable)

Read a variable while preserving its native dimensionality.

read_system_netcdf_variable(path, varname)

Read a flat variable through the system libnetcdf. This is primarily a fallback for HDF4/HDF-EOS2 files when the Julia NetCDF_jll build lacks HDF4 support.

create_output_dataset(outfile, grid_spec, time_spec, grid_vars, compute_std) -> Dataset, Dict

Create the output NetCDF4 file with dimensions, coordinate variables, and data variable definitions. Returns the dataset and a Dict mapping variable names to NCDatasets variables.

default_cache_dir() -> String

Return the base cache directory for generated auxiliary data. Override with SATELLITEGRIDDING_CACHE.

default_modis_cache_dir(; pixels=2400) -> String

Return the default per-tile MODIS sinusoidal geolocation cache directory. Override with SATELLITEGRIDDING_MODIS_GEO_CACHE.

parse_modis_tile(path_or_tile) -> (h, v)

Parse a MODIS tile identifier such as h08v04 from a filename or string.

modis_sinusoidal_latlon(h, v, line, sample; pixels=2400)

Compute the latitude/longitude at a MODLAND sinusoidal pixel center. h and v are zero-based tile indices; line and sample are one-based pixel indices within the tile.

generate_modis_tile_geolocation(h, v; pixels=2400)

Generate latitude and longitude arrays for a MODLAND sinusoidal tile from the official projection definition. Arrays are indexed as [line, sample].

write_modis_tile_geolocation(path, h, v; pixels=2400, overwrite=false)

Generate and write one MODIS sinusoidal tile geolocation file. The output is a small per-tile NetCDF file with latitude(y, x) and longitude(y, x).

load_or_generate_modis_tile_geolocation(h, v; cache_dir=default_modis_cache_dir(), pixels=2400)

Read a per-tile MODIS geolocation cache file, generating it first if needed. Returns (latitude, longitude) arrays indexed as [line, sample].

apply_filters(dataset, config::DataSourceConfig, lat_bounds, lon_bounds,
              grid_spec::GridSpec) -> Vector{Int}

Apply spatial bounding box and quality filters, returning indices of pixels that pass all criteria.

The spatial filter requires all 4 corners to be within the grid bounds and the footprint to not cross the dateline (lon span < 50°). Quality filters from config.filters are applied as additional criteria.

compute_evi(red, nir, blue; G=2.5f0, C1=6.0f0, C2=7.5f0, L=10000.0f0)

Compute the Enhanced Vegetation Index (EVI).

EVI = G × (NIR - RED) / (NIR + C1×RED - C2×BLUE + L)

compute_ndvi(red, nir)

Compute the Normalized Difference Vegetation Index (NDVI).

NDVI = (NIR - RED) / (NIR + RED)

compute_nirv(red, nir)

Compute the Near-Infrared Reflectance of Vegetation (NIRv).

NIRv = NDVI × NIR

compute_ndwi(nir, swir)

Compute the Normalized Difference Water Index (NDWI).

NDWI = (NIR - SWIR) / (NIR + SWIR)

parse_l2_args(args=ARGS) -> Dict

Parse command-line arguments for L2 gridding. Matches the interface of the original gridL2_Dates.jl script.

parse_center_args(args=ARGS) -> Dict

Parse command-line arguments for center-coordinate gridding (MODIS-style).

args_to_grid_spec(args::Dict) -> GridSpec{Float32}

Convert parsed CLI arguments to a GridSpec.

args_to_time_spec(args::Dict) -> TimeSpec

Convert parsed CLI arguments to a TimeSpec.