geocif

MIKES EDITS

Models to visualize and forecast crop conditions and yields

Generate Climatic Impact-Drivers (CIDs) from Earth Observation (EO) data, build ML yield forecasting models, and produce agmet condition monitoring plots.

Climatic Impact-Drivers for Crop Yield Assessment at NASA Harvest

• Free software: MIT license
• Documentation: https://ritviksahajpal.github.io/yield_forecasting/

Setup

Requirements

• pixi

Pixi provisions Python 3.11 and the compiled geospatial stack (GDAL, rasterio, fiona, shapely, pyproj, rtree, cartopy) from conda-forge, so no system GDAL and no pre-built Windows wheels are needed.

Install

cd geocif                   # project root (where pyproject.toml lives)
pixi install                # resolves pixi.lock into .pixi/

Verify:

pixi run python -c "from osgeo import gdal; import geocif; print(geocif.__version__)"

Run commands inside the environment with pixi run <cmd>, or open a shell in it with pixi shell.

Extra environments

pixi run -e dev pytest      # dev tooling
pixi install -e full        # all pixi features

Fresh reinstall

rm -rf .pixi && pixi install

Installing as a library

pip install geocif still works and resolves from [project.dependencies]; this is how the cluster env pulls the package. Note that four optional extras (shap_fast, powershap, geospann, desreg) pin dependencies that conflict with geocif's own core pins and cannot currently be installed alongside it — see the comments in pyproject.toml.

Config files

Usage

Order matters: Config files are loaded left-to-right. When the same key appears in multiple files, the last file wins. The tool-specific file (geoextract.txt or geocif.txt) must be last so its [DEFAULT] values (countries, method, etc.) override the shared defaults in countries.txt.

config_dir = "/path/to/config"  # full path to your config directory

cfg_geoprepare = [f"{config_dir}/geobase.txt", f"{config_dir}/countries.txt", f"{config_dir}/crops.txt", f"{config_dir}/geoextract.txt"]
cfg_geocif = [f"{config_dir}/geobase.txt", f"{config_dir}/countries.txt", f"{config_dir}/crops.txt", f"{config_dir}/geocif.txt"]

geoprepare (download, extract, merge)

from geoprepare import geodownload
geodownload.run([f"{config_dir}/geobase.txt"])

from geoprepare import geoextract
geoextract.run(cfg_geoprepare)

from geoprepare import geomerge
geomerge.run(cfg_geoprepare)

geocif (indices, ML, agmet, analysis, experiments)

from geocif import indices_runner
indices_runner.run(cfg_geocif)

from geocif import geocif_runner
geocif_runner.run(cfg_geocif)

from geocif.agmet import geoagmet
geoagmet.run(cfg_geocif)

from geocif import season_monitor
season_monitor.run(cfg_geocif)  # pixel-level growing-season onset monitor (see below)

from geocif import calendar_validator
calendar_validator.run(cfg_cropcal)  # do NDVI transitions match the crop calendar?

from geocif import analysis
analysis.run(cfg_geocif)

from geocif import experiments
experiments.run(cfg_geocif, n_trials=30)

from geocif import yield_outlook
yield_outlook.run(cfg_geocif)  # uses config defaults (10 years, mean)
# yield_outlook.run(cfg_geocif, current_year=2026, n_years=10, aggregation="median")
# Set [ML] report_lite = True (with make_maps = True) to also emit one
# lightweight per-country PDF per run: cover + clickable TOC + "About GEOCIF" +
# per-crop accuracy (rRMSEp scorecard, scatter), predicted-yield table, and
# predicted/outlook maps for the best model by rRMSEp, plus references.

Cone-of-uncertainty figures

From an outlook DB produced with [ML] run_time_steps = all and estimate_ci_for_all = True (tabpfn native quantiles at every forecast issue month for every hindcast year), geocif.viz.cone renders per crop: a per-year hindcast validation grid of forecast cones vs observed, a full-size live-year cone with an empirical hindcast-error band, per-stage coverage vs the nominal level, and per-state mini-cones. Outputs land under {dir_output}/{project}/ml/analysis/{Month_DD_YYYY_HHhmm}/cone/, split into plots/ (PNG) and csvs/, with a lookup_plots_csvs.csv manifest in both -- the same shape the outlook products use.

Configure it with a [cone] section and the whole run is one call — which is how it runs on the cluster:

[cone]
start_year = 2005          ; first harvest year drawn in the grid
crops = ['maize', 'soybean']
; explicit per-crop DBs: an outlook DB written for one crop can also hold
; another crop's table built with the wrong detrend method
dbs = {'maize': 'outlook_09_06_2026_11h31.db', 'soybean': 'outlook_09_06_2026_11h21.db'}
model = tabpfn
units = bu/ac              ; or tn/ha
ncols = 5                  ; columns in the per-year grid
states_years = [2012, 2025, 2026]   ; years to draw per-state cone grids for
min_region_frac = 1.0      ; a stage must cover the FULL region pool to be drawn
use_nass_reference = True  ; USDA finals + monthly in-season track
nass_offline = False       ; True = cached NASS extract only, never call the API

Set [ML] make_cone_plots = True to render the cones automatically at the end of a yield_outlook run, into <analysis>/cone beside <analysis>/outlook, using the DB that run just wrote. It is off by default: the cones are a publication product rather than part of forecasting, and the NASS reference calls an external API a compute node may not reach. A failure there is logged and never discards a completed run. Otherwise render them standalone:

from geocif.viz import cone
cone.run([geobase, countries, crops, geocif])       # everything from [cone]
# CLI equivalent:
#   python -m geocif.viz.cone --config geobase.txt --config ... [--start-year 2005]
# or fully explicit, without a config:
#   python -m geocif.viz.cone --db X.db --table T --crop maize --out DIR [--no-nass]

USDA NASS reference. With use_nass_reference = True (key from [NASS] api_key), geocif.viz.nass pulls QuickStats state yields: the final estimates fill missing observed values, and the monthly in-season forecasts are drawn as USDA's own track next to the cone, so each forecast issue date is compared against the USDA number published that month. Two API facts the module encodes: freq_desc = MONTHLY is rejected outright — the monthly forecasts are reference_period_desc = "YEAR - AUG FORECAST" (…SEP/OCT/NOV) — and a "YEAR" row is the final estimate only once it was loaded after its crop year, because mid-season NASS mirrors the latest forecast into that row. The response is cached beside the figures so a rerun works without network access.

August is the earliest USDA yield forecast for corn and soybeans — there is no July or earlier estimate in the archive — so cone points issued in June and July have no USDA counterpart by construction, and the live season carries at most one square until the September report (~Sep 11). Use states_years to draw per-state grids for finished seasons as well, where the full Aug/Sep/Oct track is available.

No prose on the canvas. Figures carry a title, axis labels and a legend — nothing else. Interval level, coverage, sample sizes and units provenance live in the companion CSVs (and the run log), so they can be quoted in the surrounding text rather than shrinking the plotting area.

Crop-calendar validation

geocif.calendar_validator asks one question per (country, crop, season, calendar region): do the phenological transitions implied by satellite NDVI agree with the dates the GEOGLAM Crop Monitor calendar asserts? It compares mid-greenup (the end of calendar stage 1) against the steepest rise of a Fourier-fitted NDVI climatology, and mid-greendown (the end of stage 2) against the steepest fall, and reports the signed difference in days.

It runs at crop-calendar-region scale — the 574 polygons the calendar is actually authored against, not admin 1 — fed by a geoextract run that treats each calendar region as an admin_1 unit. Two one-off preparation steps come first: data_prep/prepare_calendar_regions.py makes the shipped region shapefile usable (it has no ID column, and geoextract drops a country silently without one), and data_prep/make_cropcal_config.py generates the 145-country config.

Alongside the rule-based method, the same targets are predicted by catboost, cubist, tabpfn and tabicl through the shared ml.trainers.auto_train factory, under three cross-validation schemes reported side by side: shuffled K-fold (leaky, the optimistic reference), leave-one- country-out, and grouped K-fold over 10 deg centroid tiles. Because the target is the calendar date, the rule-based delta and the model prediction error are the same quantity, so every method sits in one table. Metrics are also split by CM_Group: the original method was only ever exercised on AMIS countries, and 327 of the 575 regions here are EW-only.

This is a port of GEOGLAM/Code/Code/CropCalendar/. It reproduces the original numerics by default; geocif/cropcal/DEVIATIONS.md lists every intentional departure, the reason, and the config knob that restores the old behaviour.

Season monitor

geocif.season_monitor answers one question about the season that is starting right now, per pixel rather than per region: has onset happened, is it early or late against the 1981-2025 record, and where it has not happened, what the record says the chance is that it still arrives in the next 14 / 28 days. It reads geoprepare's daily intermediate rasters directly (CHIRPS rainfall, etref, CHIRTS tmax/tmin on the global 0.05 deg grid) and rasterizes the GEOGLAM Crop Monitor calendar per pixel, so it never touches the region-mean EO CSVs. It is a standalone monitoring product: it writes no CID rows, no ML features and nothing the yield pipeline reads.

Each pixel gets one of six states — BEFORE_WINDOW (the search window has not opened yet), NOT_STARTED (open, no candidate), FALSE_START (every candidate invalidated, none pending), PROVISIONAL (the latest candidate is still inside its validation window), CONFIRMED (a candidate survived, onset established) and NO_ONSET (season end plus validation passed with nothing confirmed) — plus the onset anomaly, the days already past the climatological median, the conditional onset probabilities, trailing 10/30-day rain and its percentile, the current dry run, and how many candidates have already failed.

Onset is 20 mm of rain in 3 days that is not followed by a 10-day dry spell (days below 1 mm) anywhere in the next 30 days; a candidate whose 30-day look-ahead is not fully observed yet stays PROVISIONAL rather than being confirmed early. Cessation runs a real daily soil-water bucket from onset onwards — S = min(max(S + P - PET, 0), soil_whc), capped at the water holding capacity and floored at zero — and fires on the first run of 5 consecutive empty days, never earlier than min_season_days after onset and never later than the next season's planting start.

Two settings move in opposite directions on the time axis and are easy to confuse. search_start_days_before_planting (default 30) opens the observed scan 30 days before the calendar planting start, so rain that arrives early is not censored. The 16-day CHIRPS-GEFS forecast extends the same series forward past the as-of day. One looks back at the start of the window, the other looks ahead at its end; a forecast trigger is reported next to the state and never changes it.

Days are stored as days since that pixel's own calendar planting start, not day-of-year. Negative values are therefore normal (the search opens before planting), and a cross-year season — Zimbabwe maize planted in November, harvested in May — needs no wrap correction anywhere.

PET. Cessation uses reference ET where it exists, otherwise Hargreaves from CHIRTS multiplied by a per-pixel monthly calibration factor (cached beside the climatology as pet_calibration.tif). Raw Hargreaves measured -0.5 to -1.9 mm/day biased against etref over Kenya, with near-zero spatial correlation; the out-of-sample monthly factor takes that to -0.17 mm/day and a correlation of 0.66. This is not academic: etref has been stalled at source since 2026-07-07, so the current season runs entirely on the calibrated fallback. Onset is unaffected — it is rain only.

Inputs are refreshed at run start. refresh_datasets (default CHIRPS and CHIRPS-GEFS) is handed to geoprepare's own dataset modules before anything is scored, so a monitor run is one call rather than a download step plus a run. The refresh skips itself when another geodownload.run is already in flight or the refresh lock file is under 6 h old, and it never raises: a failure is logged and the monitor proceeds on whatever is on disk. Set refresh_datasets = [] to skip it entirely.

Configure under [SEASON_MONITOR] in geocif.txt. Every key is optional — the defaults below are what you get with no section at all:

[SEASON_MONITOR]
precip_threshold = 20.0        ; mm in the accumulation window
window_days = 3                ; length of that window
dry_spell_days = 10            ; consecutive dry days that invalidate a candidate
dry_day_threshold = 1.0        ; mm; a day below this is dry
validation_days = 30           ; look-ahead a candidate must survive
search_start_days_before_planting = 30   ; observed scan opens this far BEFORE planting
soil_whc = 100.0               ; mm root-zone water holding capacity (the bucket cap)
min_season_days = 60           ; earliest cessation after onset; per-crop: min_season_days_maize = 90
empty_persist_days = 5         ; consecutive empty-bucket days that call cessation
cessation_grace_days = 60      ; scan past harvest; also capped at the next planting start
forecast_days = 16             ; CHIRPS-GEFS horizon, extending the series FORWARD
use_forecast = True
refresh_datasets = ['CHIRPS', 'CHIRPS-GEFS']   ; [] = no refresh
bbox_buffer = 0.5              ; degrees added around the country bounds
mask_to_cropland = True        ; computed on every land pixel, masked at display time
climatology_start_year = 1981
climatology_end_year = 2025
min_valid_years = 20           ; years needed before a climatology statistic is kept
rebuild_climatology = False    ; True = rebuild even when the params hash matches
n_workers = 8                  ; process / thread pool size
countries = ['kenya']          ; unset = [DEFAULT] countries
seasons_kenya = [1, 2]         ; per-country; unset = [country] seasons

The per-zone start-of-season report is a stage of the same run. For each crop-calendar zone it answers, over a window centred on the calendar planting date: how often has the season started historically, when, and how does the one now running compare. It is a pure aggregation of the cached per-year onset rasters, so it adds no computation, and it writes a zones/ subtree with a six-panel chart of every year against the window, a bar chart of how often each zone starts on time, and three tables.

One rule matters more than the rest: keep zone_report_window_days well below search_start_days_before_planting. When the two are equal, an onset that arrived before the search opened is reported AT the search boundary and cannot be told apart from one that happened exactly there, so the in-window share becomes a floor rather than a measurement. Measured on Kenya maize season 1 at a 30-day lead with a 30-day window: 82 % of West and 75 % of Rift Valley cropland was censored that way, while the four eastern zones were at 0 %. The shipped default is a 60-day lead against a 30-day window, leaving 30 days of margin. share_censored appears in every table and the search boundary is drawn on the chart, so the case is visible whenever it arises.

seasons_<country> exists because [country] seasons is shared with the CID, extract and merge pipelines, which need a merged CSV per season; adding a secondary season there would send them looking for files that do not exist. The monitor reads rasters and the crop calendar directly, so it needs no such file. Kenya season 2 is the October-February short rains.

from geocif import season_monitor
season_monitor.run(cfg_geocif)                       # everything from [SEASON_MONITOR]
season_monitor.run(cfg_geocif, as_of="2026-09-14")   # any config key, plus as_of, overrides the file
# CLI equivalent:
#   python -m geocif.season_monitor --config geobase.txt --config countries.txt \
#     --config crops.txt --config geocif.txt [--as-of YYYY-MM-DD] [--country kenya] \
#     [--no-refresh] [--no-forecast] [--rebuild-climatology] [--workers 8]

The as-of date is the last CHIRPS day actually on disk unless --as-of overrides it, and every raster tag, map title and table column states it. Each active (country, crop, season) runs in its own try/except, so one bad combination is logged and the rest still publish.

{dir_output}/{project}/ml/analysis/{Month_DD_YYYY}/season_monitor/{country}/{crop}/s{season}_hy{harvest_year}/
├── rasters/   one GeoTIFF per monitor layer, {layer}_{country}_{crop}_s{season}_hy{YEAR}_asof{YYYYMMDD}.tif
├── maps/      season_state, onset_anomaly_days, days_past_median, p_onset_28d,
│              rain_30d_percentile, fcst_trigger_days (+ lookup_plots_csvs.csv)
├── csvs/      one companion CSV per map, plus status_area_... (cropland-area share per
│              state, by Admin 1) and onset_summary_... (by Admin 1 and by calendar zone)
└── climatology/{maps,csvs}/   the 1981-2025 reference set: onset_median, onset_p75_minus_p25,
                               eos_median, lgs_median, false_start_rate, onset_n_valid

{dir_output}/{project}/phenology/climatology/{country}/{crop}/s{season}/   cache, NOT date-stamped

The harvest year is in the directory and in every file name because two of them can be active on the same day: a season stays open until harvest plus the validation window while the next one's search has already begun, which overlaps for any calendar spanning 305 days or more (DRC maize spans 472, Mexico maize 639). Without the year the later season silently overwrote the earlier one.

The cache is keyed by a params hash over the algorithm parameters, the year list and the grid, so a changed setting rebuilds it and an unchanged one is a no-op — there is nothing to invalidate by hand. A cache is only reused when it also holds every year it was asked for, so years lost to a transient outage are retried rather than frozen out, and it is refused outright if it was built with different parameters (an anomaly measured against the wrong reference is a systematic bias with nothing on the map to reveal it).

Status: nothing has run on the cluster yet. The numbers here come from an offline prototype comparison against the existing region-mean path. Aggregated back up to Admin 1, per-pixel onset matched the region-mean series to within a median of 1.6 days for Kenya maize 2024 (89 % of regions within 7 days, correlation 0.93). The gain shows up where the region mean flattens: for three Zimbabwe provinces the region series reports a flat "no onset" while the pixel field reports 0 %, 21 % and 21 % of cropland already started.

No prose on the canvas. Each map title carries the country, the crop, the layer, the admin level of the outlines and the data-through date — nothing else. Pixel distributions, per-state area shares, sample sizes and units live in the companion CSVs.

Cropmask optimizers

Two consumers of geoprepare extraction outputs that tune the cropland mask used downstream. Run after the corresponding geoprepare extractor has written its outputs.

# Uniform threshold T over the region (single absolute or rank-based knob).
# Reads geoprepare.extract_sweep output:
#   ${PATHS:dir_output}/threshold_sweep/{country}/{crop}/{country}_{crop}_s{season}_sweep.csv
from geocif import threshold_optimizer
threshold_optimizer.run(cfg_geocif)

# Per-cell binary mask — independent in/out decision per cropland cell.
# Reads geoprepare.extract_cells output:
#   ${PATHS:dir_output}/cell_optimizer/{country}/{crop}/{country}_{crop}_s{season}_cells.parquet
# Writes a production-mask parquet at the same location that geoextract picks up.
from geocif import cell_optimizer
cell_optimizer.run(cfg_geocif)

Configure under [THRESHOLD_OPTIMIZER] and [CELL_OPTIMIZER] in geocif.txt. Outputs land under ${PATHS:dir_output}/ml/analysis/{date}/{threshold_sweep_summary|cell_optimizer}/.

Using the optimized cell mask in production extraction

geoprepare 0.6.273+ can apply the per-cell mask produced by cell_optimizer during EO extraction. Opt in per country (or in [DEFAULT]) in geoextract.txt:

[DEFAULT]
use_optimized_mask = True

When the flag is on, geoprepare.extract_EO reads ${PATHS:dir_output}/cell_optimizer/{country}/{crop}/{country}_{crop}_s{season}_optimized_mask.parquet for every configured (country, crop, season) and AND-s it with the existing floor/ceiling AFI mask. Cells the optimizer marked included=False are dropped from the per-region aggregate even if they pass the floor/ceiling rule. Multi-season countries get the union across seasons — a cell is kept if any season's optimizer selected it.

Pipeline order with the optimized mask:

geoprepare.extract_cells.run(cfg_geoprepare)   # writes per-cell parquets
geocif.cell_optimizer.run(cfg_geocif)          # writes optimized_mask.parquet
geoprepare.geoextract.run(cfg_geoprepare)      # reads optimized_mask.parquet

extract_EO aborts at startup with a missing-parquet list if use_optimized_mask = True for any country whose mask hasn't been produced yet — silent fallback to the floor/ceiling rule when the operator asked for the overlay would be a confusing footgun, so it doesn't.

Currently wired in process_aef, process_fldas, process_chirps_mfc, process_soilgrids (the static + monthly-forecast EO paths). The daily-EO path through geom_extract (NDVI, daily CHIRPS, ESI, etc.) is not yet wired — track via a future change in geoprepare.

Annual (leave-one-out) masks

Enable annual_mask = True under [CELL_OPTIMIZER] in geocif.txt to produce one mask per historical year instead of a single pooled mask. For each year Y, the GA trains on every OTHER year — year Y's yield never sees the cell selection — and that mask is written to a _y{year}_optimized_mask.parquet file alongside the pooled one. geoprepare.extract_EO prefers the year-specific file when extracting year Y (FLDAS / CHIRPS-MFC, which are per-year datasets) and falls back to the pooled file for forecast / current years. AEF and SoilGrids (static) always use the pooled mask.

This closes the overfitting failure mode where the pooled mask was selected with year Y's yield as part of the training data — visible in pre-0.4.747 runs as regions whose Pearson r between yield and NDVI flipped sign after selection (the GA found anti-correlated cells because R² is sign-blind).

Cost. Roughly (n_years + 1) × the pooled-only default per region. On a country with 25 yield years that's ~26× more GA runs; expect runtime to scale accordingly. Opt in only when the data span justifies it.

Off by default. Existing configs without annual_mask continue to write the single pooled parquet.

ML models

geocif supports the following model types (configured via models in [DEFAULT]):

The spatial models (pygrf, tabpfn_gsa, tabfm_gsa, tabicl_gsa, and geospaNN) read region centroids from the lat/lon feature columns and therefore require [ML] include_lat_lon_as_feature = True — geocif fails fast at setup if it is unset, because otherwise every region's fit raises inside the per-region catch and the run "succeeds" having stored zero predictions. Each new model needs its own [<model>] config section with ML_model = True. Hyperparameter overrides live in [ML]: pygrf_band_width / pygrf_local_weight / pygrf_n_estimators / pygrf_max_features / pygrf_kernel / pygrf_resampled (defaults: adaptive 15%-of-n bandwidth, floor 20; local weight = global Moran's I of y), tabpfn_gsa_K / tabpfn_gsa_s (grid cells — rounded to a perfect square — and distant-sampling rate), george_kernel (expsquared | matern32 | matern52) / george_jitter, and bnn_epochs / bnn_batch_size / bnn_lr / bnn_prior_sigma / bnn_kl_weight / bnn_warmup_epochs / bnn_n_mc / bnn_calibrate_sigma, and mitra_hf_model / mitra_device / mitra_n_estimators / mitra_max_samples_support / mitra_max_samples_query / mitra_max_features / mitra_precision (plus mitra_fine_tune_steps / mitra_lr / mitra_warmup_steps / mitra_weight_decay / mitra_val_frac for mitra_ft).

tabfm_gsa and tabicl_gsa share tabpfn_gsa_K / tabpfn_gsa_s with tabpfn_gsa so the arms grid identically, and each takes its own inner-ensemble width — tabfm_gsa_n_estimators and tabicl_gsa_n_estimators (both default 8). The two keys are deliberately distinct: one config can run all three GSA arms side by side and no arm inherits another's width. causilo takes causilo_n_estimators / causilo_device / causilo_use_kv_cache / causilo_retain_preprocessing; its random_state is not configurable — it is forced from geocif's fold seed so causilo permutes in step with every other model.

TabPFN checkpoint selection

tabpfn pins TabPFN-3 by default, rather than following whatever the installed tabpfn package defaults to. This matters because tabpfn 9.0.0 moved its own default to TabPFN-3.5, which would silently have changed every archived result on upgrade. Override per config:

[ML]
tabpfn_model_version = v3.5      ; v2 | v2.5 | v2.6 | v3 | v3.5 | v3.5-fast, or "default" to follow the installed library
tabpfn_n_estimators = auto       ; an int, or "auto" to use the count the checkpoint declares for itself

Both are optional and both reproduce pre-9.0.0 behaviour when unset. Note that TabPFN-3.5 ships from a separate gated HuggingFace repo from TabPFN-3 and needs its own one-time licence acceptance at https://ux.priorlabs.ai; the repo's own LICENSE file says nothing about this.

Causilo licensing

causilo is kept out of geocif's core dependencies on purpose. Its code is Apache-2.0, but its model weights ship under the Causilo License v1.0, which permits non-commercial research, testing and evaluation only — commercial or production use, and hosted/API/SaaS services whether paid or free, require a separate licence from Nums AI Inc. Institutional affiliation and grant funding do not by themselves make a use commercial, and publishing papers and results (including model outputs) is explicitly permitted; routing predictions into an operational product is not. Install deliberately via pip install geocif[causilo].

The bnn model reports its own predictive uncertainty: total sigma = sqrt(aleatoric² + epistemic²) from Monte-Carlo weight sampling, rescaled by a factor fitted on the held-out most recent training year (two-pass: fit on earlier years to estimate the scale, then refit on all years). kl_weight defaults to 0.05 — the paper's full-KL objective collapses the sigma head to a near-constant. Its intervals bypass the conformal wrapper and land in the same "lower CI" / "upper CI" columns.

mitra / mitra_ft run Mitra-v2 (Tao et al. 2026, arXiv:2609.04540), an in-context tabular foundation model pretrained only on synthetic data: a 12-layer 2D Transformer (~76.7 M parameters) that attends across both rows and columns, with Apache-2.0 weights on the HuggingFace Hub (autogluon/mitra-regressor-2). The backbone ships inside autogluon.tabular[mitra], installed via the geocif[mitra] extra; on the cluster it is side-installed with --no-deps and reached through PYTHONPATH so the production environment's pinned pandas/pyarrow are untouched.

Two things make geocif/ml/mitra.py necessary rather than a thin call into AutoGluon. AutoGluon's public MitraRegressor hardcodes a scalar MSE head (dim_output = 1), which cannot drive the v2 regression checkpoint's 1,000-bin cross-entropy head; and stock AutoGluon decodes that head with argmax over bins, whereas the released Mitra-v2 recipe uses the softmax-weighted mean over bin centres, which its authors identify as the only quality-critical difference. geocif implements both.

The upside of the binned head is that Mitra's output is a genuine predictive distribution over the target range, so mitra gets native prediction intervals from its own quantiles with no conformal wrapper, like bnn. Keep folds inside the pretraining envelope — 160–5,120 support rows and 1–50 features, i.e. feature_selection = gOMP_medium — since above 256 features the model falls back to a train-only truncated-SVD reduction.

Feature selection methods

Configured via feature_selection in [ML]:

none, SelectKBest, BorutaPy, Leshy, gOMP, RFECV, RFE, lasso, mrmr, SHAP, stabl, PowerShap, BorutaShap, Genetic, feature_engine, multi

The gOMP tiers cap how many features survive: gOMP_high (500), gOMP_medium (50), gOMP_low (15). A county-scale sweep (usa_admin2 maize, cubist, matched county-years) found gOMP_medium best on every metric and 2.2× faster than both gOMP_high and none — none does not save time, it moves the cost into the model. Feature selection is the dominant cost of a fold (~74% at gOMP_high), so this setting matters more than the model choice.

Selection results are cached on disk and shared across the models of a fold (see Performance below), since the selection depends only on the training data and method.

Performance

Two settings govern how a run uses the machine:

[DEFAULT]
do_parallel_ml = True
fraction_cpus = 0.15          ; concurrent fold tasks, as a fraction of cores

[ML]
cache_feature_selection = True  ; reuse selections across models/runs (default)
threads_per_worker = auto       ; auto | <int> | -1 (unlimited)

• Feature-selection cache — selection is model-independent, so the first model of a fold computes it and the rest reuse it, on disk (workers are separate processes). Measured at county scale: catboost 32 min/fold computing it, cubist 21 min/fold reusing it. Keys are content-addressed, so any change to the data, method, or selector code is a miss rather than a stale hit.
• Thread budget — each worker gets cores / workers threads instead of every worker's model grabbing every core (which measured a load average of 940 on a 128-core node). Note fraction_cpus sets concurrency, not total footprint: the budget deliberately fills the machine without oversubscribing it. Set threads_per_worker explicitly to leave headroom for other users.

Cluster analysis

Optional analysis that clusters regions by their CID profiles and identifies which CIDs discriminate each cluster. Works with or without yield data — falls back to a proxy CID (e.g., AUC_NDVI) when yield is unavailable. Enabled via [ML]:

run_cluster_analysis = True
cluster_analysis_proxy = AUC_NDVI   ; proxy CID when yield is unavailable
cluster_analysis_max_k = 8          ; maximum clusters for silhouette selection
cluster_analysis_top_n = 20         ; top N CIDs in discrimination heatmap
cluster_analysis_variance = 0.85    ; cumulative PCA variance to retain

Pipeline: PCA dimensionality reduction → Ward's hierarchical clustering (silhouette-selected k) → Kruskal-Wallis + Cohen's d for CID discrimination → mutual information for CID-target association. Outputs: cluster map (choropleth), dendrogram, PCA biplot, discrimination heatmap with significance stars, target boxplot, and per-CID maps for top discriminating indices.

Spatial neighbor features

Optional GraphSAGE-style preprocessing that computes yield-correlation-weighted averages of neighboring regions' features. Enabled via [ML]:

use_spatial_neighbors = True
spatial_neighbor_method = knn   ; knn or full
spatial_neighbor_k = 5          ; number of nearest neighbors

For each admin region, the neighbor graph is built from training data using haversine distances and Pearson yield correlations as edge weights. Neighbor-aggregated features are added as nbr_* columns and flow through standard feature selection.

Experiments

The experiments runner (geocif.experiments) provides 6 experiments for model selection, feature importance, and hyperparameter tuning:

Dependencies: Experiments 1 and 2 require experiment 0 first. Experiments 3–5 are independent.

Configure in geocif.txt:

[experiments]
run_experiments = ["model_comparison", "cid_ablation"]
comparison_models = ["catboost", "tabpfn", "tabicl"]
n_trials = 30
n_trials_cid_types = 30
n_trials_cid_indices = 60
max_cid_indices = 25

Run:

from geocif import experiments
experiments.run(cfg_geocif)

Standalone diagnostics (geocif.experiments package)

Separate from the sweep runner above, geocif/experiments/ holds read-only diagnostics that answer one question against data already on disk — no model training, no writes to the outlook DB. Each writes a date-stamped ml/analysis/<stamp>/explore/<name>/ tree of plots/, tables/ and a README.md carrying the verdict and the knobs that change it.

from geocif.experiments import irrigation_kansas
irrigation_kansas.run(cfg)                          # kansas maize, ESI + heat + rain candidates
irrigation_kansas.run(cfg, state="nebraska")
irrigation_kansas.run(cfg, stress_cids=["TX90p", "KDD"])

irrigation_kansas standardises yields within county (Theil-Sen detrended residual over that county's own residual sd) before correlating with irr_share. Raw yield levels would only re-measure the irrigated/dryland level gap (171.3 vs 69.8 bu/ac in Kansas maize); the open question is whether irrigation buffers year-to-year variability, which is what the anomaly isolates. It ranks candidate stress CIDs by the R² their irr_share x stress interaction adds over a no-interaction model, so [ML] irrigation_stress_cid can be chosen from evidence rather than assumed.

Experiments output

The experiments runner writes to a dedicated DB and analysis folder under dir_output:

{dir_output}/
└── ml/
    ├── db/
    │   └── experiments_{MMMM_DD_YYYY_HH}H.db
    │
    └── analysis/
        └── {MMMM_DD_YYYY}/
            ├── experiments/                            # Experiment 0 (model comparison)
            │   ├── experiment_metrics.csv
            │   ├── heatmap_models.png
            │   ├── boxplot_models.png
            │   ├── regional_mape_models_{country}.png
            │   ├── error_distribution_models.png
            │   └── metric_comparison.png
            │
            └── optimization/                           # Optuna hyperparameter search
                ├── optuna_trials.csv
                ├── best_params.csv
                ├── convergence.png
                ├── optimization_history.png
                ├── param_importances.png
                └── parallel_coordinate.png

Outlook output

The yield outlook runner produces a diverging choropleth map showing current forecast yield as a percentage of the historical mean/median prediction per region, plus a combined CSV. Per (country, crop, model, stage) it also renders skill choropleths — mape_map, rmse_map, and r2_map (per-region MAPE, RMSE, and R²) — and metric-progression plots across forecast stages, each in both a full per-region version and a decluttered *_national version (single national line with a gray ±1 std-across-regions band).

In addition, model performance is aggregated to every admin level above the run level ([ML] plot_parent_aggregations, default True): an admin_2 (county) run also gets admin_1 (state) and national outputs; an admin_1 run gets national outputs. Observed/predicted yields are area-weighted means per (parent, year) — groups with missing/zero Area (ha) fall back to unweighted means — with the county→state mapping reused from the production-statistics file (ml.stats.admin1_lookup). Per level this writes a pooled obs-vs-pred scatter (+ hexbin), a per-parent metrics CSV (r², rRMSE%, MAPE% + pooled row), a national obs-vs-pred time series, per-parent rRMSE%/r² choropleths (admin_1 only; boundary from [ML] parent_boundary_admin_1, falling back to Level_1.shp — skipped with a warning when missing or when <50% of parent names match), and a lookup_plots_csvs.csv mapping each plot to its companion CSV.

{dir_output}/
└── ml/
    └── analysis/
        └── {MMMM_DD_YYYY}/
            └── outlook/
                ├── yield_outlook_{country}_{crop}_{model}_{stage}_{year}.png
                ├── yield_outlook_{year}.csv
                ├── plots/{model}/{country}/{crop}/{admin_1|national}/
                │   ├── scatter_{level}_{country}_{crop}_{model}.png (+ _hexbin)
                │   ├── timeseries_national_{country}_{crop}_{model}.png  (national only)
                │   ├── rrmse_map_admin_1_...png / r2_map_admin_1_...png  (admin_1 only)
                │   └── lookup_plots_csvs.csv
                └── csvs/{model}/{country}/{crop}/{admin_1|national}/
                    ├── aggregated_predictions_{level}_{country}_{crop}_{model}.csv
                    ├── metrics_{level}_{country}_{crop}_{model}.csv
                    ├── scatter/timeseries companion CSVs
                    └── lookup_plots_csvs.csv

Config file documentation

geobase.txt

Shared paths and dataset settings. All directory paths are derived from dir_base.

[PATHS]
dir_base = /gpfs/data1/cmongp1/GEO

dir_inputs = ${dir_base}/inputs
dir_logs = ${dir_base}/logs
dir_download = ${dir_inputs}/download
dir_intermed = ${dir_inputs}/intermed
dir_metadata = ${dir_inputs}/metadata
dir_condition = ${dir_inputs}/crop_condition
dir_crop_inputs = ${dir_condition}/crop_t20

dir_boundary_files = ${dir_metadata}/boundary_files
dir_crop_calendars = ${dir_metadata}/crop_calendars
dir_crop_masks = ${dir_metadata}/crop_masks
dir_images = ${dir_metadata}/images
dir_production_statistics = ${dir_metadata}/production_statistics

dir_output = ${dir_base}/outputs

[DATASETS]
datasets = ['CHIRPS', 'CPC', 'NDVI', 'ESI', 'NSIDC', 'AEF']

countries.txt

Single source of truth for per-country config. Shared by both geoprepare and geocif.

[DEFAULT]
boundary_file = gaul1_asap_v04.shp
admin_level = admin_1
seasons = [1]
crops = ['maize']
category = AMIS
use_cropland_mask = False
calendar_file = crop_calendar.csv

; AMIS countries (inherit from DEFAULT, override crops if needed)
[argentina]
crops = ['soybean', 'winter_wheat', 'maize']

; EWCM countries (full per-country config)
[kenya]
category = EWCM
admin_level = admin_1
seasons = [1, 2]
use_cropland_mask = True
boundary_file = adm_shapefile.gpkg
calendar_file = EWCM_2025-04-21.xlsx
crops = ['maize']

[malawi]
category = EWCM
admin_level = admin_2
use_cropland_mask = True
boundary_file = adm_shapefile.gpkg
calendar_file = EWCM_2025-04-21.xlsx
crops = ['maize']

Restricting a run to specific regions (run_regions)

Optional per-country knob (falls back to [DEFAULT]). Unset = run every region, which is the default behaviour.

[usa]
admin_level = admin_2

; names below are admin_1 (state) even though the run is admin_2 (county)
run_regions_level = admin_1
; flat list  -> applies to every crop
run_regions = ["illinois", "iowa"]
; or per crop -> a crop missing from the dict runs all its regions
run_regions = {"maize": ["illinois", "iowa"], "soybean": ["arkansas"]}

• run_regions_level declares the admin level the names are given at; it defaults to the run's own admin_level.
• An admin_1 selection on an admin_2 run keeps every county inside the selected states (parent lookup via the production-statistics file, same normalization as the yield join). An admin_2 selection on an admin_1 run maps the county names up to their parent states.
• Names are matched case-insensitively with _ treated as a space, so South Dakota, south_dakota and SOUTH DAKOTA are the same region.
• Selected names that match nothing are listed in a warning; a selection that matches no region raises a ValueError instead of running empty.
• The filter is applied in _prepare_ml_dataframe, so training, LOOCV, the SQLite output, plots and parent aggregations all see only the selected regions.

crops.txt

Crop mask filenames and calendar category definitions.

; Crop masks
[maize]
mask = Percent_Maize.tif

[winter_wheat]
mask = Percent_Winter_Wheat.tif

[sorghum]
mask = cropland_v9.tif

; Calendar categories
[EWCM]
use_cropland_mask = True
calendar_file = EWCM_2026-01-05.xlsx
crops = ['maize', 'sorghum', 'millet', 'rice', 'winter_wheat', 'teff']
eo_model = ['aef', 'nsidc_surface', 'nsidc_rootzone', 'ndvi', 'cpc_tmax', 'cpc_tmin', 'chirps', 'chirps_gefs', 'esi_4wk']

[AMIS]
calendar_file = AMISCM_2026-01-05.xlsx

geoextract.txt

Extraction-only settings for geoprepare. Loaded last so its [DEFAULT] overrides shared defaults.

[DEFAULT]
method = JRC
redo = False
threshold = True
floor = 20
ceil = 90
countries = ["malawi"]
forecast_seasons = [2022]

[PROJECT]
parallel_extract = True
parallel_merge = False

geocif.txt

Indices, ML, and agmet settings for geocif. Country overrides go here when geocif needs different values than countries.txt (e.g., a subset of crops).

[AGMET]
eo_plot = ['ndvi', 'chirts_era5_tmax', 'chirts_era5_tmin', 'chirps', 'esi_4wk', 'nsidc_surface', 'nsidc_rootzone']
logo_harvest = harvest.png
logo_geoglam = geoglam.png

; Country overrides (only where geocif differs from countries.txt)
[ethiopia]
crops = ['winter_wheat']

[bangladesh]
crops = ['rice']
admin_level = admin_2
boundary_file = bangladesh.shp

; ML model definitions
[catboost]
ML_model = True

[analog]
ML_model = False

[ML]
model_type = REGRESSION
target = Yield (tn per ha)
feature_selection = gOMP
cluster_strategy = single
check_yield_trend = False
use_spatial_neighbors = True
spatial_neighbor_method = knn
spatial_neighbor_k = 5
lag_yield_as_feature = True
lag_years = 3
median_yield_as_feature = False
median_years = 5
include_lat_lon_as_feature = False
panel_model = True
cat_features = ["Harvest Year", "Region_ID", "Region"]
outlook_n_years = 10        ; Number of historical years for yield outlook comparison
outlook_aggregation = mean  ; mean or median
run_time_steps = latest         ; latest, current, all, N, pre_season, or auto
                                ; pre_season: forecast BEFORE the season starts from
                                ; S2S/FLDAS leads only — forecast_seasons may then
                                ; include a season that has not started yet
                                ; (e.g. south_africa maize 2027 launched Sep 2026)
run_cluster_analysis = False
cluster_analysis_proxy = AUC_NDVI
cluster_analysis_max_k = 8
cluster_analysis_top_n = 20
cluster_analysis_variance = 0.85

[LOGGING]
log_level = INFO

[DEFAULT]
data_source = harvest
method = monthly_r
project_name = geocif
countries = ["kenya"]
crops = ['maize']
admin_level = admin_1
models = ['catboost']
seasons = [1]
threshold = True
floor = 20

FLDAS forecast overlay

When FLDAS columns are present in the merged data (e.g. fldas_tair_tavg_lead0 through _lead5), agmet plots automatically overlay forecast dots on matching panels:

Each lead time (0–5) appears as a diamond marker with decreasing opacity (lead 0 = most opaque). Dots beyond the harvest date are suppressed. No config changes are needed — detection is automatic.

Release

To publish a new version to PyPI:

1. Bump __version__ in geocif/__init__.py and version in pyproject.toml
2. Build and upload:

   uv build
   uvx twine upload dist/geocif-<version>*
   

3. Commit:

   git add geocif/__init__.py pyproject.toml
   git commit -m "Bump to <version>"
   

Credits

This project was supported by NASA Applied Sciences Grant No. 80NSSC17K0625 through the NASA Harvest Consortium, and the NASA Acres Consortium under NASA Grant #80NSSC23M0034.