GitHub

Tutorial: Gridding TROPOMI SIF Data

This tutorial walks through gridding TROPOMI Solar-Induced Fluorescence (SIF) data onto a regular lat/lon grid using SatelliteGridding.jl.

Setup

First, load the package and define the time range:

using SatelliteGridding
using Dates

Step 1: Load the Configuration

Configuration files (TOML format) tell the package where to find data files, which NetCDF variables to read, and what quality filters to apply.

Here's what a typical TROPOMI SIF config looks like:

filePattern = "S5P_PAL__L2B_SIF____YYYYMMDD*.nc"
folder = "/path/to/tropomi/sif/"

[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"
config = load_config(joinpath(pkgdir(SatelliteGridding), "examples", "tropomi_sif.toml"))
DataSourceConfig(Dict("lat" => "PRODUCT/latitude", "lon_bnd" => "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds", "lat_bnd" => "PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds", "lon" => "PRODUCT/longitude"), OrderedCollections.OrderedDict("sif_735_corr" => "PRODUCT/SIF_Corr_735", "sif_743" => "PRODUCT/SIF_743", "sif_735" => "PRODUCT/SIF_735", "cloud_fraction" => "PRODUCT/SUPPORT_DATA/INPUT_DATA/cloud_fraction_L2", "sif_743_corr" => "PRODUCT/SIF_Corr_743", "sif_743_err" => "PRODUCT/SIF_ERROR_743"), FilterRule[FilterRule("PRODUCT/SUPPORT_DATA/GEOLOCATIONS/solar_zenith_angle", :lt, 80.0, NaN)], "S5P_PAL__L2B_SIF____YYYYMMDD*.nc", "/kiwi-data/Data/satellite/TROPOMI/TROPOMI_SIF_S5P-PAL/original/", Dict{String, Any}())

The config object contains:

  • basic: Coordinate variable paths (lat, lon, and corner bounds)
  • grid_vars: Variables to be gridded (e.g., SIF at 743nm and 735nm)
  • filters: Quality filter rules (here: solar zenith angle < 80°)
  • file_pattern: Glob pattern with date placeholders
  • folder: Root data directory

Step 2: Define the Output Grid

A GridSpec defines the spatial extent and resolution of the output grid. All values are in degrees. The package uses Float32 for memory efficiency.

grid_spec = GridSpec(
    lat_min = -60f0,    # Southern bound
    lat_max =  80f0,    # Northern bound
    lon_min = -180f0,   # Western bound
    lon_max =  180f0,   # Eastern bound
    dlat = 0.5f0,       # Latitude resolution (degrees)
    dlon = 0.5f0        # Longitude resolution (degrees)
)
RectangularGridSpec{Float32}(-60.0f0, 80.0f0, -180.0f0, 180.0f0, 0.5f0, 0.5f0, Float32[-59.75, -59.25, -58.75, -58.25, -57.75, -57.25, -56.75, -56.25, -55.75, -55.25  …  75.25, 75.75, 76.25, 76.75, 77.25, 77.75, 78.25, 78.75, 79.25, 79.75], Float32[-179.75, -179.25, -178.75, -178.25, -177.75, -177.25, -176.75, -176.25, -175.75, -175.25  …  175.25, 175.75, 176.25, 176.75, 177.25, 177.75, 178.25, 178.75, 179.25, 179.75])

This creates a 720×280 grid. The grid_spec.lat and grid_spec.lon vectors contain cell center coordinates.

println("Grid dimensions: $(length(grid_spec.lon)) × $(length(grid_spec.lat))")
Grid dimensions: 720 × 280

Step 3: Define the Time Specification

A TimeSpec controls temporal binning. You can use Dates.Day or Dates.Month for the time step.

time_spec = TimeSpec(
    DateTime("2019-07-01"),     # Start date
    DateTime("2019-07-31"),     # Stop date
    Dates.Day(16)               # 16-day composites
)
TimeSpec(Dates.DateTime("2019-07-01T00:00:00"), Dates.DateTime("2019-07-31T00:00:00"), Dates.Day(16), 1.0f0, nothing, nothing)

For monthly composites:

time_spec = TimeSpec(DateTime("2019-01-01"), DateTime("2019-12-31"),
                     Dates.Month(1))

Step 4: Run the Gridding

The grid_l2 function processes all matching files, applies filters, subdivides footprints, and writes the output NetCDF file.

# grid_l2(config, grid_spec, time_spec;
#         outfile = "tropomi_sif_july2019.nc",
#         n_oversample = 10)    # 10×10 sub-pixels per footprint

Key Options

  • n_oversample: Sub-pixel subdivision factor. Higher = more accurate spatial distribution but slower. Default: auto-computed from footprint/grid ratio.
  • compute_std: Set to true to also compute per-cell standard deviation.
  • backend: Compute backend — nothing (sequential Welford), CPU() (KA parallel), CUDABackend() (NVIDIA GPU), MetalBackend() (Apple GPU), or another compatible KernelAbstractions backend.

Using the KA Backend

For large datasets, the KernelAbstractions backend parallelizes the computation:

using KernelAbstractions
grid_l2(config, grid_spec, time_spec;
        outfile = "output.nc",
        backend = CPU())

The KA backend uses sum-based accumulation (instead of Welford's incremental mean), which is fully parallelizable. The mean is computed at the end: mean = sum / weight.

Circular Footprints (GOSAT)

Some products, such as GOSAT SIF, describe a circular footprint with four bounding coordinates or with center coordinates plus a radius. Use CircularFootprintGridding for those products:

gosat_config = load_config("examples/gosat_sif_center_radius.toml")
grid(gosat_config, grid_spec, time_spec,
     CircularFootprintGridding(n_oversample=80);
     backend=resolve_backend("cpu"),
     outfile="gosat_sif.nc")

For center-plus-radius configs, set [basic] lat, [basic] lon, and either [basic] radius or scalar [circle] radius.

Step 5: Inspect the Output

The output NetCDF file contains:

  • lat, lon: Grid cell center coordinates
  • time: Time step dates
  • n: Number of observations per cell
  • One variable per entry in [grid] (e.g., sif_743, sif_735)
  • Optional _std suffix variables if compute_std=true
using NCDatasets
ds = Dataset("tropomi_sif_july2019.nc")
sif = ds["sif_743"][:, :, 1]   # First time step
weights = ds["n"][:, :, 1]
close(ds)

Center-Coordinate Gridding (MODIS)

For instruments without footprint bounds (e.g., MODIS), use grid_center:

config = load_config("examples/modis_reflectance.toml")
grid_spec = GridSpec(dlat=0.05f0, dlon=0.05f0)
time_spec = TimeSpec(DateTime("2019-01-01"), DateTime("2019-12-31"),
                     Dates.Day(1))
grid_center(config, grid_spec, time_spec;
            geo_provider=:modis,
            veg_indices=true,
            outfile="modis_2019.nc")

Cubed-Sphere (GEOS/GCHP) Output Grid

To compare TROPOMI directly against GEOS-Chem (GCHP) output on its native grid, swap the rectangular GridSpec for a CubedSphereGridSpec — everything else (the config, the grid call) is identical. The default geometry is the GEOS/GCHP convention (equal-distance gnomonic, GEOS-native panels, -10° shift), so Nc = 360 gives a GCHP-style C360 grid. Here we build a small C24 grid:

cs_spec = CubedSphereGridSpec(; Nc = 24)   # GMAOCubedSphereDefinition() by default
CubedSphereGridSpec{Float32, CubedSphereDefinition{GMAOEqualDistanceGnomonic, FourCornerNormalizedCenter, GEOSNativePanelConvention}}(24, CubedSphereDefinition{GMAOEqualDistanceGnomonic, FourCornerNormalizedCenter, GEOSNativePanelConvention}(GMAOEqualDistanceGnomonic(), FourCornerNormalizedCenter(), GEOSNativePanelConvention(), -10.0, :gmao_equal_distance), Float32[306.57373 306.57568 … 306.57568 306.57373; 309.80566 309.80765 … 309.80765 309.80566; … ; 30.194336 30.192371 … 30.192371 30.194336; 33.42628 33.424313 … 33.424313 33.42628;;; 36.573715 36.57568 … 36.57568 36.573715; 39.805668 39.807632 … 39.807632 39.805668; … ; 120.19433 120.19236 … 120.19236 120.19433; 123.42628 123.42431 … 123.42431 123.42628;;; 34.999996 31.755072 … 308.24493 305.0; 38.244923 34.999996 … 305.0 301.75507; … ; 121.755066 124.99999 … 215.00002 218.24493; 124.99999 128.24492 … 211.75508 215.00002;;; 126.573715 129.80566 … 210.19434 213.4263; 126.575676 129.80763 … 210.19238 213.42433; … ; 126.575676 129.80763 … 210.19238 213.42433; 126.573715 129.80566 … 210.19434 213.4263;;; 216.57373 219.80568 … 300.19434 303.4263; 216.5757 219.80765 … 300.19238 303.42432; … ; 216.5757 219.80765 … 300.19238 303.42432; 216.57373 219.80568 … 300.19434 303.4263;;; 215.00002 218.24493 … 301.75507 305.0; 211.75508 215.00002 … 305.0 308.24493; … ; 128.24492 124.99999 … 34.999996 31.755072; 124.99999 121.755066 … 38.244926 34.999996], Float32[-34.50544 -31.534674 … 31.534674 34.505436; -35.865456 -32.836414 … 32.836414 35.865456; … ; -35.865456 -32.836414 … 32.836414 35.865456; -34.505444 -31.534674 … 31.534674 34.50544;;; -34.50544 -31.534674 … 31.534674 34.505436; -35.865456 -32.836414 … 32.836414 35.865456; … ; -35.865456 -32.836414 … 32.836414 35.865456; -34.505444 -31.534674 … 31.534674 34.50544;;; 36.761852 38.245205 … 38.245205 36.761852; 38.245205 39.92845 … 39.92845 38.245205; … ; 38.245205 39.92845 … 39.92845 38.245205; 36.761852 38.245205 … 38.245205 36.761852;;; 34.50544 35.865456 … 35.865456 34.505444; 31.534678 32.836414 … 32.836414 31.534674; … ; -31.534678 -32.836414 … -32.836414 -31.534674; -34.505436 -35.865456 … -35.865456 -34.50544;;; 34.50544 35.865456 … 35.865456 34.505444; 31.534678 32.836414 … 32.836414 31.534674; … ; -31.534678 -32.836414 … -32.836414 -31.534674; -34.505436 -35.865456 … -35.865456 -34.50544;;; -36.761852 -38.245205 … -38.245205 -36.761852; -38.245213 -39.92845 … -39.92845 -38.245205; … ; -38.245205 -39.92845 … -39.92845 -38.245205; -36.761852 -38.245205 … -38.245205 -36.761852], Float32[305.0 305.0 … 305.0 305.0; 308.1741 308.17413 … 308.17413 308.1741; … ; 31.82589 31.825886 … 31.825886 31.82589; 34.999996 34.999996 … 34.999996 34.999996;;; 34.999996 34.999996 … 34.999996 34.999996; 38.17411 38.174114 … 38.174114 38.17411; … ; 121.82588 121.82587 … 121.82587 121.82588; 124.99999 124.99999 … 124.99999 124.99999;;; 34.999996 31.82589 … 308.1741 305.0; 38.17411 34.999996 … 305.0 301.8259; … ; 121.82588 124.99999 … 215.00002 218.17412; 124.99999 128.1741 … 211.8259 215.00002;;; 124.99999 128.1741 … 211.8259 215.00002; 124.99999 128.1741 … 211.8259 215.00002; … ; 124.99999 128.1741 … 211.8259 215.00002; 124.99999 128.1741 … 211.8259 215.00002;;; 215.00002 218.17412 … 301.8259 305.0; 215.00002 218.17412 … 301.8259 305.0; … ; 215.00002 218.17412 … 301.8259 305.0; 215.00002 218.17412 … 301.8259 305.0;;; 215.00002 218.17412 … 301.8259 305.0; 211.8259 215.00002 … 305.0 308.1741; … ; 128.1741 124.99999 … 34.999996 31.82589; 124.99999 121.82588 … 38.17411 34.999996], Float32[-35.264393 -32.325695 … 32.325695 35.264393; -36.69255 -33.6981 … 33.6981 36.69255; … ; -36.692554 -33.6981 … 33.6981 36.692554; -35.264393 -32.325695 … 32.325695 35.264393;;; -35.264393 -32.325695 … 32.325695 35.264393; -36.69255 -33.6981 … 33.6981 36.69255; … ; -36.692554 -33.6981 … 33.6981 36.692554; -35.264393 -32.325695 … 32.325695 35.264393;;; 35.264393 36.692554 … 36.69255 35.264393; 36.69255 38.313572 … 38.313572 36.692554; … ; 36.692554 38.313572 … 38.313572 36.69255; 35.264393 36.69255 … 36.692554 35.264393;;; 35.264393 36.69255 … 36.692554 35.264393; 32.325695 33.6981 … 33.6981 32.325695; … ; -32.325695 -33.6981 … -33.6981 -32.325695; -35.264393 -36.69255 … -36.692554 -35.264393;;; 35.264393 36.69255 … 36.692554 35.264393; 32.325695 33.6981 … 33.6981 32.325695; … ; -32.325695 -33.6981 … -33.6981 -32.325695; -35.264393 -36.69255 … -36.692554 -35.264393;;; -35.264393 -36.69255 … -36.692554 -35.264393; -36.692554 -38.313572 … -38.313572 -36.69255; … ; -36.69255 -38.313572 … -38.313572 -36.692554; -35.264393 -36.692554 … -36.69255 -35.264393])

The cube has Nc cells per panel edge across 6 panels. Internally it is folded into a 2D (Nc, 6·Nc) accumulator and un-folded to (Xdim, Ydim, nf) on write:

println("C$(cs_spec.Nc): grid_shape = $(grid_shape(cs_spec)), centers = $(size(cs_spec.centers_lon))")
C24: grid_shape = (24, 144), centers = (24, 24, 6)

Grid exactly as before — only the grid spec changed (sequential CPU path):

grid(config, cs_spec, time_spec, SubpixelGridding();
     outfile = "tropomi_sif_c24.nc")

Or from the CLI:

julia --project=. bin/grid.jl l2 --config examples/tropomi_sif.toml \
    --gridType cs --Nc 360 --startDate 2020-07-01 --stopDate 2020-07-16 \
    --dDays 16 -o tropomi_sif_c360.nc

The output matches native GCHP NetCDF (dims nf/Ydim/Xdim, 2D lons/lats and corner_lons/corner_lats, coordinates = "lons lats"), so it opens directly in Panoply and xarray. Plot one panel on its 2D coordinates:

import xarray as xr
ds = xr.open_dataset("tropomi_sif_c360.nc")
ds["sif_743"].isel(time=0, nf=2).plot(x="lons", y="lats")

See the Cubed-Sphere Grids page for the full convention, output layout, and the bin/run_tropomi_c360.sh fan-out runner.