Metadata-Version: 2.4
Name: geocif
Version: 0.4.1040
Summary: Models to visualize and forecast crop conditions and yields
Author-email: Ritvik Sahajpal <ritvik@umd.edu>
License: MIT
Project-URL: Homepage, https://ritviksahajpal.github.io/yield_forecasting/
Keywords: geocif
Classifier: Development Status :: 2 - Pre-Alpha
Classifier: Intended Audience :: Developers
Classifier: License :: OSI Approved :: MIT License
Classifier: Natural Language :: English
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: boruta>=0.4.3
Requires-Dist: catboost>=1.2.8
Requires-Dist: fiona
Requires-Dist: gdal>=3.10
Requires-Dist: pyproj
Requires-Dist: rasterio
Requires-Dist: rtree
Requires-Dist: shap>=0.48.0
Requires-Dist: numba>=0.59
Requires-Dist: shapiq>=1.0
Requires-Dist: shapely
Requires-Dist: optuna
Requires-Dist: xarray>=2026.2.0
Requires-Dist: pooch>=1.8.0
Requires-Dist: arrow>=1.4.0
Requires-Dist: icclim>=7.0.4
Requires-Dist: geoprepare>=0.6.275
Requires-Dist: logzero>=1.7.0
Requires-Dist: geopandas>=1.1.2
Requires-Dist: tabpfn>=6.4.1
Requires-Dist: tabicl>=2.0.2
Requires-Dist: statsmodels>=0.14.6
Requires-Dist: palettable>=3.3.3
Requires-Dist: seaborn>=0.13.2
Requires-Dist: scikit-misc>=0.5.2
Requires-Dist: setuptools<81
Requires-Dist: choix>=0.3.4
Requires-Dist: scienceplots>=2.0.0
Requires-Dist: cartopy>=0.22
Requires-Dist: Rbeast>=0.1.20
Requires-Dist: scikit-learn>=1.4
Requires-Dist: bottleneck>=1.3
Requires-Dist: mapclassify>=2.5
Requires-Dist: pymannkendall>=1.4
Requires-Dist: pangres>=4.0
Requires-Dist: kneed>=0.8
Requires-Dist: lifelines>=0.27
Requires-Dist: esda>=2.5
Requires-Dist: libpysal>=4.10
Requires-Dist: crepes>=0.6
Requires-Dist: cubist>=1.0
Requires-Dist: mapie>=0.8
Requires-Dist: merf>=1.0
Requires-Dist: ngboost>=0.5
Requires-Dist: PyGRF>=0.0.12
Requires-Dist: tabpfn-extensions>=0.4
Requires-Dist: BorutaShap>=1.0
Requires-Dist: arfs>=2.0
Requires-Dist: feature-engine>=1.6
Requires-Dist: mrmr-selection>=0.2
Requires-Dist: sklearn-genetic-opt>=0.10
Requires-Dist: stabl>=0.0.1
Requires-Dist: cachetools>=5.0
Requires-Dist: geopy>=2.0
Requires-Dist: scikit-image>=0.21
Requires-Dist: matplotlib<3.11
Requires-Dist: pandas<3
Requires-Dist: openpyxl
Provides-Extra: dashboard
Requires-Dist: panel>=1.4.0; extra == "dashboard"
Requires-Dist: hvplot>=0.10.0; extra == "dashboard"
Requires-Dist: holoviews>=1.18; extra == "dashboard"
Provides-Extra: aquacrop
Requires-Dist: aquacrop>=3.0; extra == "aquacrop"
Provides-Extra: aquacrop-calibration
Requires-Dist: aquacrop>=3.0; extra == "aquacrop-calibration"
Requires-Dist: pygmo>=2.19; extra == "aquacrop-calibration"
Requires-Dist: SALib>=1.5; extra == "aquacrop-calibration"
Provides-Extra: spatial
Requires-Dist: pysal>=2.6; extra == "spatial"
Provides-Extra: causilo
Requires-Dist: causilo>=1.0.1; extra == "causilo"
Provides-Extra: gee
Requires-Dist: earthengine-api>=1.0; extra == "gee"
Requires-Dist: geemap>=0.30; extra == "gee"
Provides-Extra: narrative
Requires-Dist: pymupdf>=1.23; extra == "narrative"
Requires-Dist: pdfplumber>=0.10; extra == "narrative"
Requires-Dist: reportlab>=4.0; extra == "narrative"
Requires-Dist: anthropic>=0.30; extra == "narrative"
Provides-Extra: shap-fast
Requires-Dist: fasttreeshap>=0.1; extra == "shap-fast"
Provides-Extra: powershap
Requires-Dist: powershap>=0.0.10; extra == "powershap"
Provides-Extra: mitra
Requires-Dist: autogluon.tabular[mitra]<1.7,>=1.6; extra == "mitra"
Provides-Extra: geospann
Requires-Dist: geospaNN>=0.1; extra == "geospann"
Provides-Extra: george
Requires-Dist: george>=0.4.4; extra == "george"
Provides-Extra: desreg
Requires-Dist: desReg>=0.1; extra == "desreg"
Dynamic: license-file

# geocif

MIKES EDITS

[![image](https://img.shields.io/pypi/v/geocif.svg)](https://pypi.python.org/pypi/geocif)

**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](https://www.loom.com/share/5c2dc62356c6406193cd9d9725c2a6a9)

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


## Setup

### Requirements

- [pixi](https://pixi.sh/latest/#installation)

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

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

Verify:

```bash
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

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

### Fresh reinstall

```bash
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

| File | Purpose | Used by |
|------|---------|---------|
| [`geobase.txt`](#geobasetxt) | Paths, shapefile column mappings | both |
| [`countries.txt`](#countriestxt) | Per-country config (boundary files, admin levels, seasons, crops) | both |
| [`crops.txt`](#cropstxt) | Crop masks, calendar categories (EWCM, AMIS) | both |
| [`geoextract.txt`](#geoextracttxt) | Extraction-only settings (method, threshold, parallelism) | geoprepare |
| [`geocif.txt`](#geociftxt) | Indices/ML/agmet settings, country overrides, runtime selections | geocif |

## 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`.

```python
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)

```python
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)

```python
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:

```ini
[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:

```python
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:

```ini
[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.


```python
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.

```python
# 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`:

```ini
[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:**
```python
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]`):

| Model | Key | Type |
|-------|-----|------|
| CatBoost | `catboost` | Gradient boosting |
| XGBoost | `xgboost` | Gradient boosting |
| TabPFN | `tabpfn` | Prior-fitted network |
| TabICL | `tabicl` | In-context learning |
| NGBoost | `ngboost` | Natural gradient boosting |
| YDF | `ydf` | Yggdrasil decision forests |
| Oblique RF | `oblique` | Oblique random forest |
| Cubist | `cubist` | Rule-based regression |
| MERF | `merf` | Mixed effects random forest |
| Linear | `linear` | LassoCV / LogisticRegressionCV |
| GAM | `gam` | Generalized additive model |
| GeoSpaNN | `geospaNN` | Geospatial neural network |
| PyGRF | `pygrf` | Geographical random forest (global + per-sample local RFs, Moran's-I blend) |
| TabPFN-GSA | `tabpfn_gsa` | TabPFN with geospatial sparse attention (local in-context sampling; git-only dep) |
| george | `george` | Exact GP regression (dfm/george, fitted kernel + noise; `geocif[george]` extra) |
| BNN | `bnn` | Bayesian NN (Ma et al. 2021 RSE): variational two-headed torch net with native heteroscedastic prediction intervals |
| Mitra-v2 | `mitra` | In-context tabular foundation model (arXiv:2609.04540, Apache-2.0 weights); 1,000-bin regression head gives native prediction intervals; `geocif[mitra]` extra |
| Mitra-v2 fine-tuned | `mitra_ft` | Same checkpoint under the released 50-step fine-tune — minutes per fold, GPU recommended |
| TabFM-GSA | `tabfm_gsa` | GSA sampler with Google Research TabFM as the local in-context estimator; GPU strongly recommended |
| TabICL-GSA | `tabicl_gsa` | GSA sampler with TabICL as the local in-context estimator. Note it gets **conformal** intervals, not TabICL's native quantiles — GSA's function backend returns point predictions only |
| Causilo | `causilo` | Pretrained tabular foundation model (Nums AI); 999 native quantiles give native prediction intervals; `geocif[causilo]` extra. **Weights are non-commercial research only — see the licensing note below** |
| GPR | `gpr` | sklearn Gaussian process regression |
| Median | `median` | Median baseline |
| Analog | `analog` | Analogous year baseline |

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:

```ini
[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:

```ini
[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]`:

```ini
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]`:

```ini
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:

| # | Config name | Internal name | What it does |
|---|-------------|---------------|--------------|
| 0 | `model_comparison` | `models` | Runs each model in `comparison_models` head-to-head. Produces Bradley-Terry ranking, scatter plots, MAPE bars. Identifies best model per country (required by experiments 1 & 2). |
| 1 | `cid_ablation` | `cids` | Runs the best model once per CID Type in isolation (Cold alone, FLDAS alone, etc.). Shows which climate driver category contributes most. Produces MAPE-by-CID bar chart, region×CID heatmap, year×CID chart, CID rank over time. |
| 2 | `region_filter` | `region_filter` | Drops low-production regions and re-runs the best model to test if excluding noisy regions improves national accuracy. |
| 3 | `optuna` | `optuna` | Bayesian (TPE) search over ML hyperparameters (learning rate, depth, regularization, etc.). Produces convergence, parameter importance, and parallel coordinate plots. |
| 4 | `optuna_cid_types` | `optuna_cid_types` | Bayesian search for the best combination of CID Type categories (e.g. Rain+VI+ESI may beat using all 8 types). |
| 5 | `optuna_cid_indices` | `optuna_cid_indices` | Bayesian search for the best subset of individual CID indices (e.g. PRCPTOT + AUC_NDVI + TG90p). Capped at `max_cid_indices` per trial. |

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

Configure in `geocif.txt`:

```ini
[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:

```python
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.

| Module | Question it answers |
|---|---|
| `residuals_vs_cid` | Which CID anomalies does the yield model fail to absorb? |
| `drought_audit` | Do drought years actually look anomalous in the inputs, or did the pipeline flatten them? |
| `irrigation_kansas` | Does irrigated share interact with in-season stress to move county yield anomalies — and which stress CID and month carry it? |

```python
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`.

```ini
[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.

```ini
[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.

```ini
[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.

```ini
; 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.

```ini
[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).

```ini
[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:

| FLDAS variable | Target panel |
|---|---|
| `fldas_tair_tavg` | Temperature |
| `fldas_totalprecip_tavg` | Daily precipitation |
| `fldas_soilmoist_tavg` | Soil moisture (surface) |

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:
   ```bash
   uv build
   uvx twine upload dist/geocif-<version>*
   ```
3. Commit:
   ```bash
   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.
