feat: add functionality to skip empty tiles in TIFF processing and clean up orphaned files

Code improvements

.claude/settings.local.json

@@ -9,7 +9,8 @@
       "Bash(/c/Users/timon/AppData/Local/r-miniconda/python.exe -c \":*)",
       "Bash(python3 -c \":*)",
       "Bash(Rscript -e \":*)",
-      "Bash(\"/c/Program Files/R/R-4.4.3/bin/x64/Rscript.exe\" -e \":*)"
+      "Bash(\"/c/Program Files/R/R-4.4.3/bin/x64/Rscript.exe\" -e \":*)",
+      "Read(//c/Users/timon/Documents/SmartCane_code/r_app/**)"
     ]
   }
 }

create_field_checklist.R

@@ -0,0 +1,186 @@
+# Creates an Excel checklist from pivot.geojson
+# Fields sorted largest to smallest, split across Timon/Joey/Dimitra side-by-side
+
+# Install packages if needed
+if (!requireNamespace("jsonlite", quietly = TRUE)) install.packages("jsonlite", repos = "https://cloud.r-project.org")
+if (!requireNamespace("openxlsx", quietly = TRUE)) install.packages("openxlsx", repos = "https://cloud.r-project.org")
+
+library(jsonlite)
+library(openxlsx)
+
+# ---- Load GeoJSON ----
+geojson_path <- "laravel_app/storage/app/angata/pivot.geojson"
+gj <- fromJSON(geojson_path, simplifyVector = FALSE)
+features <- gj$features
+cat(sprintf("Total features: %d\n", length(features)))
+
+# ---- Shoelace area (degrees²) ----
+shoelace <- function(ring) {
+  n <- length(ring)
+  lons <- sapply(ring, `[[`, 1)
+  lats <- sapply(ring, `[[`, 2)
+  area <- 0
+  for (i in seq_len(n)) {
+    j <- (i %% n) + 1
+    area <- area + lons[i] * lats[j] - lons[j] * lats[i]
+  }
+  abs(area) / 2
+}
+
+# ---- Approx area in m² ----
+area_m2 <- function(ring) {
+  R <- 6371000
+  lats <- sapply(ring, `[[`, 2)
+  mean_lat <- mean(lats)
+  lat_rad  <- mean_lat * pi / 180
+  m_per_deg_lat <- R * pi / 180
+  m_per_deg_lon <- R * cos(lat_rad) * pi / 180
+  shoelace(ring) * m_per_deg_lat * m_per_deg_lon
+}
+
+# ---- Compute feature areas ----
+compute_area <- function(feat) {
+  geom <- feat$geometry
+  total <- 0
+  if (geom$type == "MultiPolygon") {
+    for (polygon in geom$coordinates) {
+      total <- total + area_m2(polygon[[1]])  # outer ring
+    }
+  } else if (geom$type == "Polygon") {
+    total <- total + area_m2(geom$coordinates[[1]])
+  }
+  total
+}
+
+field_names <- sapply(features, function(f) f$properties$field)
+areas_m2    <- sapply(features, compute_area)
+areas_ha    <- areas_m2 / 10000
+
+df <- data.frame(
+  field   = field_names,
+  area_ha = round(areas_ha, 2),
+  stringsAsFactors = FALSE
+)
+
+# Sort largest to smallest
+df <- df[order(df$area_ha, decreasing = TRUE), ]
+df$rank <- seq_len(nrow(df))
+
+cat("\nTop 10 fields by area:\n")
+print(head(df[, c("rank", "field", "area_ha")], 10))
+
+# ---- Split: Timon=1st, Joey=2nd, Dimitra=3rd ----
+idx <- seq_len(nrow(df))
+timon   <- df[idx %% 3 == 1, ]
+joey    <- df[idx %% 3 == 2, ]
+dimitra <- df[idx %% 3 == 0, ]
+
+cat(sprintf("\nSplit: Timon=%d, Joey=%d, Dimitra=%d\n",
+            nrow(timon), nrow(joey), nrow(dimitra)))
+
+# ---- Build Excel ----
+wb <- createWorkbook()
+addWorksheet(wb, "Field Checklist")
+
+# Header colors
+col_timon   <- "1F6AA5"
+col_joey    <- "2E7D32"
+col_dimitra <- "7B1FA2"
+alt_timon   <- "D6E4F0"
+alt_joey    <- "D7F0D8"
+alt_dimitra <- "EDD7F0"
+
+header_font <- createStyle(fontName = "Calibri", fontSize = 11, fontColour = "FFFFFF",
+                            halign = "CENTER", valign = "center", textDecoration = "bold",
+                            border = "TopBottomLeftRight")
+sub_font    <- createStyle(fontName = "Calibri", fontSize = 10, fontColour = "FFFFFF",
+                            halign = "CENTER", valign = "center", textDecoration = "bold",
+                            border = "TopBottomLeftRight")
+
+# Title row
+writeData(wb, "Field Checklist",
+          "Angata Pivot Field Checklist — sorted largest to smallest",
+          startRow = 1, startCol = 1)
+mergeCells(wb, "Field Checklist", cols = 1:14, rows = 1)
+addStyle(wb, "Field Checklist",
+         createStyle(fontName = "Calibri", fontSize = 13, textDecoration = "bold",
+                     halign = "CENTER", valign = "center",
+                     fgFill = "F0F0F0"),
+         rows = 1, cols = 1)
+setRowHeights(wb, "Field Checklist", rows = 1, heights = 28)
+
+# Person block writer
+write_person_block <- function(wb, ws_name, data, start_col, hdr_color, alt_color, person_name) {
+  end_col <- start_col + 3
+
+  # Person name header (row 2)
+  mergeCells(wb, ws_name, cols = start_col:end_col, rows = 2)
+  writeData(wb, ws_name, person_name, startRow = 2, startCol = start_col)
+  addStyle(wb, ws_name,
+           createStyle(fontName = "Calibri", fontSize = 12, fontColour = "FFFFFF",
+                       textDecoration = "bold", halign = "CENTER", valign = "center",
+                       fgFill = hdr_color, border = "TopBottomLeftRight",
+                       borderColour = "999999"),
+           rows = 2, cols = start_col:end_col)
+
+  # Sub-headers (row 3)
+  sub_headers <- c("#", "Field", "Area (ha)", "Checked \u2713")
+  writeData(wb, ws_name, as.data.frame(t(sub_headers)),
+            startRow = 3, startCol = start_col, colNames = FALSE)
+  addStyle(wb, ws_name,
+           createStyle(fontName = "Calibri", fontSize = 10, fontColour = "FFFFFF",
+                       textDecoration = "bold", halign = "CENTER", valign = "center",
+                       fgFill = hdr_color, border = "TopBottomLeftRight",
+                       borderColour = "999999"),
+           rows = 3, cols = start_col:end_col)
+
+  # Data rows (starting row 4)
+  n <- nrow(data)
+  for (i in seq_len(n)) {
+    row_num <- i + 3
+    bg <- if (i %% 2 == 0) alt_color else "FFFFFF"
+    # Rank
+    writeData(wb, ws_name, i, startRow = row_num, startCol = start_col)
+    # Field name
+    writeData(wb, ws_name, data$field[i], startRow = row_num, startCol = start_col + 1)
+    # Area
+    writeData(wb, ws_name, data$area_ha[i], startRow = row_num, startCol = start_col + 2)
+    # Checked (empty)
+    writeData(wb, ws_name, "", startRow = row_num, startCol = start_col + 3)
+
+    row_style <- createStyle(fontName = "Calibri", fontSize = 10, halign = "center",
+                              fgFill = bg, border = "TopBottomLeftRight",
+                              borderColour = "CCCCCC")
+    field_style <- createStyle(fontName = "Calibri", fontSize = 10, halign = "left",
+                                fgFill = bg, border = "TopBottomLeftRight",
+                                borderColour = "CCCCCC")
+    addStyle(wb, ws_name, row_style,   rows = row_num, cols = start_col)
+    addStyle(wb, ws_name, field_style, rows = row_num, cols = start_col + 1)
+    addStyle(wb, ws_name, row_style,   rows = row_num, cols = start_col + 2)
+    addStyle(wb, ws_name, row_style,   rows = row_num, cols = start_col + 3)
+  }
+}
+
+write_person_block(wb, "Field Checklist", timon,   1,  col_timon,   alt_timon,   "Timon")
+write_person_block(wb, "Field Checklist", joey,    6,  col_joey,    alt_joey,    "Joey")
+write_person_block(wb, "Field Checklist", dimitra, 11, col_dimitra, alt_dimitra, "Dimitra")
+
+# Column widths (col 5 and 10 = spacers)
+setColWidths(wb, "Field Checklist",
+             cols  = c(1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12, 13, 14),
+             widths = c(5, 14, 10, 12,  2,  5, 14, 10, 12,   2,  5, 14, 10, 12))
+
+# Row heights
+setRowHeights(wb, "Field Checklist", rows = 2:3, heights = c(22, 18))
+max_rows <- max(nrow(timon), nrow(joey), nrow(dimitra))
+setRowHeights(wb, "Field Checklist", rows = 4:(max_rows + 3), heights = 16)
+
+# Freeze panes below header
+freezePane(wb, "Field Checklist", firstActiveRow = 4)
+
+# Save
+out_path <- "angata_field_checklist.xlsx"
+saveWorkbook(wb, out_path, overwrite = TRUE)
+cat(sprintf("\nExcel saved to: %s\n", out_path))
+cat(sprintf("Total: %d fields — Timon: %d, Joey: %d, Dimitra: %d\n",
+            nrow(df), nrow(timon), nrow(joey), nrow(dimitra)))

create_field_checklist.py

@@ -0,0 +1,194 @@
+"""
+Creates an Excel checklist from pivot.geojson, with fields sorted largest to smallest,
+split across Timon (1st), Joey (2nd), Dimitra (3rd) in a single sheet side-by-side.
+"""
+import json
+import math
+import os
+
+try:
+    import openpyxl
+    from openpyxl.styles import Font, PatternFill, Alignment, Border, Side
+    from openpyxl.utils import get_column_letter
+except ImportError:
+    print("Installing openpyxl...")
+    import subprocess, sys
+    subprocess.check_call([sys.executable, "-m", "pip", "install", "openpyxl"])
+    import openpyxl
+    from openpyxl.styles import Font, PatternFill, Alignment, Border, Side
+    from openpyxl.utils import get_column_letter
+
+
+def shoelace_area(coords):
+    """Compute area in square degrees using the shoelace formula."""
+    n = len(coords)
+    area = 0.0
+    for i in range(n):
+        j = (i + 1) % n
+        area += coords[i][0] * coords[j][1]
+        area -= coords[j][0] * coords[i][1]
+    return abs(area) / 2.0
+
+
+def polygon_area_m2(ring):
+    """
+    Approximate polygon area in m² using spherical coordinates.
+    ring: list of [lon, lat] pairs
+    """
+    R = 6371000  # Earth radius in meters
+    area_deg2 = shoelace_area(ring)
+    # Convert to m² using mean latitude
+    mean_lat = sum(p[1] for p in ring) / len(ring)
+    lat_rad = math.radians(mean_lat)
+    # 1 degree lat ≈ R * pi/180 meters, 1 degree lon ≈ R * cos(lat) * pi/180 meters
+    m_per_deg_lat = R * math.pi / 180
+    m_per_deg_lon = R * math.cos(lat_rad) * math.pi / 180
+    return area_deg2 * m_per_deg_lat * m_per_deg_lon
+
+
+def feature_area(feature):
+    """Compute total area of a feature (MultiPolygon or Polygon) in m²."""
+    geom = feature["geometry"]
+    total = 0.0
+    if geom["type"] == "MultiPolygon":
+        for polygon in geom["coordinates"]:
+            # First ring is outer boundary
+            total += polygon_area_m2(polygon[0])
+    elif geom["type"] == "Polygon":
+        total += polygon_area_m2(geom["coordinates"][0])
+    return total
+
+
+# Load GeoJSON
+geojson_path = r"C:\Users\timon\Documents\SmartCane_code\laravel_app\storage\app\angata\pivot.geojson"
+with open(geojson_path, "r", encoding="utf-8") as f:
+    gj = json.load(f)
+
+features = gj["features"]
+print(f"Total features: {len(features)}")
+
+# Compute areas and sort
+fields = []
+for feat in features:
+    field_name = feat["properties"].get("field", "?")
+    area = feature_area(feat)
+    fields.append({"field": field_name, "area_m2": area, "area_ha": area / 10000})
+
+fields.sort(key=lambda x: x["area_m2"], reverse=True)
+
+# Print top 10 for verification
+print("\nTop 10 fields by area:")
+for i, f in enumerate(fields[:10]):
+    print(f"  {i+1:3d}. Field {f['field']:15s} {f['area_ha']:.2f} ha")
+
+# Split: index 0,3,6,... → Timon; 1,4,7,... → Joey; 2,5,8,... → Dimitra
+timon   = [f for i, f in enumerate(fields) if i % 3 == 0]
+joey    = [f for i, f in enumerate(fields) if i % 3 == 1]
+dimitra = [f for i, f in enumerate(fields) if i % 3 == 2]
+
+print(f"\nSplit: Timon={len(timon)}, Joey={len(joey)}, Dimitra={len(dimitra)}")
+
+# Create Excel
+wb = openpyxl.Workbook()
+ws = wb.active
+ws.title = "Field Checklist"
+
+# Color palette
+colors = {
+    "timon":   {"header_bg": "1F6AA5", "alt_bg": "D6E4F0"},
+    "joey":    {"header_bg": "2E7D32", "alt_bg": "D7F0D8"},
+    "dimitra": {"header_bg": "7B1FA2", "alt_bg": "EDD7F0"},
+}
+
+white_font = Font(name="Calibri", bold=True, color="FFFFFF", size=11)
+black_font = Font(name="Calibri", size=10)
+bold_black  = Font(name="Calibri", bold=True, size=10)
+center_align = Alignment(horizontal="center", vertical="center", wrap_text=True)
+
+thin = Side(style="thin", color="CCCCCC")
+border = Border(left=thin, right=thin, top=thin, bottom=thin)
+
+# Column layout: each person gets 3 columns (Rank, Field, Area ha, Checked)
+# Spacing: col 1-4 = Timon, col 5 = spacer, col 6-9 = Joey, col 10 = spacer, col 11-14 = Dimitra
+persons = [
+    ("Timon",   timon,   1),
+    ("Joey",    joey,    6),
+    ("Dimitra", dimitra, 11),
+]
+
+# Row 1: Person name header (merged across 4 cols)
+ws.row_dimensions[1].height = 25
+ws.row_dimensions[2].height = 20
+
+for name, data, start_col in persons:
+    c = colors[name.lower()]
+    # Merge cells for person name
+    end_col = start_col + 3
+    ws.merge_cells(
+        start_row=1, start_column=start_col,
+        end_row=1, end_column=end_col
+    )
+    cell = ws.cell(row=1, column=start_col, value=name)
+    cell.font = white_font
+    cell.fill = PatternFill("solid", fgColor=c["header_bg"])
+    cell.alignment = center_align
+    cell.border = border
+
+    # Sub-headers
+    sub_headers = ["#", "Field", "Area (ha)", "Checked ✓"]
+    for i, hdr in enumerate(sub_headers):
+        cell = ws.cell(row=2, column=start_col + i, value=hdr)
+        cell.font = white_font
+        cell.fill = PatternFill("solid", fgColor=c["header_bg"])
+        cell.alignment = center_align
+        cell.border = border
+
+    # Data rows
+    for row_i, field in enumerate(data):
+        row_num = row_i + 3
+        ws.row_dimensions[row_num].height = 16
+        alt = row_i % 2 == 1
+        bg = c["alt_bg"] if alt else "FFFFFF"
+        fill = PatternFill("solid", fgColor=bg)
+
+        values = [row_i + 1, field["field"], round(field["area_ha"], 2), ""]
+        for col_i, val in enumerate(values):
+            cell = ws.cell(row=row_num, column=start_col + col_i, value=val)
+            cell.font = black_font
+            cell.fill = fill
+            cell.border = border
+            cell.alignment = Alignment(horizontal="center" if col_i != 1 else "left",
+                                       vertical="center")
+
+# Column widths
+col_widths = {1: 5, 2: 14, 3: 10, 4: 12,   # Timon
+              5: 2,                            # spacer
+              6: 5, 7: 14, 8: 10, 9: 12,     # Joey
+              10: 2,                           # spacer
+              11: 5, 12: 14, 13: 10, 14: 12} # Dimitra
+
+for col, width in col_widths.items():
+    ws.column_dimensions[get_column_letter(col)].width = width
+
+# Freeze header rows
+ws.freeze_panes = "A3"
+
+# Title row above everything — insert a title row
+ws.insert_rows(1)
+ws.merge_cells("A1:N1")
+title_cell = ws.cell(row=1, column=1, value="Angata Pivot Field Checklist — sorted largest to smallest")
+title_cell.font = Font(name="Calibri", bold=True, size=13, color="1A1A1A")
+title_cell.alignment = Alignment(horizontal="center", vertical="center")
+title_cell.fill = PatternFill("solid", fgColor="F0F0F0")
+ws.row_dimensions[1].height = 28
+
+# Re-freeze after insert
+ws.freeze_panes = "A4"
+
+out_path = r"C:\Users\timon\Documents\SmartCane_code\angata_field_checklist.xlsx"
+wb.save(out_path)
+print(f"\nExcel saved to: {out_path}")
+print(f"Total fields: {len(fields)}")
+print(f"  Timon:   {len(timon)} fields")
+print(f"  Joey:    {len(joey)} fields")
+print(f"  Dimitra: {len(dimitra)} fields")

python_app/clean_empty_tiles.py

@@ -0,0 +1,154 @@
+"""
+Clean empty field-tile TIFFs and orphaned RDS files
+====================================================
+Scans field_tiles/ and/or field_tiles_CI/ directories and identifies TIF files
+where ALL pixels have RGBNIR == 0 (no satellite data collected).
+
+Partially-covered tiles (some valid pixels present) are kept.
+
+When deleting from field_tiles_CI/, also deletes the paired
+daily_ci_vals/{FIELD}/{DATE}.rds file if it exists.
+
+USAGE:
+  # Dry run — list empty files (default, scans both dirs):
+  & "C:\\Users\\timon\\anaconda3\\envs\\pytorch_gpu\\python.exe" python_app/clean_empty_tiles.py
+
+  # Actually delete:
+  & "...\\python.exe" python_app/clean_empty_tiles.py --delete
+
+  # Only one directory type:
+  & "...\\python.exe" python_app/clean_empty_tiles.py --dirs field_tiles
+  & "...\\python.exe" python_app/clean_empty_tiles.py --dirs field_tiles_CI
+
+  # Specific projects or fields:
+  & "...\\python.exe" python_app/clean_empty_tiles.py --projects angata aura --delete
+  & "...\\python.exe" python_app/clean_empty_tiles.py --fields 544 301 --delete
+"""
+
+import argparse
+from pathlib import Path
+
+import numpy as np
+import rasterio
+
+ROOT = Path(__file__).resolve().parent.parent
+DEFAULT_DIRS = ["field_tiles", "field_tiles_CI"]
+
+
+def is_empty_tif(path: Path) -> bool:
+    """Return True if ALL pixels in RGBNIR bands are 0 or NaN (no satellite data).
+
+    Cloud-masked pixels are stored as 0 in uint16 (NaN is not representable).
+    A tile is considered empty only when every pixel across bands 1-4 is 0 or NaN,
+    meaning no valid satellite data was captured for that field on that date.
+    Partially-covered tiles (some pixels valid) return False and are left alone.
+    """
+    try:
+        with rasterio.open(path) as src:
+            if src.count < 4:
+                return False  # unexpected band count — leave it alone
+            rgbnir = src.read([1, 2, 3, 4]).astype(np.float32)
+    except Exception as e:
+        print(f"  WARNING: could not open {path.name}: {e}")
+        return False
+
+    return bool(np.all((rgbnir == 0) | np.isnan(rgbnir)))
+
+
+def scan_directory(storage_root: Path, dir_name: str, delete: bool, fields: list = None) -> dict:
+    """Scan one tile directory within a project storage root.
+
+    When dir_name == 'field_tiles_CI' and delete=True, also removes the paired
+    daily_ci_vals/{FIELD}/{DATE}.rds file for each deleted TIF.
+
+    Returns:
+        dict mapping field_id -> list of empty Path objects
+    """
+    tiff_root = storage_root / dir_name
+    # Paired RDS files only exist for field_tiles_CI output
+    rds_root = storage_root / "daily_ci_vals" if dir_name == "field_tiles_CI" else None
+
+    if not tiff_root.exists():
+        print(f"  [{dir_name}] Directory not found: {tiff_root}")
+        return {}
+
+    field_dirs = sorted(d for d in tiff_root.iterdir() if d.is_dir())
+    if fields:
+        field_dirs = [d for d in field_dirs if d.name in fields]
+
+    print(f"\n  [{dir_name}] Scanning {len(field_dirs)} fields ...")
+
+    results = {}
+    for field_dir in field_dirs:
+        tif_files = sorted(field_dir.glob("*.tif"))
+        empty = [f for f in tif_files if is_empty_tif(f)]
+        if empty:
+            results[field_dir.name] = empty
+            print(f"    Field {field_dir.name:>6}: {len(empty)}/{len(tif_files)} empty"
+                  f"  ({', '.join(f.stem for f in empty)})")
+
+    total_empty = sum(len(v) for v in results.values())
+    total_tifs = sum(len(list(d.glob("*.tif"))) for d in field_dirs)
+    print(f"\n  [{dir_name}] Summary: {total_empty} empty / {total_tifs} total TIFs"
+          f" across {len(results)} fields")
+
+    if delete and total_empty > 0:
+        print(f"\n  [{dir_name}] Deleting {total_empty} empty TIFs ...")
+        rds_deleted = 0
+        for field_id, files in results.items():
+            for f in files:
+                f.unlink()
+                print(f"    Deleted TIF: {f.relative_to(ROOT)}")
+                # Also remove the paired RDS from daily_ci_vals/ (Script 20 output)
+                if rds_root is not None:
+                    paired_rds = rds_root / field_id / f"{f.stem}.rds"
+                    if paired_rds.exists():
+                        paired_rds.unlink()
+                        print(f"    Deleted RDS: {paired_rds.relative_to(ROOT)}")
+                        rds_deleted += 1
+        print(f"  [{dir_name}] Done. ({rds_deleted} paired RDS files also removed)")
+    elif not delete and total_empty > 0:
+        print(f"\n  [{dir_name}] Dry run — pass --delete to remove these files.")
+
+    return results
+
+
+def scan_project(project: str, delete: bool, fields: list = None, dirs: list = None) -> None:
+    storage_root = ROOT / "laravel_app" / "storage" / "app" / project
+    if not storage_root.exists():
+        print(f"[{project}] Project directory not found: {storage_root}")
+        return
+    print(f"\n[{project}] ========================================")
+    for dir_name in (dirs or DEFAULT_DIRS):
+        scan_directory(storage_root, dir_name, delete, fields)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Remove empty field-tile TIFFs and paired RDS files"
+    )
+    parser.add_argument(
+        "--delete", action="store_true",
+        help="Actually delete empty files (default: dry run)"
+    )
+    parser.add_argument(
+        "--projects", nargs="+", default=["angata"],
+        help="Project names to scan (default: angata)"
+    )
+    parser.add_argument(
+        "--fields", nargs="+", default=None,
+        help="Limit to specific field IDs, e.g. --fields 544 301"
+    )
+    parser.add_argument(
+        "--dirs", nargs="+", default=None, choices=DEFAULT_DIRS,
+        help=f"Which subdirs to scan (default: both {DEFAULT_DIRS})"
+    )
+    args = parser.parse_args()
+
+    print("=== Mode: DELETE ===" if args.delete else "=== Mode: DRY RUN (use --delete to remove) ===")
+    for project in args.projects:
+        scan_project(project, args.delete, args.fields, args.dirs)
+
+
+if __name__ == "__main__":
+    main()

python_app/weather_api_comparison.py

@@ -0,0 +1,973 @@
+"""
+weather_api_comparison.py
+=========================
+Compare daily precipitation from multiple free weather APIs across two locations:
+  - Arnhem, Netherlands  (51.985°N, 5.899°E)  — European climate
+  - Angata, Kenya        ( 1.330°S, 34.738°E) — tropical / sugarcane context
+
+ARCHIVE providers (no API key required):
+  1. Open-Meteo ERA5          — current SmartCane provider  (0.25°, global)
+  2. Open-Meteo ERA5-Land     — higher resolution variant   (0.10°, global)
+  3. Open-Meteo CERRA         — EU regional reanalysis      (0.05°, EU only)
+  4. NASA POWER               — completely independent source (0.50°, global)
+
+FORECAST providers (no API key required):
+  5. Open-Meteo Forecast      — deterministic NWP (global)
+  6. Open-Meteo Ensemble      — ECMWF IFS 51-member ensemble; gives probability bands
+  7. YR.no LocationForecast   — Norwegian Met Institute (~10 days, global)
+
+FORECAST providers (API key required — set in CONFIG below, leave "" to skip):
+  8. OpenWeatherMap            — free tier, 1000 calls/day
+  9. WeatherAPI.com            — free tier
+
+OUTPUT:
+  Plots saved to: weather_comparison_plots/
+    archive_<loc>.png              — daily lines + multi-source spread band + agreement signal
+    archive_rolling_<loc>.png      — 30-day rolling mean comparison (original style)
+    cumulative_<loc>.png           — cumulative annual precipitation
+    vs_era5_<loc>.png              — each provider vs ERA5 scatter (note: ERA5 is not ground truth)
+    pairwise_<loc>.png             — pairwise Pearson r and RMSE heatmaps (unbiased)
+    wetdry_agreement_<loc>.png     — % of days providers agree on wet vs dry
+    forecast_ensemble_<loc>.png    — ensemble uncertainty bands + exceedance probability
+    forecast_<loc>.png             — deterministic forecast bars (original style)
+
+Usage:
+  python weather_api_comparison.py
+"""
+
+import datetime
+import time
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.dates as mdates
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import requests
+from pathlib import Path
+
+
+# ============================================================
+# CONFIG
+# ============================================================
+
+LOCATIONS = {
+    "Arnhem_NL": {"lat": 51.985, "lon": 5.899},
+    "Angata_KE": {"lat": -1.330, "lon": 34.738},
+}
+
+# Archive: last 12 months
+ARCHIVE_END   = datetime.date.today() - datetime.timedelta(days=2)  # ERA5 lags ~2 days
+ARCHIVE_START = ARCHIVE_END - datetime.timedelta(days=365)
+
+# Forecast: today + 7 days
+FORECAST_START = datetime.date.today()
+FORECAST_END   = FORECAST_START + datetime.timedelta(days=7)
+
+# Optional API keys — leave "" to skip that provider
+OPENWEATHERMAP_KEY = ""   # https://openweathermap.org/api
+WEATHERAPI_KEY     = ""   # https://www.weatherapi.com/
+
+OUTPUT_DIR = Path("weather_comparison_plots")
+OUTPUT_DIR.mkdir(exist_ok=True)
+
+USER_AGENT = "SmartCane-WeatherComparison/1.0 (research; contact via github)"
+
+# ============================================================
+# ARCHIVE FETCHERS
+# ============================================================
+
+def fetch_openmeteo_archive(lat, lon, start, end, model="era5"):
+    """Open-Meteo ERA5 / ERA5-Land / CERRA archive.
+    ERA5 is the default (no models param needed). ERA5-Land and CERRA use lowercase names.
+    """
+    model_suffix = "" if model == "era5" else f"&models={model}"
+    url = (
+        f"https://archive-api.open-meteo.com/v1/archive"
+        f"?latitude={lat}&longitude={lon}"
+        f"&daily=precipitation_sum"
+        f"&start_date={start}&end_date={end}"
+        f"{model_suffix}"
+        f"&timezone=UTC"
+    )
+    r = requests.get(url, timeout=30)
+    r.raise_for_status()
+    body = r.json()
+    df = pd.DataFrame({
+        "date":    pd.to_datetime(body["daily"]["time"]),
+        "rain_mm": body["daily"]["precipitation_sum"],
+    })
+    df["rain_mm"] = pd.to_numeric(df["rain_mm"], errors="coerce").clip(lower=0).fillna(0)
+    # ERA5-Land sometimes returns values in meters (Open-Meteo API quirk).
+    # Auto-detect: if annual total < 50mm for any non-polar location, assume m → convert.
+    if df["rain_mm"].sum() < 50 and len(df) > 30:
+        df["rain_mm"] = df["rain_mm"] * 1000
+        print(f"    ⚠ Unit auto-converted m→mm (values were implausibly small)")
+    return df
+
+
+def fetch_nasa_power(lat, lon, start, end):
+    """NASA POWER — daily PRECTOTCORR (precipitation corrected), 0.5° grid."""
+    url = (
+        "https://power.larc.nasa.gov/api/temporal/daily/point"
+        f"?parameters=PRECTOTCORR"
+        f"&community=AG"
+        f"&longitude={lon}&latitude={lat}"
+        f"&start={start.strftime('%Y%m%d')}&end={end.strftime('%Y%m%d')}"
+        f"&format=JSON"
+    )
+    r = requests.get(url, timeout=60)
+    r.raise_for_status()
+    body = r.json()
+    raw = body["properties"]["parameter"]["PRECTOTCORR"]
+    df = pd.DataFrame([
+        {"date": pd.to_datetime(k, format="%Y%m%d"), "rain_mm": max(v, 0)}
+        for k, v in raw.items()
+        if v != -999  # NASA POWER fill value
+    ])
+    return df.sort_values("date").reset_index(drop=True)
+
+
+# ============================================================
+# FORECAST FETCHERS
+# ============================================================
+
+def fetch_openmeteo_forecast(lat, lon, days=8):
+    """Open-Meteo NWP forecast — default best_match model."""
+    url = (
+        f"https://api.open-meteo.com/v1/forecast"
+        f"?latitude={lat}&longitude={lon}"
+        f"&daily=precipitation_sum"
+        f"&forecast_days={days + 1}"
+        f"&timezone=UTC"
+    )
+    r = requests.get(url, timeout=30)
+    r.raise_for_status()
+    body = r.json()
+    df = pd.DataFrame({
+        "date":    pd.to_datetime(body["daily"]["time"]),
+        "rain_mm": body["daily"]["precipitation_sum"],
+    })
+    df["rain_mm"] = pd.to_numeric(df["rain_mm"], errors="coerce").fillna(0)
+    return df
+
+
+def fetch_openmeteo_ensemble_forecast(lat, lon, days=8):
+    """Open-Meteo ECMWF IFS ensemble forecast — 51 members at 0.4° resolution.
+
+    Returns a DataFrame with daily percentile bands and exceedance probabilities:
+      p10, p25, p50 (median), p75, p90  — mm/day
+      prob_gt_1mm, prob_gt_5mm          — % of members exceeding threshold
+      n_members                          — number of members in response
+    """
+    url = (
+        f"https://ensemble-api.open-meteo.com/v1/ensemble"
+        f"?latitude={lat}&longitude={lon}"
+        f"&daily=precipitation_sum"
+        f"&models=ecmwf_ifs04"
+        f"&forecast_days={days}"
+        f"&timezone=UTC"
+    )
+    r = requests.get(url, timeout=60)
+    r.raise_for_status()
+    body = r.json()
+    daily = body["daily"]
+    dates = pd.to_datetime(daily["time"])
+
+    # Member columns are named like "precipitation_sum_member01", etc.
+    member_keys = [k for k in daily.keys() if k.startswith("precipitation_sum")]
+    if not member_keys:
+        print("    ⚠ No member columns found in ensemble response")
+        return None
+
+    # Shape: (n_members, n_days)
+    members = np.array([daily[k] for k in member_keys], dtype=float)
+    members = np.where(members < 0, 0, members)  # clip negatives
+
+    df = pd.DataFrame({
+        "date":        dates,
+        "p10":         np.nanpercentile(members, 10, axis=0),
+        "p25":         np.nanpercentile(members, 25, axis=0),
+        "p50":         np.nanpercentile(members, 50, axis=0),
+        "p75":         np.nanpercentile(members, 75, axis=0),
+        "p90":         np.nanpercentile(members, 90, axis=0),
+        "mean":        np.nanmean(members, axis=0),
+        "prob_gt_1mm": np.mean(members > 1.0, axis=0) * 100,
+        "prob_gt_5mm": np.mean(members > 5.0, axis=0) * 100,
+        "n_members":   members.shape[0],
+    })
+    return df
+
+
+def fetch_yr_forecast(lat, lon):
+    """YR.no LocationForecast 2.0 — hourly precip aggregated to daily.
+    Note: forecast-only service; no historical archive is available.
+    """
+    url = f"https://api.met.no/weatherapi/locationforecast/2.0/compact?lat={lat}&lon={lon}"
+    headers = {"User-Agent": USER_AGENT}
+    r = requests.get(url, headers=headers, timeout=30)
+    r.raise_for_status()
+    body = r.json()
+
+    records = []
+    for entry in body["properties"]["timeseries"]:
+        ts   = pd.to_datetime(entry["time"])
+        data = entry["data"]
+        precip = 0.0
+        if "next_1_hours" in data:
+            precip = data["next_1_hours"]["details"].get("precipitation_amount", 0.0)
+        elif "next_6_hours" in data:
+            precip = data["next_6_hours"]["details"].get("precipitation_amount", 0.0) / 6
+        records.append({"datetime": ts, "precip_hour": precip})
+
+    hourly = pd.DataFrame(records)
+    hourly["date"] = hourly["datetime"].dt.date
+    daily = (
+        hourly.groupby("date")["precip_hour"]
+        .sum()
+        .reset_index()
+        .rename(columns={"precip_hour": "rain_mm"})
+    )
+    daily["date"] = pd.to_datetime(daily["date"])
+    return daily
+
+
+def fetch_openweathermap_forecast(lat, lon, api_key):
+    """OpenWeatherMap One Call 3.0 — daily forecast (needs paid/free key)."""
+    url = (
+        f"https://api.openweathermap.org/data/3.0/onecall"
+        f"?lat={lat}&lon={lon}"
+        f"&exclude=current,minutely,hourly,alerts"
+        f"&appid={api_key}&units=metric"
+    )
+    r = requests.get(url, timeout=30)
+    r.raise_for_status()
+    body = r.json()
+    records = []
+    for day in body.get("daily", []):
+        records.append({
+            "date":    pd.to_datetime(day["dt"], unit="s").normalize(),
+            "rain_mm": day.get("rain", 0.0),
+        })
+    return pd.DataFrame(records)
+
+
+def fetch_weatherapi_forecast(lat, lon, api_key, days=7):
+    """WeatherAPI.com free forecast (up to 3 days on free tier, 14 on paid)."""
+    url = (
+        f"https://api.weatherapi.com/v1/forecast.json"
+        f"?key={api_key}&q={lat},{lon}&days={days}&aqi=no&alerts=no"
+    )
+    r = requests.get(url, timeout=30)
+    r.raise_for_status()
+    body = r.json()
+    records = []
+    for day in body["forecast"]["forecastday"]:
+        records.append({
+            "date":    pd.to_datetime(day["date"]),
+            "rain_mm": day["day"].get("totalprecip_mm", 0.0),
+        })
+    return pd.DataFrame(records)
+
+
+# ============================================================
+# STATS
+# ============================================================
+
+def compare_stats(df, ref_col, other_col):
+    """Compute MAE, RMSE, bias, Pearson r between two columns."""
+    valid = df[[ref_col, other_col]].dropna()
+    if len(valid) < 5:
+        return {"n": len(valid), "MAE": None, "RMSE": None, "Bias": None, "r": None}
+    diff  = valid[other_col] - valid[ref_col]
+    mae   = diff.abs().mean()
+    rmse  = (diff**2).mean()**0.5
+    bias  = diff.mean()
+    r     = valid[ref_col].corr(valid[other_col])
+    return {"n": len(valid), "MAE": round(mae,2), "RMSE": round(rmse,2),
+            "Bias": round(bias,2), "r": round(r,3)}
+
+
+def compute_pairwise_stats(data_dict):
+    """Compute pairwise Pearson r and RMSE for all archive provider pairs.
+
+    Returns:
+        names       — list of provider names (same order as matrices)
+        r_matrix    — (n x n) Pearson r array, NaN where not computable
+        rmse_matrix — (n x n) RMSE array (mm/day), 0 on diagonal
+    """
+    providers = [(name, df) for name, df in data_dict.items()
+                 if df is not None and len(df) > 5]
+    names = [p[0] for p in providers]
+    n = len(names)
+
+    r_matrix    = np.full((n, n), np.nan)
+    rmse_matrix = np.full((n, n), np.nan)
+
+    for i in range(n):
+        r_matrix[i, i]    = 1.0
+        rmse_matrix[i, i] = 0.0
+        _, df_i = providers[i]
+        for j in range(i + 1, n):
+            _, df_j = providers[j]
+            merged = df_i.merge(df_j, on="date", suffixes=("_i", "_j"))
+            valid = merged[["rain_mm_i", "rain_mm_j"]].dropna()
+            if len(valid) < 5:
+                continue
+            r    = valid["rain_mm_i"].corr(valid["rain_mm_j"])
+            rmse = ((valid["rain_mm_i"] - valid["rain_mm_j"]) ** 2).mean() ** 0.5
+            r_matrix[i, j]    = r_matrix[j, i]    = round(r, 3)
+            rmse_matrix[i, j] = rmse_matrix[j, i] = round(rmse, 2)
+
+    return names, r_matrix, rmse_matrix
+
+
+def wetdry_agreement(data_dict, threshold=1.0):
+    """For each provider pair, compute % of days both-dry / both-wet / disagree.
+
+    A day is 'dry' if rain_mm < threshold (default 1 mm).
+    Returns a list of dicts: pair, both_dry, both_wet, disagree, n.
+    """
+    providers = [(name, df) for name, df in data_dict.items()
+                 if df is not None and len(df) > 5]
+    results = []
+    for i, (name_i, df_i) in enumerate(providers):
+        for j in range(i + 1, len(providers)):
+            name_j, df_j = providers[j]
+            merged = df_i.merge(df_j, on="date", suffixes=("_i", "_j"))
+            valid = merged[["rain_mm_i", "rain_mm_j"]].dropna()
+            if len(valid) < 5:
+                continue
+            dry_i = valid["rain_mm_i"] < threshold
+            dry_j = valid["rain_mm_j"] < threshold
+            both_dry = (dry_i  &  dry_j).mean() * 100
+            both_wet = (~dry_i & ~dry_j).mean() * 100
+            disagree = 100 - both_dry - both_wet
+            results.append({
+                "pair":      f"{name_i}\nvs\n{name_j}",
+                "both_dry":  round(both_dry, 1),
+                "both_wet":  round(both_wet, 1),
+                "disagree":  round(disagree, 1),
+                "n":         len(valid),
+            })
+    return results
+
+
+# ============================================================
+# PLOTTING — ARCHIVE
+# ============================================================
+
+ARCHIVE_COLORS = {
+    "ERA5 (Open-Meteo)":      "#1f77b4",
+    "ERA5-Land (Open-Meteo)": "#ff7f0e",
+    "CERRA (Open-Meteo)":     "#2ca02c",
+    "NASA POWER":             "#d62728",
+}
+
+
+def _build_spread_frame(data_dict):
+    """Merge all valid archive providers onto a common date axis.
+    Returns a DataFrame with one column per provider + _mean and _std columns.
+    """
+    valid = {name: df for name, df in data_dict.items()
+             if df is not None and len(df) > 0}
+    if not valid:
+        return None, []
+
+    merged = None
+    for name, df in valid.items():
+        tmp = df[["date", "rain_mm"]].rename(columns={"rain_mm": name})
+        merged = tmp if merged is None else merged.merge(tmp, on="date", how="outer")
+
+    cols = list(valid.keys())
+    merged[cols]   = merged[cols].clip(lower=0)
+    merged["_mean"] = merged[cols].mean(axis=1)
+    merged["_std"]  = merged[cols].std(axis=1)
+    return merged.sort_values("date").reset_index(drop=True), cols
+
+
+def plot_archive_with_spread(data_dict, location_name, start, end, output_dir):
+    """Three-panel archive plot:
+      Top:    Individual provider lines + multi-source mean±std shading
+      Middle: 30-day rolling mean
+      Bottom: Inter-source std (agreement signal — used to flag uncertain periods)
+    """
+    spread, _ = _build_spread_frame(data_dict)
+    if spread is None:
+        return
+
+    valid_providers = {name: df for name, df in data_dict.items()
+                       if df is not None and len(df) > 0}
+
+    fig, axes = plt.subplots(3, 1, figsize=(14, 12), sharex=True)
+
+    # --- Top: daily raw + spread band ---
+    ax1 = axes[0]
+    ax1.fill_between(
+        spread["date"],
+        (spread["_mean"] - spread["_std"]).clip(lower=0),
+        spread["_mean"] + spread["_std"],
+        alpha=0.15, color="gray", label="Multi-source ±1 std"
+    )
+    ax1.plot(spread["date"], spread["_mean"], color="black", linewidth=1.0,
+             linestyle="--", alpha=0.6, label="Multi-source mean", zorder=5)
+    for name, df in valid_providers.items():
+        ax1.plot(df["date"], df["rain_mm"],
+                 label=name, color=ARCHIVE_COLORS.get(name, "gray"),
+                 linewidth=0.7, alpha=0.85)
+    ax1.set_ylabel("Precipitation (mm/day)")
+    ax1.set_title(
+        f"{location_name} — Daily Precipitation with Multi-Source Spread\n"
+        f"{start} → {end}  |  Grey band = ±1 std across all providers"
+    )
+    ax1.legend(fontsize=8, ncol=2)
+    ax1.grid(True, alpha=0.3)
+
+    # --- Middle: 30-day rolling mean ---
+    ax2 = axes[1]
+    roll_mean = spread.set_index("date")["_mean"].rolling(30, min_periods=15).mean()
+    roll_std  = spread.set_index("date")["_std"].rolling(30, min_periods=15).mean()
+    ax2.fill_between(
+        roll_mean.index,
+        (roll_mean - roll_std).clip(lower=0),
+        roll_mean + roll_std,
+        alpha=0.15, color="gray"
+    )
+    for name, df in valid_providers.items():
+        rolled = df.set_index("date")["rain_mm"].rolling(30, min_periods=15).mean()
+        ax2.plot(rolled.index, rolled.values,
+                 label=name, color=ARCHIVE_COLORS.get(name, "gray"), linewidth=1.5)
+    ax2.set_ylabel("30-day rolling mean (mm/day)")
+    ax2.legend(fontsize=8)
+    ax2.grid(True, alpha=0.3)
+
+    # --- Bottom: inter-source std (agreement signal) ---
+    ax3 = axes[2]
+    ax3.fill_between(
+        spread["date"], 0, spread["_std"],
+        color="purple", alpha=0.35,
+        label="Std across providers (higher = less agreement)"
+    )
+    median_std = spread["_std"].median()
+    ax3.axhline(median_std, color="purple", linestyle=":", linewidth=1,
+                label=f"Median std = {median_std:.2f} mm/day")
+    ax3.set_ylabel("Std dev across\nproviders (mm)")
+    ax3.set_xlabel("Date")
+    ax3.legend(fontsize=8)
+    ax3.grid(True, alpha=0.3)
+    ax3.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
+
+    fig.autofmt_xdate()
+    plt.tight_layout()
+    path = output_dir / f"archive_{location_name}.png"
+    plt.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+def plot_cumulative(data_dict, location_name, output_dir):
+    """Cumulative annual precipitation — most relevant for crop/irrigation context."""
+    fig, ax = plt.subplots(figsize=(14, 5))
+
+    for name, df in data_dict.items():
+        if df is None or len(df) == 0:
+            continue
+        s = df.set_index("date")["rain_mm"].sort_index().cumsum()
+        total = s.iloc[-1]
+        ax.plot(s.index, s.values,
+                label=f"{name}  (total: {total:.0f} mm)",
+                color=ARCHIVE_COLORS.get(name, "gray"), linewidth=1.8)
+
+    ax.set_ylabel("Cumulative precipitation (mm)")
+    ax.set_xlabel("Date")
+    ax.set_title(
+        f"{location_name} — Cumulative Annual Precipitation by Provider\n"
+        "Divergence = sources disagree on total seasonal rainfall"
+    )
+    ax.legend(fontsize=9)
+    ax.grid(True, alpha=0.3)
+    ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
+    fig.autofmt_xdate()
+    plt.tight_layout()
+
+    path = output_dir / f"cumulative_{location_name}.png"
+    plt.savefig(path, dpi=150)
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+def plot_vs_era5(data_dict, location_name, output_dir):
+    """Each provider vs ERA5 reference: scatter + regression line.
+
+    NOTE: ERA5 is used as reference here for visual comparison only.
+    It is a reanalysis product, not a ground-truth measurement.
+    Use the pairwise heatmap (plot_pairwise_heatmap) for an unbiased comparison.
+
+    How to read:
+      - Each panel shows one provider (y-axis) vs ERA5 (x-axis) for daily precip.
+      - Points on the red diagonal = perfect agreement.
+      - Points above = provider wetter than ERA5 on that day.
+      - r = Pearson correlation (1 = perfect). MAE = mean absolute error in mm/day.
+      - Bias = provider minus ERA5 on average (positive = provider wetter).
+    """
+    ref_name = "ERA5 (Open-Meteo)"
+    ref_df = data_dict.get(ref_name)
+    if ref_df is None:
+        return
+
+    others = [(n, df) for n, df in data_dict.items()
+              if n != ref_name and df is not None and len(df) > 0]
+    if not others:
+        return
+
+    n = len(others)
+    fig, axes = plt.subplots(1, n, figsize=(5 * n, 5), squeeze=False)
+
+    for i, (name, df) in enumerate(others):
+        ax = axes[0][i]
+        merged = ref_df.merge(df, on="date", suffixes=("_ref", "_cmp"))
+        valid  = merged[["rain_mm_ref", "rain_mm_cmp"]].dropna()
+
+        color = ARCHIVE_COLORS.get(name, "steelblue")
+        ax.scatter(valid["rain_mm_ref"], valid["rain_mm_cmp"],
+                   s=4, alpha=0.35, color=color)
+
+        lim = max(valid.max().max(), 1) * 1.05
+        ax.plot([0, lim], [0, lim], "r--", linewidth=1, label="Perfect agreement")
+
+        if len(valid) > 5:
+            coeffs = np.polyfit(valid["rain_mm_ref"], valid["rain_mm_cmp"], 1)
+            x_fit  = np.linspace(0, lim, 100)
+            ax.plot(x_fit, np.polyval(coeffs, x_fit), "k-", linewidth=1,
+                    alpha=0.6, label=f"Regression (slope={coeffs[0]:.2f})")
+
+        stats = compare_stats(merged, "rain_mm_ref", "rain_mm_cmp")
+        ax.set_xlim(0, lim); ax.set_ylim(0, lim)
+        ax.set_xlabel("ERA5 (Open-Meteo) mm/day", fontsize=9)
+        ax.set_ylabel(f"{name} mm/day", fontsize=9)
+        ax.set_title(
+            f"{name}\nr={stats['r']}  MAE={stats['MAE']} mm  Bias={stats['Bias']:+.2f} mm",
+            fontsize=9
+        )
+        ax.legend(fontsize=7)
+        ax.grid(True, alpha=0.3)
+
+    fig.suptitle(
+        f"{location_name} — Daily Precip vs ERA5 (reference only — ERA5 is not ground truth)\n"
+        "Red dashed = perfect agreement. See pairwise heatmap for unbiased comparison.",
+        fontsize=10
+    )
+    plt.tight_layout()
+    path = output_dir / f"vs_era5_{location_name}.png"
+    plt.savefig(path, dpi=150)
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+def plot_pairwise_heatmap(data_dict, location_name, output_dir):
+    """Two heatmaps side by side: pairwise Pearson r and pairwise RMSE.
+
+    No provider is treated as reference — all pairs compared equally.
+    High r + low RMSE = sources agree. Where they diverge, neither is ground truth.
+    """
+    names, r_matrix, rmse_matrix = compute_pairwise_stats(data_dict)
+    if len(names) < 2:
+        return
+
+    n = len(names)
+    short = [nm.replace(" (Open-Meteo)", "\n(OM)") for nm in names]
+
+    fig, axes = plt.subplots(1, 2, figsize=(max(10, n * 2.5), max(5, n * 1.5 + 1)))
+
+    # Pearson r
+    ax1 = axes[0]
+    im1 = ax1.imshow(r_matrix, vmin=0, vmax=1, cmap="YlGn", aspect="auto")
+    plt.colorbar(im1, ax=ax1, label="Pearson r")
+    ax1.set_xticks(range(n)); ax1.set_yticks(range(n))
+    ax1.set_xticklabels(short, fontsize=8, rotation=45, ha="right")
+    ax1.set_yticklabels(short, fontsize=8)
+    for i in range(n):
+        for j in range(n):
+            val = r_matrix[i, j]
+            if not np.isnan(val):
+                text_color = "white" if val > 0.8 else "black"
+                ax1.text(j, i, f"{val:.2f}", ha="center", va="center",
+                         fontsize=9, color=text_color)
+    ax1.set_title("Pairwise Pearson r\n(1.0 = perfect agreement, no reference)")
+
+    # RMSE
+    ax2 = axes[1]
+    off_diag = rmse_matrix[rmse_matrix > 0]
+    vmax_rmse = float(np.nanmax(off_diag)) if len(off_diag) > 0 else 10
+    im2 = ax2.imshow(rmse_matrix, vmin=0, vmax=vmax_rmse, cmap="YlOrRd_r", aspect="auto")
+    plt.colorbar(im2, ax=ax2, label="RMSE (mm/day)")
+    ax2.set_xticks(range(n)); ax2.set_yticks(range(n))
+    ax2.set_xticklabels(short, fontsize=8, rotation=45, ha="right")
+    ax2.set_yticklabels(short, fontsize=8)
+    for i in range(n):
+        for j in range(n):
+            val = rmse_matrix[i, j]
+            if not np.isnan(val):
+                ax2.text(j, i, f"{val:.1f}", ha="center", va="center", fontsize=9)
+    ax2.set_title("Pairwise RMSE (mm/day)\n(0 = perfect agreement)")
+
+    fig.suptitle(
+        f"{location_name} — Archive Provider Agreement (unbiased — no single reference)\n"
+        "Pairs with high r + low RMSE are consistent. Divergent pairs reveal dataset uncertainty.",
+        fontsize=10
+    )
+    plt.tight_layout()
+    path = output_dir / f"pairwise_{location_name}.png"
+    plt.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+def plot_wetdry_agreement(data_dict, location_name, output_dir, threshold=1.0):
+    """Stacked bar chart: for each provider pair, % of days both-dry / both-wet / disagree.
+
+    High 'disagree' % means one source says rain, the other says dry on the same day.
+    This is the most practically relevant divergence for SmartCane crop monitoring.
+    """
+    results = wetdry_agreement(data_dict, threshold)
+    if not results:
+        return
+
+    pairs    = [r["pair"]     for r in results]
+    both_dry = [r["both_dry"] for r in results]
+    both_wet = [r["both_wet"] for r in results]
+    disagree = [r["disagree"] for r in results]
+
+    fig, ax = plt.subplots(figsize=(max(8, len(pairs) * 2.8), 6))
+    x = np.arange(len(pairs))
+
+    ax.bar(x, both_dry, label=f"Both dry (<{threshold} mm)", color="#a8d5e2", edgecolor="white")
+    ax.bar(x, both_wet, bottom=both_dry,
+           label=f"Both wet (≥{threshold} mm)", color="#3a7ebf", edgecolor="white")
+    bottom2 = [d + w for d, w in zip(both_dry, both_wet)]
+    ax.bar(x, disagree, bottom=bottom2,
+           label="Disagree (one wet, one dry)", color="#e07b54", edgecolor="white")
+
+    # Add % labels inside bars
+    for i, r in enumerate(results):
+        if r["both_dry"] > 4:
+            ax.text(i, r["both_dry"] / 2, f"{r['both_dry']:.0f}%",
+                    ha="center", va="center", fontsize=7, color="black")
+        if r["both_wet"] > 4:
+            ax.text(i, r["both_dry"] + r["both_wet"] / 2, f"{r['both_wet']:.0f}%",
+                    ha="center", va="center", fontsize=7, color="white")
+        if r["disagree"] > 4:
+            ax.text(i, bottom2[i] + r["disagree"] / 2, f"{r['disagree']:.0f}%",
+                    ha="center", va="center", fontsize=7, color="white")
+
+    ax.set_xticks(x)
+    ax.set_xticklabels(pairs, fontsize=8)
+    ax.set_ylim(0, 108)
+    ax.set_ylabel("% of days in archive period")
+    ax.set_title(
+        f"{location_name} — Provider Agreement on Wet vs Dry Days\n"
+        f"Threshold = {threshold} mm/day  |  Orange = providers disagree on whether it rained"
+    )
+    ax.legend(loc="upper right", fontsize=9)
+    ax.grid(True, axis="y", alpha=0.3)
+    plt.tight_layout()
+
+    path = output_dir / f"wetdry_agreement_{location_name}.png"
+    plt.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+# ============================================================
+# PLOTTING — FORECAST
+# ============================================================
+
+def plot_forecast_with_ensemble(forecast_data, ensemble_df, location_name, output_dir):
+    """Two-panel forecast comparison:
+      Top:    Deterministic bars (YR.no, Open-Meteo) + ECMWF ensemble percentile shading
+      Bottom: Exceedance probabilities P(rain > 1 mm) and P(rain > 5 mm)
+
+    Reading guide:
+      - Dark blue shading = where 50% of ensemble members agree (25th–75th %ile)
+      - Light blue shading = where 80% of ensemble members agree (10th–90th %ile)
+      - Bars = deterministic (single-value) forecast from each provider
+      - Bottom panel: at 50% on the dashed line = coin-flip uncertainty about rain event
+    """
+    has_ensemble = ensemble_df is not None and len(ensemble_df) > 0
+    det_providers = [(name, df) for name, df in forecast_data.items()
+                     if df is not None and len(df) > 0]
+
+    if not has_ensemble and not det_providers:
+        return
+
+    fig, axes = plt.subplots(
+        2, 1, figsize=(12, 8), sharex=True,
+        gridspec_kw={"height_ratios": [3, 1.5]}
+    )
+
+    # ---- TOP: deterministic bars + ensemble shading ----
+    ax1 = axes[0]
+
+    if has_ensemble:
+        ax1.fill_between(
+            ensemble_df["date"], ensemble_df["p10"], ensemble_df["p90"],
+            alpha=0.12, color="steelblue", label="Ensemble 10th–90th %ile"
+        )
+        ax1.fill_between(
+            ensemble_df["date"], ensemble_df["p25"], ensemble_df["p75"],
+            alpha=0.28, color="steelblue", label="Ensemble 25th–75th %ile"
+        )
+        ax1.plot(
+            ensemble_df["date"], ensemble_df["p50"],
+            color="steelblue", linewidth=1.8, linestyle="-",
+            label=f"Ensemble median ({int(ensemble_df['n_members'].iloc[0])} members)",
+            zorder=5
+        )
+
+    det_colors = {"Open-Meteo Forecast": "#1f77b4", "YR.no": "#e8882a"}
+    bar_half   = datetime.timedelta(hours=10)  # offset so bars don't overlap
+    n_det = len(det_providers)
+    offsets = np.linspace(-bar_half.total_seconds() / 3600 * (n_det - 1) / 2,
+                           bar_half.total_seconds() / 3600 * (n_det - 1) / 2,
+                           n_det)
+    bar_width_days = 0.3
+
+    for k, (name, df) in enumerate(det_providers):
+        offset_td = datetime.timedelta(hours=offsets[k])
+        shifted_dates = df["date"] + offset_td
+        ax1.bar(
+            shifted_dates, df["rain_mm"],
+            width=bar_width_days,
+            label=name,
+            color=det_colors.get(name, f"C{k}"),
+            alpha=0.75,
+            zorder=4,
+        )
+
+    ax1.set_ylabel("Precipitation (mm/day)")
+    ax1.set_title(
+        f"{location_name} — 7-Day Forecast\n"
+        "Shading = Open-Meteo ECMWF ensemble spread  |  Bars = deterministic forecasts"
+    )
+    ax1.legend(fontsize=8, loc="upper right")
+    ax1.grid(True, axis="y", alpha=0.3)
+    ax1.set_ylim(bottom=0)
+
+    # ---- BOTTOM: exceedance probabilities ----
+    ax2 = axes[1]
+    if has_ensemble:
+        ax2.step(
+            ensemble_df["date"], ensemble_df["prob_gt_1mm"],
+            where="mid", color="steelblue", linewidth=2.0,
+            label="P(rain > 1 mm)"
+        )
+        ax2.step(
+            ensemble_df["date"], ensemble_df["prob_gt_5mm"],
+            where="mid", color="darkblue", linewidth=2.0, linestyle="--",
+            label="P(rain > 5 mm)"
+        )
+        ax2.fill_between(
+            ensemble_df["date"], 0, ensemble_df["prob_gt_5mm"],
+            step="mid", alpha=0.18, color="darkblue"
+        )
+        ax2.axhline(50, color="gray", linestyle=":", linewidth=0.9, alpha=0.8,
+                    label="50% (coin-flip)")
+        ax2.set_ylim(0, 108)
+        ax2.set_ylabel("Exceedance\nprobability (%)")
+        ax2.legend(fontsize=8, loc="upper right")
+        ax2.grid(True, alpha=0.3)
+    ax2.set_xlabel("Date")
+    ax2.xaxis.set_major_formatter(mdates.DateFormatter("%d %b"))
+
+    fig.autofmt_xdate()
+    plt.tight_layout()
+
+    path = output_dir / f"forecast_ensemble_{location_name}.png"
+    plt.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+def plot_forecast(data_dict, location_name, output_dir):
+    """Bar chart comparing deterministic 7-day forecasts (original style, kept for reference)."""
+    _, ax = plt.subplots(figsize=(12, 5))
+
+    providers = [(name, df) for name, df in data_dict.items() if df is not None and len(df) > 0]
+    n = len(providers)
+    if n == 0:
+        plt.close()
+        return
+
+    all_dates = sorted(set(
+        d for _, df in providers
+        for d in df["date"].dt.date.tolist()
+    ))
+    x = np.arange(len(all_dates))
+    width = 0.8 / n
+
+    cmap = matplotlib.colormaps["tab10"].resampled(n)
+    for i, (name, df) in enumerate(providers):
+        date_map = dict(zip(df["date"].dt.date, df["rain_mm"]))
+        vals = [date_map.get(d, 0.0) for d in all_dates]
+        ax.bar(x + i * width, vals, width, label=name, color=cmap(i), alpha=0.85)
+
+    ax.set_xticks(x + width * (n - 1) / 2)
+    ax.set_xticklabels([d.strftime("%d %b") for d in all_dates], rotation=45, ha="right")
+    ax.set_ylabel("Precipitation (mm/day)")
+    ax.set_title(f"{location_name} — 7-Day Deterministic Forecast Comparison")
+    ax.legend(fontsize=9)
+    ax.grid(True, axis="y", alpha=0.3)
+    plt.tight_layout()
+
+    path = output_dir / f"forecast_{location_name}.png"
+    plt.savefig(path, dpi=150)
+    plt.close()
+    print(f"  Saved: {path}")
+
+
+# ============================================================
+# MAIN
+# ============================================================
+
+def run_location(loc_name, lat, lon):
+    print(f"\n{'='*60}")
+    print(f"  {loc_name}  ({lat}°, {lon}°)")
+    print(f"{'='*60}")
+
+    # ---- ARCHIVE ----
+    print("\n[Archive]")
+    archive_data = {}
+
+    print("  Fetching Open-Meteo ERA5...")
+    try:
+        archive_data["ERA5 (Open-Meteo)"] = fetch_openmeteo_archive(
+            lat, lon, ARCHIVE_START, ARCHIVE_END, model="era5"
+        )
+        print(f"    → {len(archive_data['ERA5 (Open-Meteo)'])} days")
+    except Exception as e:
+        print(f"    ✗ ERA5 failed: {e}")
+        archive_data["ERA5 (Open-Meteo)"] = None
+
+    time.sleep(0.5)
+    print("  Fetching Open-Meteo ERA5-Land...")
+    try:
+        archive_data["ERA5-Land (Open-Meteo)"] = fetch_openmeteo_archive(
+            lat, lon, ARCHIVE_START, ARCHIVE_END, model="era5_land"
+        )
+        print(f"    → {len(archive_data['ERA5-Land (Open-Meteo)'])} days")
+    except Exception as e:
+        print(f"    ✗ ERA5-Land failed: {e}")
+        archive_data["ERA5-Land (Open-Meteo)"] = None
+
+    time.sleep(0.5)
+    # CERRA only covers Europe (roughly 20°W–45°E, 30°N–80°N)
+    if -20 <= lon <= 45 and 30 <= lat <= 80:
+        print("  Fetching Open-Meteo CERRA (EU only)...")
+        try:
+            archive_data["CERRA (Open-Meteo)"] = fetch_openmeteo_archive(
+                lat, lon, ARCHIVE_START, ARCHIVE_END, model="cerra"
+            )
+            print(f"    → {len(archive_data['CERRA (Open-Meteo)'])} days")
+        except Exception as e:
+            print(f"    ✗ CERRA failed: {e}")
+            archive_data["CERRA (Open-Meteo)"] = None
+    else:
+        print("  Skipping CERRA (outside EU coverage)")
+        archive_data["CERRA (Open-Meteo)"] = None
+
+    time.sleep(0.5)
+    print("  Fetching NASA POWER...")
+    try:
+        archive_data["NASA POWER"] = fetch_nasa_power(lat, lon, ARCHIVE_START, ARCHIVE_END)
+        print(f"    → {len(archive_data['NASA POWER'])} days")
+    except Exception as e:
+        print(f"    ✗ NASA POWER failed: {e}")
+        archive_data["NASA POWER"] = None
+
+    # Pairwise stats (no ERA5 reference bias)
+    print("\n  Pairwise archive stats (Pearson r | RMSE mm/day | Bias mm/day):")
+    names, r_mat, rmse_mat = compute_pairwise_stats(archive_data)
+    for i in range(len(names)):
+        for j in range(i + 1, len(names)):
+            if not np.isnan(r_mat[i, j]):
+                print(f"    {names[i]:30s} vs {names[j]:30s}  "
+                      f"r={r_mat[i,j]:.3f}  RMSE={rmse_mat[i,j]:.2f} mm")
+
+    plot_archive_with_spread(archive_data, loc_name, ARCHIVE_START, ARCHIVE_END, OUTPUT_DIR)
+    plot_cumulative(archive_data, loc_name, OUTPUT_DIR)
+    plot_vs_era5(archive_data, loc_name, OUTPUT_DIR)
+    plot_pairwise_heatmap(archive_data, loc_name, OUTPUT_DIR)
+    plot_wetdry_agreement(archive_data, loc_name, OUTPUT_DIR)
+
+    # ---- FORECAST ----
+    print("\n[Forecast]")
+    forecast_data = {}
+
+    print("  Fetching Open-Meteo deterministic forecast...")
+    try:
+        forecast_data["Open-Meteo Forecast"] = fetch_openmeteo_forecast(lat, lon)
+        print(f"    → {len(forecast_data['Open-Meteo Forecast'])} days")
+    except Exception as e:
+        print(f"    ✗ Open-Meteo forecast failed: {e}")
+        forecast_data["Open-Meteo Forecast"] = None
+
+    time.sleep(0.5)
+    print("  Fetching Open-Meteo ECMWF ensemble forecast...")
+    ensemble_df = None
+    try:
+        ensemble_df = fetch_openmeteo_ensemble_forecast(lat, lon)
+        if ensemble_df is not None:
+            print(f"    → {len(ensemble_df)} days  "
+                  f"({int(ensemble_df['n_members'].iloc[0])} ensemble members)")
+    except Exception as e:
+        print(f"    ✗ Open-Meteo ensemble failed: {e}")
+
+    time.sleep(0.5)
+    print("  Fetching YR.no LocationForecast...")
+    try:
+        forecast_data["YR.no"] = fetch_yr_forecast(lat, lon)
+        print(f"    → {len(forecast_data['YR.no'])} days")
+    except Exception as e:
+        print(f"    ✗ YR.no failed: {e}")
+        forecast_data["YR.no"] = None
+
+    if OPENWEATHERMAP_KEY:
+        time.sleep(0.5)
+        print("  Fetching OpenWeatherMap forecast...")
+        try:
+            forecast_data["OpenWeatherMap"] = fetch_openweathermap_forecast(
+                lat, lon, OPENWEATHERMAP_KEY
+            )
+            print(f"    → {len(forecast_data['OpenWeatherMap'])} days")
+        except Exception as e:
+            print(f"    ✗ OpenWeatherMap failed: {e}")
+            forecast_data["OpenWeatherMap"] = None
+
+    if WEATHERAPI_KEY:
+        time.sleep(0.5)
+        print("  Fetching WeatherAPI.com forecast...")
+        try:
+            forecast_data["WeatherAPI.com"] = fetch_weatherapi_forecast(
+                lat, lon, WEATHERAPI_KEY
+            )
+            print(f"    → {len(forecast_data['WeatherAPI.com'])} days")
+        except Exception as e:
+            print(f"    ✗ WeatherAPI.com failed: {e}")
+            forecast_data["WeatherAPI.com"] = None
+
+    plot_forecast_with_ensemble(forecast_data, ensemble_df, loc_name, OUTPUT_DIR)
+    plot_forecast(forecast_data, loc_name, OUTPUT_DIR)
+
+
+if __name__ == "__main__":
+    print(f"Weather API Comparison — {datetime.date.today()}")
+    print(f"Archive: {ARCHIVE_START} → {ARCHIVE_END}")
+    print(f"Forecast: {FORECAST_START} → {FORECAST_END}")
+    print(f"Output: {OUTPUT_DIR.resolve()}")
+
+    for loc_name, coords in LOCATIONS.items():
+        run_location(loc_name, coords["lat"], coords["lon"])
+        time.sleep(1)
+
+    print(f"\nDone. Plots saved to: {OUTPUT_DIR.resolve()}")

r_app/10_create_per_field_tiffs_utils.R

@@ -129,8 +129,18 @@ crop_tiff_to_fields <- function(tif_path, tif_date, fields, output_base_dir) {
     # Crop raster to field boundary
     tryCatch({
       field_rast <- crop(rast, field_geom)
-      writeRaster(field_rast, output_path, overwrite = TRUE)
-      created <- created + 1
+
+      # Skip empty tiles: cloud-masked pixels are stored as 0 in uint16
+      # (NaN cannot be represented in that format). A band sum of 0 means
+      # no satellite data was captured for this field on this date.
+      band_sums <- terra::global(field_rast, fun = "sum", na.rm = TRUE)
+      if (sum(band_sums$sum, na.rm = TRUE) == 0) {
+        safe_log(paste("SKIP (no data):", field_name, tif_date), "WARNING")
+        skipped <- skipped + 1
+      } else {
+        writeRaster(field_rast, output_path, overwrite = TRUE)
+        created <- created + 1
+      }
     }, error = function(e) {
       safe_log(paste("ERROR cropping field", field_name, ":", e$message), "ERROR")
       errors <<- errors + 1

r_app/20_ci_extraction_per_field.R

@@ -128,23 +128,60 @@ main <- function() {
   }
   
   # Process each DATE (load merged TIFF once, extract all fields from it)
-  total_success <- 0
+  total_processed_dates <- 0
+  total_skipped_dates <- 0
+  total_already_complete_dates <- 0
   total_error <- 0
   
   for (date_str in dates_filter) {
-    # Load the MERGED TIFF (farm-wide) ONCE for this date
     input_tif_merged <- file.path(setup$merged_tif_folder, sprintf("%s.tif", date_str))
+    output_tifs <- file.path(setup$field_tiles_ci_dir, fields, sprintf("%s.tif", date_str))
+    output_rds <- file.path(setup$daily_ci_vals_dir, fields, sprintf("%s.rds", date_str))
+    names(output_tifs) <- fields
+    names(output_rds) <- fields
+    tif_exists <- file.exists(output_tifs)
+    rds_exists <- file.exists(output_rds)
+    fields_need_rds_only <- fields[tif_exists & !rds_exists]
+    fields_need_raster <- fields[!tif_exists]
+    
+    if (length(fields_need_rds_only) == 0 && length(fields_need_raster) == 0) {
+      total_already_complete_dates <- total_already_complete_dates + 1
+      safe_log(sprintf("  %s: All field outputs already exist (skipping)", date_str))
+      next
+    }
+    
+    fields_processed_this_date <- 0
+    rds_only_processed <- 0
+    raster_processed_this_date <- 0
+    
+    if (length(fields_need_rds_only) > 0) {
+      for (field in fields_need_rds_only) {
+        tryCatch({
+          extract_rds_from_ci_tiff(output_tifs[[field]], output_rds[[field]], field_boundaries_sf, field)
+          fields_processed_this_date <- fields_processed_this_date + 1
+          rds_only_processed <- rds_only_processed + 1
+        }, error = function(e) {
+          safe_log(sprintf("    Error regenerating RDS for field %s: %s", field, e$message), "WARNING")
+        })
+      }
+    }
+    
+    if (length(fields_need_raster) == 0) {
+      total_processed_dates <- total_processed_dates + 1
+      safe_log(sprintf("  %s: Regenerated %d RDS files from existing CI TIFFs", date_str, rds_only_processed))
+      next
+    }
     
     if (!file.exists(input_tif_merged)) {
       safe_log(sprintf("  %s: merged_tif not found (skipping)", date_str))
-      total_error <<- total_error + 1
+      total_skipped_dates <- total_skipped_dates + 1
       next
     }
     
     tryCatch({
-      # Load 4-band TIFF ONCE
+      # Load the merged TIFF only when at least one field still needs a CI TIFF.
       raster_4band <- terra::rast(input_tif_merged)
-      safe_log(sprintf("  %s: Loaded merged TIFF, processing %d fields...", date_str, length(fields)))
+      safe_log(sprintf("  %s: Loaded merged TIFF, processing %d fields...", date_str, length(fields_need_raster)))
       
       # Calculate CI from 4-band
       ci_raster <- calc_ci_from_raster(raster_4band)
@@ -154,80 +191,75 @@ main <- function() {
       five_band <- c(raster_4band, ci_raster)
       names(five_band) <- c("Red", "Green", "Blue", "NIR", "CI")
       
-      # Now process all fields from this single merged TIFF
-      fields_processed_this_date <- 0
-      
-      for (field in fields) {
-        field_ci_path <- file.path(setup$field_tiles_ci_dir, field)
-        field_daily_vals_path <- file.path(setup$daily_ci_vals_dir, field)
+      # Now process only the fields that still need CI TIFF output for this date.
+      for (field in fields_need_raster) {
+        output_tif_path <- output_tifs[[field]]
+        output_rds_path <- output_rds[[field]]
         
-        # Pre-create output directories
-        dir.create(field_ci_path, showWarnings = FALSE, recursive = TRUE)
-        dir.create(field_daily_vals_path, showWarnings = FALSE, recursive = TRUE)
-        
-        output_tif <- file.path(field_ci_path, sprintf("%s.tif", date_str))
-        output_rds <- file.path(field_daily_vals_path, sprintf("%s.rds", date_str))
-        
-        # MODE 3: Skip if both outputs already exist
-        if (file.exists(output_tif) && file.exists(output_rds)) {
-          next
-        }
-        
-        # MODE 2: Regeneration mode - RDS missing but CI TIFF exists
-        if (file.exists(output_tif) && !file.exists(output_rds)) {
-          tryCatch({
-            extract_rds_from_ci_tiff(output_tif, output_rds, field_boundaries_sf, field)
-            fields_processed_this_date <- fields_processed_this_date + 1
-          }, error = function(e) {
-            # Continue to next field
-          })
-          next
-        }
-        
-        # MODE 1: Normal mode - crop 5-band TIFF to field boundary and save
         tryCatch({
           # Crop 5-band TIFF to field boundary
           field_geom <- field_boundaries_sf %>% filter(field == !!field)
           five_band_cropped <- terra::crop(five_band, field_geom, mask = TRUE)
-          
-          # Save 5-band field TIFF
-          terra::writeRaster(five_band_cropped, output_tif, overwrite = TRUE)
-          
-          # Extract CI statistics by sub_field (from cropped CI raster)
-          ci_cropped <- five_band_cropped[[5]]  # 5th band is CI
-          ci_stats <- extract_ci_by_subfield(ci_cropped, field_boundaries_sf, field)
-          
-          # Save RDS
-          if (!is.null(ci_stats) && nrow(ci_stats) > 0) {
-            saveRDS(ci_stats, output_rds)
+
+          # Skip empty tiles: cloud-masked pixels are stored as 0 in uint16
+          # (NaN cannot be represented in that format). Sum of RGBNIR bands == 0
+          # means no valid satellite data was captured for this field on this date.
+          rgbnir_sum <- sum(
+            terra::global(five_band_cropped[[1:4]], fun = "sum", na.rm = TRUE)$sum,
+            na.rm = TRUE
+          )
+          if (rgbnir_sum == 0) {
+            safe_log(sprintf("    SKIP (no data): field %s on %s", field, date_str), "WARNING")
+          } else {
+            # Save 5-band field TIFF
+            terra::writeRaster(five_band_cropped, output_tif_path, overwrite = TRUE)
+
+            # Extract CI statistics by sub_field (from cropped CI raster)
+            ci_cropped <- five_band_cropped[[5]]  # 5th band is CI
+            ci_stats <- extract_ci_by_subfield(ci_cropped, field_boundaries_sf, field)
+
+            # Save RDS
+            if (!is.null(ci_stats) && nrow(ci_stats) > 0) {
+              saveRDS(ci_stats, output_rds_path)
+            }
+
+            fields_processed_this_date <- fields_processed_this_date + 1
+            raster_processed_this_date <- raster_processed_this_date + 1
           }
-          
-          fields_processed_this_date <- fields_processed_this_date + 1
-          
         }, error = function(e) {
           # Error in individual field, continue to next
           safe_log(sprintf("    Error processing field %s: %s", field, e$message), "WARNING")
         })
       }
       
-      # Increment success counter if at least one field succeeded
       if (fields_processed_this_date > 0) {
-        total_success <<- total_success + 1
-        safe_log(sprintf("  %s: Processed %d fields", date_str, fields_processed_this_date))
+        total_processed_dates <- total_processed_dates + 1
+        safe_log(sprintf(
+          "  %s: Processed %d fields (%d CI TIFFs created, %d RDS-only regenerated)",
+          date_str,
+          fields_processed_this_date,
+          raster_processed_this_date,
+          rds_only_processed
+        ))
+      } else {
+        total_error <- total_error + 1
+        safe_log(sprintf("  %s: No field outputs were created; see warnings above", date_str), "ERROR")
       }
       
     }, error = function(e) {
-      total_error <<- total_error + 1
+      total_error <- total_error + 1
       safe_log(sprintf("  %s: Error loading or processing merged TIFF - %s", date_str, e$message), "ERROR")
     })
   }
   
   # Summary
   safe_log(sprintf("\n=== Processing Complete ==="))
-  safe_log(sprintf("Successfully processed: %d", total_success))
+  safe_log(sprintf("Dates processed: %d", total_processed_dates))
+  safe_log(sprintf("Dates skipped (missing merged_tif): %d", total_skipped_dates))
+  safe_log(sprintf("Dates already complete: %d", total_already_complete_dates))
   safe_log(sprintf("Errors encountered: %d", total_error))
   
-  if (total_success > 0) {
+  if (total_processed_dates > 0) {
     safe_log("Output files created in:")
     safe_log(sprintf("  TIFFs: %s", setup$field_tiles_ci_dir))
     safe_log(sprintf("  RDS:   %s", setup$daily_ci_vals_dir))

r_app/80_utils_cane_supply.R

@@ -170,7 +170,7 @@ calculate_status_alert <- function(imminent_prob, age_week, mean_ci,
   
   # Priority 1: HARVEST READY - highest business priority
   # Field is mature (≥12 months) AND harvest model predicts imminent harvest
-  if (!is.na(imminent_prob) && imminent_prob > 0.5 && !is.na(age_week) && age_week >= 52) {
+  if (!is.na(imminent_prob) && imminent_prob > 0.3 && !is.na(age_week) && age_week >= 52) {
     return("harvest_ready")
   }
   

r_app/80_utils_common.R

@@ -650,48 +650,6 @@ get_phase_by_age <- function(age_weeks) {
   return("Unknown")
 }
 
-#' Get status trigger based on CI values and field age
-get_status_trigger <- function(ci_values, ci_change, age_weeks) {
-  if (is.na(age_weeks) || length(ci_values) == 0) return(NA_character_)
-  
-  ci_values <- ci_values[!is.na(ci_values)]
-  if (length(ci_values) == 0) return(NA_character_)
-  
-  pct_above_2 <- sum(ci_values > 2) / length(ci_values) * 100
-  pct_at_or_above_2 <- sum(ci_values >= 2) / length(ci_values) * 100
-  ci_cv <- if (mean(ci_values, na.rm = TRUE) > 0) sd(ci_values) / mean(ci_values, na.rm = TRUE) else 0
-  mean_ci <- mean(ci_values, na.rm = TRUE)
-  
-  if (age_weeks >= 0 && age_weeks <= 6) {
-    if (pct_at_or_above_2 >= 70) {
-      return("germination_complete")
-    } else if (pct_above_2 > 10) {
-      return("germination_started")
-    }
-  }
-  
-  if (age_weeks >= 45) {
-    return("harvest_ready")
-  }
-  
-  if (age_weeks > 6 && !is.na(ci_change) && ci_change < -1.5 && ci_cv < 0.25) {
-    return("stress_detected_whole_field")
-  }
-  
-  if (age_weeks > 6 && !is.na(ci_change) && ci_change > 1.5) {
-    return("strong_recovery")
-  }
-  
-  if (age_weeks >= 4 && age_weeks < 39 && !is.na(ci_change) && ci_change > 0.2) {
-    return("growth_on_track")
-  }
-  
-  if (age_weeks >= 39 && age_weeks < 45 && mean_ci > 3.5) {
-    return("maturation_progressing")
-  }
-  
-  return(NA_character_)
-}
 
 #' Extract planting dates from harvesting data
 extract_planting_dates <- function(harvesting_data, field_boundaries_sf = NULL) {

r_app/90_CI_report_with_kpis_agronomic_support.Rmd

@@ -56,7 +56,6 @@ suppressPackageStartupMessages({
   # Visualization
   library(tmap)       # For interactive maps (field boundary visualization)
   library(ggspatial)  # For basemap tiles and spatial annotations (OSM basemap with ggplot2)
-  
   # Reporting
   library(knitr)      # For R Markdown document generation (code execution and output)
   library(flextable)  # For formatted tables in Word output (professional table styling)
@@ -84,6 +83,16 @@ tryCatch({
   })
 })
 
+tryCatch({
+  source("90_rainfall_utils.R")
+}, error = function(e) {
+  tryCatch({
+    source(here::here("r_app", "90_rainfall_utils.R"))
+  }, error = function(e) {
+    message("Could not load 90_rainfall_utils.R - rain overlay disabled: ", e$message)
+  })
+})
+
 ```
 
 ```{r initialize_project_config, message=FALSE, warning=FALSE, include=FALSE}
@@ -590,7 +599,17 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
   if (!is.na(total_fields)) {
     cat("\n\n", tr_key("tot_fields_analyzed"))
   }
- 
+
+  # Total area analyzed
+  if (exists("field_details_table") && !is.null(field_details_table)) {
+    area_col_name <- paste0("Area_", get_area_unit_label(AREA_UNIT_PREFERENCE))
+    unit_label    <- get_area_unit_label(AREA_UNIT_PREFERENCE)
+    if (area_col_name %in% names(field_details_table)) {
+      total_area <- sum(field_details_table[[area_col_name]], na.rm = TRUE)
+      cat("\n\n", tr_key("tot_area_analyzed"))
+    }
+  }
+
 } else {
   cat(tr_key("kpi_na"))
 }
@@ -609,15 +628,89 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
   tryCatch({
     # KPI metadata for display
     kpi_display_order <- list(
-      uniformity = list(display = "Field Uniformity", level_col = "interpretation", count_col = "field_count"),
-      area_change = list(display = "Area Change", level_col = "interpretation", count_col = "field_count"),
-      growth_decline = list(display = "4-Week Trend", level_col = "trend_interpretation", count_col = "field_count"),
-      patchiness = list(display = "Field Patchiness", level_col = "gini_category", count_col = "field_count", detail_col = "patchiness_risk"),
-      tch_forecast = list(display = "TCH Forecasted", level_col = "tch_category", detail_col = "range", count_col = "field_count"),
-      gap_filling = list(display = "Gaps", level_col = "gap_level", count_col = "field_count")
+      uniformity     = list(display = "Field Uniformity", level_col = "interpretation",     count_col = "field_count", area_col = "area_sum"),
+      area_change    = list(display = "Area Change",      level_col = "interpretation",     count_col = "field_count", area_col = "area_sum"),
+      growth_decline = list(display = "4-Week Trend",     level_col = "trend_interpretation", count_col = "field_count", area_col = "area_sum"),
+      patchiness     = list(display = "Field Patchiness", level_col = "gini_category",     count_col = "field_count", detail_col = "patchiness_risk", area_col = "area_sum"),
+      tch_forecast   = list(display = "TCH Forecasted",   level_col = "tch_category",      count_col = "field_count", detail_col = "range", area_col = "area_sum"),
+      gap_filling    = list(display = "Gaps",              level_col = "gap_level",         count_col = "field_count", area_col = "area_sum")
     )
 
-    standardize_kpi <- function(df, level_col, count_col, detail_col = NULL) {
+    # Enrich summary_tables with area_sum from field_details_table (mirrors script 91 pattern)
+    area_col_name <- paste0("Area_", get_area_unit_label(AREA_UNIT_PREFERENCE))
+    unit_label    <- get_area_unit_label(AREA_UNIT_PREFERENCE)
+    has_area <- exists("field_details_table") && !is.null(field_details_table) &&
+                area_col_name %in% names(field_details_table)
+    if (has_area) {
+      fdt <- field_details_table
+      if (!is.null(summary_tables$uniformity) && "Uniformity_Category" %in% names(fdt)) {
+        summary_tables$uniformity <- summary_tables$uniformity %>%
+          left_join(fdt %>% group_by(interpretation = Uniformity_Category) %>%
+                      summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop"),
+                    by = "interpretation")
+      }
+      if (!is.null(summary_tables$area_change) && "Area_Change_Interpretation" %in% names(fdt)) {
+        summary_tables$area_change <- summary_tables$area_change %>%
+          left_join(fdt %>% group_by(interpretation = Area_Change_Interpretation) %>%
+                      summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop"),
+                    by = "interpretation")
+      }
+      if (!is.null(summary_tables$growth_decline) && "Trend_Interpretation" %in% names(fdt)) {
+        summary_tables$growth_decline <- summary_tables$growth_decline %>%
+          left_join(fdt %>% group_by(trend_interpretation = Trend_Interpretation) %>%
+                      summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop"),
+                    by = "trend_interpretation")
+      }
+      if (!is.null(summary_tables$patchiness) && all(c("Patchiness_Risk", "Gini_Coefficient") %in% names(fdt))) {
+        summary_tables$patchiness <- summary_tables$patchiness %>%
+          left_join(
+            fdt %>%
+              mutate(gini_category = case_when(
+                Gini_Coefficient < 0.2 ~ "Uniform (Gini<0.2)",
+                Gini_Coefficient < 0.4 ~ "Moderate (Gini 0.2-0.4)",
+                TRUE                   ~ "High (Gini≥0.4)"
+              )) %>%
+              group_by(gini_category, patchiness_risk = Patchiness_Risk) %>%
+              summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop"),
+            by = c("gini_category", "patchiness_risk")
+          )
+      }
+      if (!is.null(summary_tables$gap_filling) && "Gap_Level" %in% names(fdt)) {
+        summary_tables$gap_filling <- summary_tables$gap_filling %>%
+          left_join(fdt %>% group_by(gap_level = Gap_Level) %>%
+                      summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop"),
+                    by = "gap_level")
+      }
+      # TCH forecast: reproduce the same quartile logic used when building summary_tables so we
+      # can assign each field to a tch_category and sum its area
+      if (!is.null(summary_tables$tch_forecast) && "TCH_Forecasted" %in% names(fdt)) {
+        tch_vals <- fdt %>% dplyr::filter(!is.na(TCH_Forecasted)) %>% dplyr::pull(TCH_Forecasted)
+        if (length(tch_vals) > 0) {
+          if (length(unique(tch_vals)) == 1) {
+            area_tch <- fdt %>%
+              dplyr::filter(!is.na(TCH_Forecasted)) %>%
+              dplyr::summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE)) %>%
+              dplyr::mutate(tch_category = "All equal")
+          } else {
+            q25 <- quantile(tch_vals, 0.25, na.rm = TRUE)
+            q75 <- quantile(tch_vals, 0.75, na.rm = TRUE)
+            area_tch <- fdt %>%
+              dplyr::filter(!is.na(TCH_Forecasted)) %>%
+              dplyr::mutate(tch_category = dplyr::case_when(
+                TCH_Forecasted >= q75 ~ "Top 25%",
+                TCH_Forecasted >= q25 ~ "Middle 50%",
+                TRUE                  ~ "Bottom 25%"
+              )) %>%
+              dplyr::group_by(tch_category) %>%
+              dplyr::summarise(area_sum = sum(.data[[area_col_name]], na.rm = TRUE), .groups = "drop")
+          }
+          summary_tables$tch_forecast <- summary_tables$tch_forecast %>%
+            left_join(area_tch, by = "tch_category")
+        }
+      }
+    }
+
+    standardize_kpi <- function(df, level_col, count_col, detail_col = NULL, area_col = NULL) {
       if (is.null(level_col) || !(level_col %in% names(df)) || is.null(count_col) || !(count_col %in% names(df))) {
         return(NULL)
       }
@@ -631,10 +724,16 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
         display_level <- df[[level_col]]
       }
 
+      area_vals <- if (!is.null(area_col) && area_col %in% names(df))
+                     round(df[[area_col]], 1)
+                   else
+                     rep(NA_real_, nrow(df))
+
       df %>%
         dplyr::transmute(
           Level = if (level_col == "trend_interpretation") map_trend_to_arrow(display_level, include_text = TRUE) else as.character(display_level),
           Count = as.integer(round(as.numeric(.data[[count_col]]))),
+          Area  = area_vals,
           Percent = if (is.na(total)) {
             NA_real_
           } else {
@@ -652,9 +751,10 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
       kpi_df <- summary_tables[[kpi_key]]
       if (is.null(kpi_df) || !is.data.frame(kpi_df) || nrow(kpi_df) == 0) next
 
-      # Pass detail_col if it exists in config
+      # Pass detail_col and area_col if present in config
       detail_col <- if (!is.null(config$detail_col)) config$detail_col else NULL
-      kpi_rows <- standardize_kpi(kpi_df, config$level_col, config$count_col, detail_col)
+      kpi_rows <- standardize_kpi(kpi_df, config$level_col, config$count_col, detail_col,
+                                  config[["area_col"]])
 
       if (!is.null(kpi_rows)) {
         kpi_rows$KPI <- config$display
@@ -676,10 +776,14 @@ if (exists("summary_tables") && !is.null(summary_tables) && length(summary_table
       kpi_group_sizes <- rle(combined_df$KPI_group)$lengths
 
       display_df <- combined_df %>%
-        dplyr::select(KPI = KPI_display, Level, Count, Percent)
+        dplyr::select(KPI = KPI_display, Level, Count, Area, Percent)
 
       # Translate the table for visualization
-      names(display_df) <- c(tr_key("KPI"), tr_key("Level"), tr_key("Count"), tr_key("Percent"))
+      names(display_df) <- c(
+        tr_key("KPI"), tr_key("Level"), tr_key("Count"),
+        paste0(tr_key("Area"), " (", unit_label, ")"),
+        tr_key("Percent")
+      )
       display_df[, 1:2] <- lapply(display_df[, 1:2], function(col) sapply(col, tr_key))      
       ft <- flextable(display_df) %>%
         merge_v(j = tr_key("KPI")) %>%
@@ -1179,11 +1283,52 @@ tryCatch({
   minus_2_ww <- get_week_year(as.Date(today) - lubridate::weeks(2))
   minus_3_ww <- get_week_year(as.Date(today) - lubridate::weeks(3))
   
-  message(paste("Processing", nrow(AllPivots_merged), "fields for weeks:", 
+  message(paste("Processing", nrow(AllPivots_merged), "fields for weeks:",
                 current_ww$week, minus_1_ww$week, minus_2_ww$week, minus_3_ww$week))
-  
+
+
   # load_per_field_mosaic() is defined in 90_report_utils.R (sourced above)
 
+  # --- Fetch rain for all fields (one API call per 0.25° tile) ---
+  all_rain_data <- NULL
+  if (exists("rain_fetch_for_fields") && !is.null(CI_quadrant) && nrow(CI_quadrant) > 0) {
+    tryCatch({
+      # Compute field centroids from AllPivots0 (already-loaded sf object)
+      field_centroids <- AllPivots0 %>%
+        dplyr::group_by(field) %>%
+        dplyr::summarise(geometry = sf::st_union(geometry), .groups = "drop") %>%
+        sf::st_centroid() %>%
+        dplyr::mutate(
+          longitude = sf::st_coordinates(geometry)[, 1],
+          latitude  = sf::st_coordinates(geometry)[, 2]
+        ) %>%
+        sf::st_drop_geometry() %>%
+        dplyr::select(field_id = field, latitude, longitude)
+
+      # Determine current-season date range across all fields
+      rain_start <- CI_quadrant %>%
+        dplyr::group_by(field) %>%
+        dplyr::filter(season == max(season)) %>%
+        dplyr::summarise(s = min(Date, na.rm = TRUE), .groups = "drop") %>%
+        dplyr::pull(s) %>%
+        min(na.rm = TRUE)
+      rain_end <- as.Date(today)
+
+      message(paste("Fetching rain data from", rain_start, "to", rain_end,
+                    "for", nrow(field_centroids), "fields"))
+
+      all_rain_data <- rain_fetch_for_fields(
+        centroids_df = field_centroids,
+        start_date   = rain_start,
+        end_date     = rain_end
+      )
+      message(paste("Rain data fetched:", nrow(all_rain_data), "rows"))
+    }, error = function(e) {
+      message("Rain fetch failed - overlay disabled: ", e$message)
+      all_rain_data <<- NULL
+    })
+  }
+
   # Iterate through fields using purrr::walk
   purrr::walk(AllPivots_merged$field, function(field_name) {
     tryCatch({
@@ -1249,6 +1394,14 @@ tryCatch({
       
       # Call cum_ci_plot for trend analysis
       if (!is.null(CI_quadrant)) {
+        field_rain <- if (!is.null(all_rain_data)) {
+          all_rain_data %>%
+            dplyr::filter(field_id == field_name) %>%
+            dplyr::select(date, rain_mm)
+        } else {
+          NULL
+        }
+
         cum_ci_plot(
           pivotName = field_name,
           ci_quadrant_data = ci_quadrant_data,
@@ -1259,7 +1412,8 @@ tryCatch({
           show_benchmarks = TRUE,
           estate_name = project_dir,
           benchmark_percentiles = c(10, 50, 90),
-          benchmark_data = benchmarks
+          benchmark_data = benchmarks,
+          rain_data = field_rain
         )
         #cat("\n")
       }
@@ -1275,6 +1429,14 @@ tryCatch({
             sprintf("**%s:** %.2f", tr_key("cv_value"), field_kpi$CV),
             sprintf("**%s:** %.2f", tr_key("mean_ci"), field_kpi$Mean_CI)
           )
+
+          # Prepend area as first item (static field attribute)
+          if (area_col_name %in% names(field_kpi) && !is.na(field_kpi[[area_col_name]])) {
+            kpi_parts <- c(
+              sprintf("**%s:** %.1f %s", tr_key("Area"), field_kpi[[area_col_name]], unit_label),
+              kpi_parts
+            )
+          }
           
           # Add Weekly_CI_Change if available (note: capital C and I)
           if (!is.null(field_kpi$Weekly_CI_Change) && !is.na(field_kpi$Weekly_CI_Change)) {

r_app/90_rainfall_utils.R

@@ -0,0 +1,376 @@
+# 90_RAINFALL_UTILS.R
+# ============================================================================
+# Rainfall data fetching and ggplot overlay utilities for SmartCane CI reports.
+#
+# Provider-agnostic design: to swap weather provider, implement a new
+# rain_fetch_<provider>() function and update the dispatch in rain_fetch().
+#
+# EXPORTS:
+#   - rain_snap_to_tile()          Snap lat/lon to Open-Meteo native 0.25° grid
+#   - rain_fetch()                 Generic fetch wrapper (dispatches to provider)
+#   - rain_fetch_for_fields()      Batch fetch with spatial tile deduplication
+#   - rain_join_to_ci()            Join rain to latest-season CI data by Date
+#   - rain_add_to_plot()           Overlay rain on single CI plot (abs or cum)
+#   - rain_add_to_faceted_plot()   Overlay rain on faceted "both" CI plot
+# ============================================================================
+
+suppressPackageStartupMessages({
+  library(dplyr)
+  library(tidyr)
+  library(httr2)
+  library(purrr)
+})
+
+# ============================================================================
+# PROVIDER LAYER — swap here to change data source
+# ============================================================================
+
+#' Fetch daily precipitation from Open-Meteo ERA5 archive
+#'
+#' ERA5 covers 1940-present at 0.25° (~28 km) resolution.
+#' No API key required.
+#'
+#' @param latitude  Numeric. WGS84 latitude.
+#' @param longitude Numeric. WGS84 longitude.
+#' @param start_date Date or "YYYY-MM-DD" string.
+#' @param end_date   Date or "YYYY-MM-DD" string.
+#' @return data.frame with columns: date (Date), rain_mm (numeric)
+rain_fetch_open_meteo_archive <- function(latitude, longitude, start_date, end_date) {
+  url <- paste0(
+    "https://archive-api.open-meteo.com/v1/archive",
+    "?latitude=",    latitude,
+    "&longitude=",   longitude,
+    "&daily=precipitation_sum",
+    "&start_date=",  format(as.Date(start_date), "%Y-%m-%d"),
+    "&end_date=",    format(as.Date(end_date),   "%Y-%m-%d"),
+    "&timezone=auto"
+  )
+
+  resp <- tryCatch(
+    httr2::request(url) %>% httr2::req_timeout(30) %>% httr2::req_perform(),
+    error = function(e) stop(paste("Rain archive request failed:", e$message), call. = FALSE)
+  )
+
+  if (httr2::resp_status(resp) != 200) {
+    stop(paste("Rain archive API returned status", httr2::resp_status(resp)), call. = FALSE)
+  }
+
+  body <- httr2::resp_body_json(resp)
+
+  if (is.null(body$daily) || is.null(body$daily$time)) {
+    stop("Unexpected Open-Meteo archive response: missing daily/time.", call. = FALSE)
+  }
+
+  data.frame(
+    date    = as.Date(unlist(body$daily$time)),
+    rain_mm = as.numeric(unlist(body$daily$precipitation_sum)),
+    stringsAsFactors = FALSE
+  )
+}
+
+# ============================================================================
+# GENERIC FETCH WRAPPER
+# ============================================================================
+
+#' Fetch daily precipitation via the configured provider
+#'
+#' @param latitude   Numeric.
+#' @param longitude  Numeric.
+#' @param start_date Date or "YYYY-MM-DD".
+#' @param end_date   Date or "YYYY-MM-DD".
+#' @param provider   Character. Currently only "open_meteo_archive".
+#' @return data.frame(date, rain_mm)
+rain_fetch <- function(latitude, longitude, start_date, end_date,
+                       provider = "open_meteo_archive") {
+  # --- ADD NEW PROVIDERS HERE ---
+  switch(provider,
+    open_meteo_archive = rain_fetch_open_meteo_archive(
+      latitude, longitude, start_date, end_date
+    ),
+    stop(paste("Unknown rain provider:", provider), call. = FALSE)
+  )
+}
+
+# ============================================================================
+# SPATIAL HELPERS
+# ============================================================================
+
+#' Snap a coordinate to the nearest ERA5 tile centre
+#'
+#' Open-Meteo ERA5 has 0.25° native resolution (~28 km). Snapping coordinates
+#' to this grid avoids redundant API calls for nearby fields.
+#'
+#' @param latitude   Numeric.
+#' @param longitude  Numeric.
+#' @param resolution Numeric. Grid resolution in degrees (default 0.25).
+#' @return Named list: tile_lat, tile_lon, tile_id
+rain_snap_to_tile <- function(latitude, longitude, resolution = 0.25) {
+  tile_lat <- round(latitude  / resolution) * resolution
+  tile_lon <- round(longitude / resolution) * resolution
+  list(
+    tile_lat = tile_lat,
+    tile_lon = tile_lon,
+    tile_id  = paste0(tile_lat, "_", tile_lon)
+  )
+}
+
+# ============================================================================
+# BATCH FETCH WITH TILE DEDUPLICATION
+# ============================================================================
+
+#' Fetch rain for multiple fields, calling the API once per unique 0.25° tile
+#'
+#' Fields that fall on the same ERA5 tile share one API call. This handles
+#' estates spread over 100+ km (e.g. Angata) without over-calling the API.
+#'
+#' @param centroids_df data.frame with columns: field_id, latitude, longitude
+#' @param start_date   Date or "YYYY-MM-DD". Start of fetch window.
+#' @param end_date     Date or "YYYY-MM-DD". End of fetch window.
+#' @param provider     Character. Passed to rain_fetch().
+#' @return data.frame with columns: field_id, date, rain_mm
+rain_fetch_for_fields <- function(centroids_df, start_date, end_date,
+                                  provider = "open_meteo_archive") {
+  required_cols <- c("field_id", "latitude", "longitude")
+  if (!all(required_cols %in% names(centroids_df))) {
+    stop(paste("centroids_df must have:", paste(required_cols, collapse = ", ")), call. = FALSE)
+  }
+
+  # Snap each field to its tile
+  centroids_tiled <- centroids_df %>%
+    dplyr::mutate(
+      tile_info = purrr::map2(latitude, longitude, rain_snap_to_tile),
+      tile_lat  = purrr::map_dbl(tile_info, "tile_lat"),
+      tile_lon  = purrr::map_dbl(tile_info, "tile_lon"),
+      tile_id   = purrr::map_chr(tile_info, "tile_id")
+    ) %>%
+    dplyr::select(-tile_info)
+
+  # Fetch once per unique tile
+  unique_tiles <- centroids_tiled %>%
+    dplyr::distinct(tile_id, tile_lat, tile_lon)
+
+  tile_rain <- unique_tiles %>%
+    dplyr::mutate(rain_data = purrr::pmap(
+      list(tile_lat, tile_lon),
+      function(lat, lon) {
+        tryCatch(
+          rain_fetch(lat, lon, start_date, end_date, provider = provider),
+          error = function(e) {
+            safe_log(paste0("Rain fetch failed for tile ", lat, "_", lon, ": ", e$message), "WARNING")
+            data.frame(date = as.Date(character(0)), rain_mm = numeric(0))
+          }
+        )
+      }
+    ))
+
+  # Join tile rain back to fields and unnest
+  centroids_tiled %>%
+    dplyr::left_join(tile_rain %>% dplyr::select(tile_id, rain_data), by = "tile_id") %>%
+    dplyr::select(field_id, rain_data) %>%
+    tidyr::unnest(rain_data)
+}
+
+# ============================================================================
+# PROCESSING
+# ============================================================================
+
+#' Join rain data to latest-season CI data by Date
+#'
+#' @param rain_df       data.frame(date, rain_mm) for the field.
+#' @param ci_season_data data.frame with at minimum: Date, DAH, week columns
+#'                       (latest season only, one row per date).
+#' @return data.frame with Date, DAH, week, rain_mm, cum_rain_mm columns,
+#'         ordered by Date. Only dates present in ci_season_data are returned.
+rain_join_to_ci <- function(rain_df, ci_season_data) {
+  if (is.null(rain_df) || nrow(rain_df) == 0) return(NULL)
+
+  ci_dates <- ci_season_data %>%
+    dplyr::distinct(Date, DAH, week) %>%
+    dplyr::arrange(Date)
+
+  joined <- ci_dates %>%
+    dplyr::left_join(
+      rain_df %>% dplyr::rename(Date = date),
+      by = "Date"
+    ) %>%
+    dplyr::mutate(
+      rain_mm     = tidyr::replace_na(rain_mm, 0),
+      cum_rain_mm = cumsum(rain_mm)
+    )
+
+  joined
+}
+
+# ============================================================================
+# PLOT OVERLAY — single plot (absolute or cumulative)
+# ============================================================================
+
+#' Add a rainfall overlay to a single-panel CI ggplot
+#'
+#' For absolute CI (mean_rolling_10_days): adds daily precipitation bars
+#' (steelblue, semi-transparent) scaled to the primary y-axis range, with a
+#' secondary right y-axis labelled in mm/day.
+#'
+#' For cumulative CI: adds a cumulative precipitation filled area scaled to the
+#' primary y-axis range, with a secondary right y-axis labelled in mm.
+#'
+#' The secondary axis is a linear transform of the primary axis (ggplot2
+#' requirement). The rain values are scaled so the maximum rain never exceeds
+#' y_max on screen. The right axis labels show the real mm values.
+#'
+#' @param g           ggplot object.
+#' @param rain_joined Result of rain_join_to_ci(). NULL returns g unchanged.
+#' @param ci_type     "mean_rolling_10_days" or "cumulative_CI".
+#' @param x_var       Character: "DAH", "week", or "Date".
+#' @param y_max       Numeric. Max value of the primary y-axis (used for scaling
+#'                    and for setting limits = c(0, y_max) on the primary axis).
+#' @return Modified ggplot object.
+rain_add_to_plot <- function(g, rain_joined, ci_type, x_var, y_max) {
+  if (is.null(rain_joined) || nrow(rain_joined) == 0) {
+    # Still enforce y limits even without rain
+    return(g + ggplot2::scale_y_continuous(limits = c(0, y_max)))
+  }
+
+  if (ci_type == "mean_rolling_10_days") {
+    # --- Daily precipitation bars ---
+    max_rain <- max(rain_joined$rain_mm, na.rm = TRUE)
+    if (is.na(max_rain) || max_rain == 0) {
+      return(g + ggplot2::scale_y_continuous(limits = c(0, y_max)))
+    }
+
+    scale_factor <- y_max / max_rain
+
+    bar_width <- switch(x_var,
+      "DAH"  = 1,
+      "week" = 0.5,
+      1      # Date: 1 day width handled by ggplot
+    )
+
+    g +
+      ggplot2::geom_col(
+        data        = rain_joined,
+        ggplot2::aes(x = .data[[x_var]], y = rain_mm * scale_factor),
+        fill        = "steelblue",
+        alpha       = 0.65,
+        width       = bar_width,
+        inherit.aes = FALSE
+      ) +
+      ggplot2::scale_y_continuous(
+        limits   = c(0, y_max),
+        sec.axis = ggplot2::sec_axis(
+          ~ . / scale_factor,
+          name = tr_key("lbl_rain_mm_day", "Precipitation (mm/day)")
+        )
+      )
+
+  } else if (ci_type == "cumulative_CI") {
+    # --- Cumulative precipitation area ---
+    max_cum_rain <- max(rain_joined$cum_rain_mm, na.rm = TRUE)
+    if (is.na(max_cum_rain) || max_cum_rain == 0) {
+      return(g)
+    }
+
+    scale_factor <- (y_max * 0.30) / max_cum_rain
+
+    g +
+      ggplot2::geom_area(
+        data        = rain_joined,
+        ggplot2::aes(x = .data[[x_var]], y = cum_rain_mm * scale_factor),
+        fill        = "steelblue",
+        alpha       = 0.15,
+        inherit.aes = FALSE
+      ) +
+      ggplot2::geom_line(
+        data        = rain_joined,
+        ggplot2::aes(x = .data[[x_var]], y = cum_rain_mm * scale_factor),
+        color       = "steelblue",
+        linewidth   = 0.8,
+        inherit.aes = FALSE
+      ) +
+      ggplot2::scale_y_continuous(
+        limits   = c(0, y_max),
+        sec.axis = ggplot2::sec_axis(
+          ~ . / scale_factor,
+          name = tr_key("lbl_cum_rain_mm", "Cumulative Rain (mm)")
+        )
+      )
+
+  } else {
+    g
+  }
+}
+
+# ============================================================================
+# PLOT OVERLAY — faceted "both" plot
+# ============================================================================
+
+#' Add rainfall overlays to a faceted (plot_type = "both") CI ggplot
+#'
+#' ggplot2 does not support per-facet secondary axes with free_y scales.
+#' This function adds the rain layers to the correct facets by setting
+#' ci_type_label to match the facet strip, but omits secondary axis labels
+#' to avoid misleading scale information across facets.
+#'
+#' Rolling mean facet receives: daily precipitation bars.
+#' Cumulative facet receives:   cumulative precipitation filled area.
+#'
+#' @param g                  ggplot object with facet_wrap(~ci_type_label).
+#' @param rain_joined        Result of rain_join_to_ci().
+#' @param rolling_mean_label Character. Must match the strip label for the
+#'                           rolling mean facet (from tr_key).
+#' @param cumulative_label   Character. Must match the strip label for the
+#'                           cumulative CI facet (from tr_key).
+#' @param x_var              Character: "DAH", "week", or "Date".
+#' @param y_max_abs          Numeric. Primary y-axis max for rolling mean facet (typically 7).
+#' @param y_max_cum          Numeric. Primary y-axis max for cumulative CI facet.
+#' @return Modified ggplot object.
+rain_add_to_faceted_plot <- function(g, rain_joined, rolling_mean_label, cumulative_label,
+                                     x_var, y_max_abs, y_max_cum) {
+  if (is.null(rain_joined) || nrow(rain_joined) == 0) return(g)
+
+  # --- Bars in rolling mean facet ---
+  max_rain <- max(rain_joined$rain_mm, na.rm = TRUE)
+  if (!is.na(max_rain) && max_rain > 0) {
+    scale_abs <- y_max_abs / max_rain
+    bar_width <- switch(x_var, "DAH" = 1, "week" = 0.5, 1)
+
+    rain_bars <- rain_joined %>%
+      dplyr::mutate(
+        ci_type_label = rolling_mean_label,
+        y_scaled      = rain_mm * scale_abs
+      )
+
+    g <- g +
+      ggplot2::geom_col(
+        data        = rain_bars,
+        ggplot2::aes(x = .data[[x_var]], y = y_scaled),
+        fill        = "steelblue",
+        alpha       = 0.65,
+        width       = bar_width,
+        inherit.aes = FALSE
+      )
+  }
+
+  # --- Area in cumulative facet ---
+  max_cum_rain <- max(rain_joined$cum_rain_mm, na.rm = TRUE)
+  if (!is.na(max_cum_rain) && max_cum_rain > 0) {
+    scale_cum <- (y_max_cum * 0.30) / max_cum_rain
+
+    rain_area <- rain_joined %>%
+      dplyr::mutate(
+        ci_type_label = cumulative_label,
+        y_scaled      = cum_rain_mm * scale_cum
+      )
+
+    g <- g +
+      ggplot2::geom_line(
+        data        = rain_area,
+        ggplot2::aes(x = .data[[x_var]], y = y_scaled),
+        color       = "steelblue",
+        linewidth   = 0.8,
+        inherit.aes = FALSE
+      )
+  }
+
+  g
+}

r_app/90_report_utils.R

@@ -13,8 +13,15 @@
 #'
 subchunkify <- function(g, fig_height=7, fig_width=5) {
   g_deparsed <- paste0(deparse(
-    function() {g}
-  ), collapse = '')
+    function() {
+      if (inherits(g, c("gtable", "grob", "gTree"))) {
+        grid::grid.newpage()
+        grid::grid.draw(g)
+      } else {
+        print(g)
+      }
+    }
+  ), collapse = '\n')
   
   sub_chunk <- paste0("
 `","``{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", dpi=300, dev='png', out.width='100%', echo=FALSE}",
@@ -394,9 +401,11 @@ ci_plot <- function(pivotName,
 #' @param show_benchmarks Whether to show historical benchmark lines (default: FALSE)
 #' @param estate_name Name of the estate for benchmark calculation (required if show_benchmarks = TRUE)
 #' @param benchmark_percentiles Vector of percentiles for benchmarks (default: c(10, 50, 90))
+#' @param rain_data data.frame(date, rain_mm) for this field's current season,
+#'   as returned by rain_fetch_for_fields() / rain_join_to_ci(). NULL disables overlay.
 #' @return NULL (adds output directly to R Markdown document)
 #'
-cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "absolute", facet_on = FALSE, x_unit = "days", colorblind_friendly = FALSE, show_benchmarks = FALSE, estate_name = NULL, benchmark_percentiles = c(10, 50, 90), benchmark_data = NULL) {
+cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "absolute", facet_on = FALSE, x_unit = "days", colorblind_friendly = FALSE, show_benchmarks = FALSE, estate_name = NULL, benchmark_percentiles = c(10, 50, 90), benchmark_data = NULL, rain_data = NULL) {
   # Input validation
   if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
     stop("pivotName is required")
@@ -524,6 +533,24 @@ cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "
                 legend.title = ggplot2::element_text(size = 8),
                 legend.text = ggplot2::element_text(size = 8)) +
           ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))
+
+        # Add benchmark lines to faceted plot (map DAH → Date per season)
+        if (!is.null(benchmark_data) && ci_type_filter %in% benchmark_data$ci_type) {
+          max_dah_clip <- max(plot_data$DAH, na.rm = TRUE) * 1.1
+          benchmark_facet <- benchmark_data %>%
+            dplyr::filter(ci_type == ci_type_filter, DAH <= max_dah_clip) %>%
+            dplyr::cross_join(
+              date_preparation_perfect_pivot %>%
+                dplyr::filter(season %in% unique_seasons) %>%
+                dplyr::select(season, min_date)
+            ) %>%
+            dplyr::mutate(Date = min_date + DAH - 1)
+          g <- g + ggplot2::geom_smooth(
+            data = benchmark_facet,
+            ggplot2::aes(x = Date, y = benchmark_value, group = factor(percentile)),
+            color = "gray70", linewidth = 0.5, se = FALSE, inherit.aes = FALSE, fullrange = FALSE
+          )
+        }
       } else {
         # Choose color palette based on colorblind_friendly flag
         color_scale <- if (colorblind_friendly) {
@@ -638,11 +665,25 @@ cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "
           ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))
       }
       
-      # Add y-axis limits for absolute CI (10-day rolling mean) to fix scale at 0-7
-      if (ci_type_filter == "mean_rolling_10_days") {
-        g <- g + ggplot2::ylim(0, 7)
+      # Y-axis limits + optional rain overlay (single-panel plots)
+      y_max <- if (ci_type_filter == "mean_rolling_10_days") {
+        7
+      } else {
+        max(plot_data$ci_value, na.rm = TRUE) * 1.05
       }
-      
+
+      if (!is.null(rain_data)) {
+        latest_ci_dates <- plot_data %>%
+          dplyr::filter(is_latest) %>%
+          dplyr::distinct(Date, DAH, week)
+        rain_joined <- rain_join_to_ci(rain_data, latest_ci_dates)
+        g <- rain_add_to_plot(g, rain_joined, ci_type_filter, x_var, y_max)
+      } else {
+        if (ci_type_filter == "mean_rolling_10_days") {
+          g <- g + ggplot2::ylim(0, y_max)
+        }
+      }
+
       return(g)
     }
     
@@ -654,174 +695,20 @@ cum_ci_plot <- function(pivotName, ci_quadrant_data = CI_quadrant, plot_type = "
       g <- create_plot("cumulative_CI", cumulative_label, "")
       subchunkify(g, 2.8, 10)
     } else if (plot_type == "both") {
-      # Create faceted plot with both CI types using pivot_longer approach
-      plot_data_both <- data_ci3 %>%
-        dplyr::filter(season %in% unique_seasons) %>%
-        dplyr::mutate(
-          ci_type_label = factor(case_when(
-            ci_type == "mean_rolling_10_days" ~ rolling_mean_label,
-            ci_type == "cumulative_CI" ~ cumulative_label,
-            TRUE ~ ci_type
-          ), levels = c(rolling_mean_label, cumulative_label)),
-          is_latest = season == latest_season  # Flag for latest season
-        )
-      
-      # Determine x-axis variable based on x_unit parameter
-      x_var <- if (x_unit == "days") {
-        if (facet_on) "Date" else "DAH"
-      } else {
-        "week"
-      }
-      
-      x_label <- switch(x_unit,
-                       "days" = if (facet_on) tr_key("lbl_date", "Date") else tr_key("lbl_age_of_crop_days", "Age of Crop (Days)"),
-                       "weeks" = tr_key("lbl_week_number", "Week Number"))
-
-      # Choose color palette based on colorblind_friendly flag
-      color_scale <- if (colorblind_friendly) {
-        ggplot2::scale_color_brewer(type = "qual", palette = "Set2")
-      } else {
-        ggplot2::scale_color_discrete()
-      }
-                 
-      # Calculate dynamic max values for breaks
-      max_dah_both <- max(plot_data_both$DAH, na.rm = TRUE) + 20
-      max_week_both <- max(as.numeric(plot_data_both$week), na.rm = TRUE) + ceiling(20 / 7)
-
-      # Pre-evaluate translated title here (not inside labs()) so {pivotName} resolves correctly
+      # Build each panel independently so each gets its own secondary rain y-axis.
+      # (facet_wrap + free_y does not support per-facet sec.axis in ggplot2.)
       both_plot_title <- tr_key("lbl_ci_analysis_title", "CI Analysis for Field {pivotName}")
 
-      # Create the faceted plot
-      g_both <- ggplot2::ggplot(data = plot_data_both) +
-        # Add benchmark lines first (behind season lines)
-        {
-          if (!is.null(benchmark_data)) {
-            # Clip benchmark to max DAH of plotted seasons + 10% buffer
-            max_dah_clip <- max(plot_data_both$DAH, na.rm = TRUE) * 1.1
-            benchmark_subset <- benchmark_data %>%
-              dplyr::filter(DAH <= max_dah_clip) %>%
-              dplyr::mutate(
-                benchmark_x = if (x_var == "DAH") {
-                  DAH
-                } else if (x_var == "week") {
-                  DAH / 7
-                } else {
-                  DAH
-                },
-                ci_type_label = factor(case_when(
-                  ci_type == "value" ~ rolling_mean_label,
-                  ci_type == "cumulative_CI" ~ cumulative_label,
-                  TRUE ~ ci_type
-                ), levels = c(rolling_mean_label, cumulative_label))
-              )
-            ggplot2::geom_smooth(
-              data = benchmark_subset,
-              ggplot2::aes(
-                x = .data[["benchmark_x"]],
-                y = .data[["benchmark_value"]],
-                group = factor(.data[["percentile"]])
-              ),
-              color = "gray70", linewidth = 0.5, se = FALSE, inherit.aes = FALSE, fullrange = FALSE
-            )
-          }
-        } +
-        ggplot2::facet_wrap(~ci_type_label, scales = "free_y") +
-        # Plot older seasons with lighter lines
-        ggplot2::geom_line(
-          data = plot_data_both %>% dplyr::filter(!is_latest),
-          ggplot2::aes(
-            x = .data[[x_var]],
-            y = .data[["ci_value"]],
-            col = .data[["season"]],
-            group = .data[["season"]]
-          ),
-          linewidth = 0.7, alpha = 0.4
-        ) +
-        # Plot latest season with thicker, more prominent line
-        ggplot2::geom_line(
-          data = plot_data_both %>% dplyr::filter(is_latest),
-          ggplot2::aes(
-            x = .data[[x_var]],
-            y = .data[["ci_value"]],
-            col = .data[["season"]],
-            group = .data[["season"]]
-          ),
-          linewidth = 1.5, alpha = 1
-        ) +
-        ggplot2::labs(title = both_plot_title,
-             color = tr_key("lbl_season", "Season"),
-             y = tr_key("lbl_ci_value", "CI Value"),
-             x = x_label) +
-        color_scale +
-        {
-          if (x_var == "DAH") {
-            # Dynamic breaks based on actual data range
-            dah_breaks <- scales::pretty_breaks(n = 5)(c(0, max_dah_both))
-            # Month breaks: in secondary axis scale (month numbers 0, 1, 2, ...)
-            n_months <- ceiling(max_dah_both / 30.44)
-            month_breaks <- seq(0, n_months, by = 1)
-            ggplot2::scale_x_continuous(
-              breaks = dah_breaks,
-              sec.axis = ggplot2::sec_axis(
-                ~ . / 30.44,
-                name = tr_key("lbl_age_in_months", "Age in Months"),
-                breaks = month_breaks,
-                labels = function(x) as.integer(x)  # Show all month labels
-              )
-            )
-          } else if (x_var == "week") {
-            # Dynamic breaks based on actual data range
-            week_breaks <- scales::pretty_breaks(n = 5)(c(0, max_week_both))
-            # Month breaks: in secondary axis scale (month numbers 0, 1, 2, ...)
-            n_months <- ceiling(max_week_both / 4.348)
-            month_breaks <- seq(0, n_months, by = 1)
-            ggplot2::scale_x_continuous(
-              breaks = week_breaks,
-              sec.axis = ggplot2::sec_axis(
-                ~ . / 4.348,
-                name = tr_key("lbl_age_in_months", "Age in Months"),
-                breaks = month_breaks,
-                labels = function(x) as.integer(x)  # Show all month labels
-              )
-            )
-          } else if (x_var == "Date") {
-            ggplot2::scale_x_date(breaks = "1 month", date_labels = "%b-%Y", sec.axis = ggplot2::sec_axis(~ ., name = tr_key("lbl_age_in_months", "Age in Months"), breaks = scales::breaks_pretty()))
-          }
-        } +
-        ggplot2::theme_minimal() +
-          ggplot2::theme(axis.text.x = ggplot2::element_text(hjust = 0.5),
-            axis.text.x.top = ggplot2::element_text(hjust = 0.5),
-            axis.title.x.top = ggplot2::element_text(size = 8),
-            legend.justification = c(1, 0),
-            legend.position = "inside",
-            legend.position.inside = c(1, 0),
-            legend.title = ggplot2::element_text(size = 8),
-            legend.text = ggplot2::element_text(size = 8)) +
-        ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = TRUE))
-      
-      # For the rolling mean data, we want to set reasonable y-axis limits
-      # Since we're using free_y scales, each facet will have its own y-axis
-      # The rolling mean will automatically scale to its data range, 
-      # but we can ensure it shows the 0-7 context by adding invisible points
-      
-      # Add invisible points to set the y-axis range for rolling mean facet
-      dummy_data <- data.frame(
-        ci_type_label = rolling_mean_label,
-        ci_value = c(0, 7),
-        stringsAsFactors = FALSE
-      )
-      dummy_data[[x_var]] <- range(plot_data_both[[x_var]], na.rm = TRUE)
-      dummy_data[["season"]] <- factor("dummy", levels = levels(plot_data_both[["season"]]))
-      
-      g_both <- g_both + 
-        ggplot2::geom_point(
-          data = dummy_data,
-          ggplot2::aes(x = .data[[x_var]], y = .data[["ci_value"]]),
-          alpha = 0, size = 0
-        )  # Invisible points to set scale
-      
-      # Display the combined faceted plot
-      subchunkify(g_both, 2.8, 10)
+      g_abs <- create_plot("mean_rolling_10_days", rolling_mean_label, "") +
+        ggplot2::labs(title = rolling_mean_label) +
+        ggplot2::theme(legend.position = "none")
+
+      g_cum <- create_plot("cumulative_CI", cumulative_label, "") +
+        ggplot2::labs(title = cumulative_label)
+
+      combined <- gridExtra::arrangeGrob(g_abs, g_cum, ncol = 2, top = both_plot_title)
+
+      subchunkify(combined, 2.8, 10)
     }
     
   }, error = function(e) {
@@ -962,9 +849,9 @@ compute_ci_benchmarks <- function(ci_quadrant_data, estate_name, percentiles = c
     # Prepare data for both CI types
     data_prepared <- data_filtered %>%
       dplyr::ungroup() %>%  # Ensure no existing groupings
-      dplyr::select(DAH, value, cumulative_CI, season) %>%
+      dplyr::select(DAH, mean_rolling_10_days = value, cumulative_CI, season) %>%
       tidyr::pivot_longer(
-        cols = c("value", "cumulative_CI"),
+        cols = c("mean_rolling_10_days", "cumulative_CI"),
         names_to = "ci_type",
         values_to = "ci_value"
       ) %>%

r_app/91_CI_report_with_kpis_cane_supply.Rmd

@@ -544,9 +544,11 @@ if (exists("summary_data") && !is.null(summary_data) && "field_analysis" %in% na
 tryCatch({
   # Use per-field field_analysis data from RDS (already loaded in load_kpi_data chunk)
   if (exists("summary_data") && !is.null(summary_data) && "field_analysis" %in% names(summary_data)) {
-    analysis_data <- summary_data$field_analysis %>% 
-      select(Field_id, Status_Alert) %>%
-      rename(Status_trigger = Status_Alert)  # Rename to Status_trigger for compatibility with hexbin logic
+    fa <- summary_data$field_analysis
+    status_col <- intersect(c("Status_Alert", "Status_trigger"), names(fa))[1]
+    analysis_data <- fa %>%
+      select(Field_id, all_of(status_col)) %>%
+      rename(Status_trigger = all_of(status_col))
   } else {
     analysis_data <- tibble(Field_id = character(), Status_trigger = character())
   }
@@ -619,10 +621,9 @@ tryCatch({
     filter(Status_trigger != "harvest_ready" | is.na(Status_trigger))
 
   # Generate breaks for color gradients
-  breaks_vec <- c(0, 5, 10, 15, 20, 30, 35)
+  breaks_vec <- c(5, 10, 15, 20, 30, 35)
   labels_vec <- as.character(breaks_vec)
   labels_vec[length(labels_vec)] <- ">30"
-  labels_vec[1] <- "0.1"
   
   # Calculate data bounds for coordinate limits (prevents basemap scale conflicts)
   # Use actual data bounds without dummy points to avoid column mismatch
@@ -635,36 +636,45 @@ tryCatch({
     ceiling(max(points_processed$Y, na.rm = TRUE) * 20) / 20
   )
 
+  # ggplot2's StatBinhex computes the hex grid origin as c(min(x), min(y)) of each layer's
+  # own data. Feeding different subsets (points_ready / points_not_ready) to the two layers
+  # gives different origins → misaligned hexes. Fix: pass points_processed to BOTH layers.
+  # For the coloured layer, use ready_area = area_value for harvest-ready fields, 0 elsewhere.
+  # The scale censors 0 → NA → transparent, so non-ready hexes are invisible in that layer.
+  points_processed <- points_processed %>%
+    mutate(ready_area = ifelse(Status_trigger == "harvest_ready", area_value, 0))
+
   # Create hexbin map with enhanced aesthetics, basemap, and proper legend
   ggplot() +
     # OpenStreetMap basemap (zoom=11 appropriate for agricultural fields)
     ggspatial::annotation_map_tile(type = "osm", zoom = 11, progress = "none", alpha = 0.5) +
-    # Hexbin for NOT ready fields (light background)
+    # Hexbin background: ALL fields as white hexes (anchors the hex grid to the full extent)
     geom_hex(
-      data = points_not_ready, 
-      aes(x = X, y = Y, weight = area_value, alpha = tr_key("hexbin_not_ready")), 
+      data = points_processed,
+      aes(x = X, y = Y, weight = area_value, alpha = tr_key("hexbin_not_ready")),
       binwidth = c(0.012, 0.012),
       fill = "#ffffff",
       colour = "#0000009a",
       linewidth = 0.1
     ) +
-    # Hexbin for READY fields (colored gradient)
+    # Hexbin overlay: same grid (same data), coloured only where ready_area > 0
     geom_hex(
-      data = points_ready, 
-      aes(x = X, y = Y, weight = area_value), 
+      data = points_processed,
+      aes(x = X, y = Y, weight = ready_area),
       binwidth = c(0.012, 0.012),
       alpha = 0.9,
       colour = "#0000009a",
       linewidth = 0.1
     ) +
-    # Color gradient scale for area
+    # Color gradient scale — values of 0 (non-ready hexes) are censored → NA → transparent
     scale_fill_viridis_b(
       option = "viridis",
       direction = -1,
       breaks = breaks_vec,
       labels = labels_vec,
-      limits = c(0, 35),
-      oob = scales::squish,
+      limits = c(0.001, 35),
+      oob = scales::oob_censor,
+      na.value = "transparent",
       name = tr_key("hexbin_legend_acres")
     ) +
     # Alpha scale for "not ready" status indication
@@ -722,7 +732,10 @@ if (exists("summary_data") && !is.null(summary_data) && "field_analysis" %in% na
   if (is.null(summary_data$field_analysis_summary) || !("field_analysis_summary" %in% names(summary_data)) || 
       !is.data.frame(summary_data$field_analysis_summary)) {
     
-    # Create summary by aggregating by Status_Alert and Phase categories
+    # Detect which status column is present (Status_Alert from cane_supply, Status_trigger from others)
+    status_col <- intersect(c("Status_Alert", "Status_trigger"), names(field_analysis_df))[1]
+
+    # Create summary by aggregating by status and Phase categories
     phase_summary <- field_analysis_df %>%
       filter(!is.na(Phase)) %>%
       group_by(Phase) %>%
@@ -736,21 +749,21 @@ if (exists("summary_data") && !is.null(summary_data) && "field_analysis" %in% na
     
     # Create Status trigger summary - includes both active alerts and "No active triggers"
     trigger_summary <- tryCatch({
-      # Active alerts (fields with non-NA Status_Alert)
+      # Active alerts (fields with non-NA status)
       active_alerts <- field_analysis_df %>%
-        filter(!is.na(Status_Alert), Status_Alert != "") %>%
-        group_by(Status_Alert) %>%
+        filter(!is.na(.data[[status_col]]), .data[[status_col]] != "") %>%
+        group_by(across(all_of(status_col))) %>%
         summarise(
           Acreage = sum(Acreage, na.rm = TRUE),
           Field_count = n_distinct(Field_id),
           .groups = "drop"
         ) %>%
-        mutate(Category = Status_Alert) %>%
+        mutate(Category = .data[[status_col]]) %>%
         select(Category, Acreage, Field_count)
-      
-      # No active triggers (fields with NA Status_Alert)
+
+      # No active triggers (fields with NA status)
       no_alerts <- field_analysis_df %>%
-        filter(is.na(Status_Alert) | Status_Alert == "") %>%
+        filter(is.na(.data[[status_col]]) | .data[[status_col]] == "") %>%
         summarise(
           Acreage = sum(Acreage, na.rm = TRUE),
           Field_count = n_distinct(Field_id),

r_app/MANUAL_PIPELINE_RUNNER.R

@@ -438,8 +438,8 @@
 #   rmarkdown::render(
 rmarkdown::render(
   "r_app/90_CI_report_with_kpis_agronomic_support.Rmd",
-  params = list(data_dir = "aura", report_date = as.Date("2026-02-18"), language = "en" ),
-  output_file = "SmartCane_Report_agronomic_support_aura_2026-02-18_en_test.docx",
+  params = list(data_dir = "aura", report_date = as.Date("2026-03-23"), language = "en" ),
+  output_file = "SmartCane_Report_agronomic_support_aura_2026-03-23_en.docx",
   output_dir = "laravel_app/storage/app/aura/reports"
 )
 
@@ -456,8 +456,8 @@ rmarkdown::render(
 #   rmarkdown::render(
 rmarkdown::render(
   "r_app/91_CI_report_with_kpis_cane_supply.Rmd",
-  params = list(data_dir = "angata", report_date = as.Date("2026-02-23")),
-  output_file = "SmartCane_Report_cane_supply_angata_2026-02-23_en.docx",
+  params = list(data_dir = "angata", report_date = as.Date("2026-03-17")),
+  output_file = "SmartCane_Report_cane_supply_angata_2026-03-17_en.docx",
   output_dir = "laravel_app/storage/app/angata/reports"
 )
 #

r_app/parameters_project.R

@@ -57,7 +57,7 @@ AREA_UNIT_PREFERENCE <- tolower(Sys.getenv("AREA_UNIT", unset = "hectare"))
 # Validate area unit value
 if (!AREA_UNIT_PREFERENCE %in% c("hectare", "acre")) {
   warning(paste0("Invalid AREA_UNIT env var: '", AREA_UNIT_PREFERENCE, "'. Using 'hectare'."))
-  AREA_UNIT_PREFERENCE <- "hectare"
+  AREA_UNIT_PREFERENCE <- "acre"
 }
 
 #' Get area unit label for display

r_app/translations/translations_90.json

@@ -512,6 +512,16 @@
       "es-mx": "**Total de parcelas analizadas:** {total_fields}",
       "pt-br": "**Total de campos analisados:** {total_fields}"
     },
+    "tot_area_analyzed": {
+      "en": "**Total Area Analyzed:** {round(total_area, 1)} {unit_label}",
+      "es-mx": "**Área total analizada:** {round(total_area, 1)} {unit_label}",
+      "pt-br": "**Área total analisada:** {round(total_area, 1)} {unit_label}"
+    },
+    "Area": {
+      "en": "Area",
+      "es-mx": "Área",
+      "pt-br": "Área"
+    },
     "Medium": {
       "en": "Medium",
       "es-mx": "Medio",

r_app/translations/translations_91.json

@@ -2,541 +2,541 @@
   "main_translations": {
     "cover_title": {
       "en": "Satellite Based Field Reporting",
-      "es-mx": "",
+      "es-mx": "Reportes de Campo Basados en Satélite",
       "pt-br": "",
       "sw": ""
     },
     "cover_subtitle": {
       "en": "Cane Supply Office\n\nChlorophyll Index (CI) Monitoring Report — {toupper(params$data_dir)} Estate (Week {if (!is.null(params$week)) params$week else format(as.Date(params$report_date), '%V')}, {format(as.Date(params$report_date), '%Y')})",
-      "es-mx": "",
+      "es-mx": "Oficina de Abastecimiento de Caña\n\nReporte de Monitoreo del Índice de Clorofila (IC) — Finca {toupper(params$data_dir)} (Semana {if (!is.null(params$week)) params$week else format(as.Date(params$report_date), '%V')}, {format(as.Date(params$report_date), '%Y')})",
       "pt-br": "",
       "sw": ""
     },
     "report_summary_heading": {
       "en": "## Report Generated",
-      "es-mx": "",
+      "es-mx": "## Reporte Generado",
       "pt-br": "",
       "sw": ""
     },
     "report_farm_location": {
       "en": "**Farm Location:**",
-      "es-mx": "",
+      "es-mx": "**Ubicación de la Finca:**",
       "pt-br": "",
       "sw": ""
     },
     "report_period_label": {
       "en": "**Report Period:** Week {current_week} of {year}",
-      "es-mx": "",
+      "es-mx": "**Período del Reporte:** Semana {current_week} de {year}",
       "pt-br": "",
       "sw": ""
     },
     "report_generated_on": {
       "en": "**Report Generated on:**",
-      "es-mx": "",
+      "es-mx": "**Reporte Generado el:**",
       "pt-br": "",
       "sw": ""
     },
     "report_farm_size": {
       "en": "**Farm Size Included in Analysis:** {formatC(total_acreage, format='f', digits=1)} acres",
-      "es-mx": "",
+      "es-mx": "**Tamaño de Finca Incluido en el Análisis:** {formatC(total_acreage, format='f', digits=1)} acres",
       "pt-br": "",
       "sw": ""
     },
     "report_data_source": {
       "en": "**Data Source:** Planet Labs Satellite Imagery",
-      "es-mx": "",
+      "es-mx": "**Fuente de Datos:** Imágenes Satelitales de Planet Labs",
       "pt-br": "",
       "sw": ""
     },
     "report_analysis_type": {
       "en": "**Analysis Type:** Chlorophyll Index (CI) Monitoring",
-      "es-mx": "",
+      "es-mx": "**Tipo de Análisis:** Monitoreo del Índice de Clorofila (IC)",
       "pt-br": "",
       "sw": ""
     },
     "key_insights": {
       "en": "## Key Insights",
-      "es-mx": "",
+      "es-mx": "## Hallazgos Clave",
       "pt-br": "",
       "sw": ""
     },
     "insight_excellent_unif": {
       "en": "- {excellent_pct}% of fields have excellent uniformity (CV < 0.08)",
-      "es-mx": "",
+      "es-mx": "- {excellent_pct}% de los campos tienen uniformidad excelente (CV < 0.08)",
       "pt-br": "",
       "sw": ""
     },
     "insight_good_unif": {
       "en": "- {good_pct}% of fields have good uniformity (CV < 0.15)",
-      "es-mx": "",
+      "es-mx": "- {good_pct}% de los campos tienen buena uniformidad (CV < 0.15)",
       "pt-br": "",
       "sw": ""
     },
     "insight_improving": {
       "en": "- {round(improving_acreage, 1)} acres ({improving_pct}%) of farm area is improving week-over-week",
-      "es-mx": "",
+      "es-mx": "- {round(improving_acreage, 1)} acres ({improving_pct}%) del área de la finca está mejorando semana a semana",
       "pt-br": "",
       "sw": ""
     },
     "insight_declining": {
       "en": "- {round(declining_acreage, 1)} acres ({declining_pct}%) of farm area is declining week-over-week",
-      "es-mx": "",
+      "es-mx": "- {round(declining_acreage, 1)} acres ({declining_pct}%) del área de la finca está en declive semana a semana",
       "pt-br": "",
       "sw": ""
     },
     "kpi_na_insights": {
       "en": "KPI data not available for key insights.",
-      "es-mx": "",
+      "es-mx": "Datos de KPI no disponibles para hallazgos clave.",
       "pt-br": "",
       "sw": ""
     },
     "report_structure_heading": {
       "en": "## Report Structure",
-      "es-mx": "",
+      "es-mx": "## Estructura del Reporte",
       "pt-br": "",
       "sw": ""
     },
     "report_structure_body": {
       "en": "**Section 1:** Cane supply zone analyses, summaries and Key Performance Indicators (KPIs)\n**Section 2:** Explanation of the report, definitions, methodology, and CSV export structure",
-      "es-mx": "",
+      "es-mx": "**Sección 1:** Análisis de zonas de abastecimiento de caña, resúmenes e Indicadores Clave de Desempeño (KPIs)\n**Sección 2:** Explicación del reporte, definiciones, metodología y estructura de exportación CSV",
       "pt-br": "",
       "sw": ""
     },
     "section_i": {
       "en": "# Section 1: Farm-wide Analyses and KPIs",
-      "es-mx": "",
+      "es-mx": "# Sección 1: Análisis General de la Finca y KPIs",
       "pt-br": "",
       "sw": ""
     },
     "section_1_1": {
       "en": "## 1.1 Overview of cane supply area, showing zones with number of acres being harvest ready",
-      "es-mx": "",
+      "es-mx": "## 1.1 Panorama del área de abastecimiento de caña, mostrando zonas con número de acres listos para cosecha",
       "pt-br": "",
       "sw": ""
     },
     "section_1_2": {
       "en": "## 1.2 Key Performance Indicators",
-      "es-mx": "",
+      "es-mx": "## 1.2 Indicadores Clave de Desempeño",
       "pt-br": "",
       "sw": ""
     },
     "kpi_empty": {
       "en": "KPI summary data available but is empty/invalid.",
-      "es-mx": "",
+      "es-mx": "Datos de resumen de KPI disponibles pero vacíos/inválidos.",
       "pt-br": "",
       "sw": ""
     },
     "kpi_unavailable": {
       "en": "KPI summary data not available.",
-      "es-mx": "",
+      "es-mx": "Datos de resumen de KPI no disponibles.",
       "pt-br": "",
       "sw": ""
     },
     "metadata": {
       "en": "## Report Metadata",
-      "es-mx": "",
+      "es-mx": "## Metadatos del Reporte",
       "pt-br": "",
       "sw": ""
     },
     "disclaimer": {
       "en": "*This report was automatically generated by the SmartCane monitoring system. For questions or additional analysis, please contact the technical team at info@smartcane.ag.*",
-      "es-mx": "",
+      "es-mx": "*Este reporte fue generado automáticamente por el sistema de monitoreo SmartCane. Para preguntas o análisis adicionales, por favor contacte al equipo técnico en info@smartcane.ag.*",
       "pt-br": "",
       "sw": ""
     },
     "section_2_heading": {
       "en": "# Section 2: Support Document for weekly SmartCane data package.",
-      "es-mx": "",
+      "es-mx": "# Sección 2: Documento de Soporte para el paquete de datos semanal de SmartCane.",
       "pt-br": "",
       "sw": ""
     },
     "about_doc_heading": {
       "en": "## 1. About This Document",
-      "es-mx": "",
+      "es-mx": "## 1. Acerca de Este Documento",
       "pt-br": "",
       "sw": ""
     },
     "about_doc_body": {
       "en": "This document is the support document to the SmartCane data file. It includes the definitions, explanatory calculations and suggestions for interpretations of the data as provided. For additional questions please feel free to contact SmartCane support, through your contact person, or via info@smartcane.ag.",
-      "es-mx": "",
+      "es-mx": "Este documento es el documento de soporte del archivo de datos SmartCane. Incluye las definiciones, cálculos explicativos y sugerencias para la interpretación de los datos proporcionados. Para preguntas adicionales, no dude en contactar al soporte de SmartCane, a través de su persona de contacto o vía info@smartcane.ag.",
       "pt-br": "",
       "sw": ""
     },
     "about_data_heading": {
       "en": "## 2. About the Data File",
-      "es-mx": "",
+      "es-mx": "## 2. Acerca del Archivo de Datos",
       "pt-br": "",
       "sw": ""
     },
     "about_data_body": {
       "en": "The data file is automatically populated based on normalized and indexed remote sensing images of provided polygons. Specific SmartCane algorithms provide tailored calculation results developed to support the sugarcane operations by:",
-      "es-mx": "",
+      "es-mx": "El archivo de datos se genera automáticamente a partir de imágenes de teledetección normalizadas e indexadas de los polígonos proporcionados. Los algoritmos específicos de SmartCane proporcionan resultados de cálculo personalizados desarrollados para apoyar las operaciones de caña de azúcar mediante:",
       "pt-br": "",
       "sw": ""
     },
     "about_data_bullet_harvest": {
       "en": "Supporting harvest planning mill-field logistics to ensure optimal tonnage and sucrose levels",
-      "es-mx": "",
+      "es-mx": "Apoyo a la planificación de cosecha y logística ingenio-campo para asegurar niveles óptimos de tonelaje y sacarosa",
       "pt-br": "",
       "sw": ""
     },
     "about_data_bullet_monitoring": {
       "en": "Monitoring of the crop growth rates on the farm, providing evidence of performance",
-      "es-mx": "",
+      "es-mx": "Monitoreo de las tasas de crecimiento del cultivo en la finca, proporcionando evidencia de rendimiento",
       "pt-br": "",
       "sw": ""
     },
     "about_data_bullet_identifying": {
       "en": "Identifying growth-related issues that are in need of attention",
-      "es-mx": "",
+      "es-mx": "Identificación de problemas relacionados con el crecimiento que requieren atención",
       "pt-br": "",
       "sw": ""
     },
     "about_data_bullet_enabling": {
       "en": "Enabling timely actions to minimize negative impact",
-      "es-mx": "",
+      "es-mx": "Habilitación de acciones oportunas para minimizar el impacto negativo",
       "pt-br": "",
       "sw": ""
     },
     "about_data_key_features": {
       "en": "Key Features of the data file: - High-resolution satellite imagery analysis - Week-over-week change detection - Individual field performance metrics - Actionable insights for crop management.",
-      "es-mx": "",
+      "es-mx": "Características principales del archivo de datos: - Análisis de imágenes satelitales de alta resolución - Detección de cambios semana a semana - Métricas de rendimiento por campo individual - Información accionable para el manejo del cultivo.",
       "pt-br": "",
       "sw": ""
     },
     "ci_section_heading": {
       "en": "#### *What is the Chlorophyll Index (CI)?*",
-      "es-mx": "",
+      "es-mx": "#### *¿Qué es el Índice de Clorofila (IC)?*",
       "pt-br": "",
       "sw": ""
     },
     "ci_intro_body": {
       "en": "The Chlorophyll Index (CI) is a vegetation index that measures the relative amount of chlorophyll in plant leaves. Chlorophyll is the green pigment responsible for photosynthesis in plants. Higher CI values indicate:",
-      "es-mx": "",
+      "es-mx": "El Índice de Clorofila (IC) es un índice de vegetación que mide la cantidad relativa de clorofila en las hojas de las plantas. La clorofila es el pigmento verde responsable de la fotosíntesis en las plantas. Valores más altos de IC indican:",
       "pt-br": "",
       "sw": ""
     },
     "ci_bullet_photosynthetic": {
       "en": "Greater photosynthetic activity",
-      "es-mx": "",
+      "es-mx": "Mayor actividad fotosintética",
       "pt-br": "",
       "sw": ""
     },
     "ci_bullet_healthy": {
       "en": "Healthier plant tissue",
-      "es-mx": "",
+      "es-mx": "Tejido vegetal más saludable",
       "pt-br": "",
       "sw": ""
     },
     "ci_bullet_nitrogen": {
       "en": "Better nitrogen uptake",
-      "es-mx": "",
+      "es-mx": "Mejor absorción de nitrógeno",
       "pt-br": "",
       "sw": ""
     },
     "ci_bullet_vigorous": {
       "en": "More vigorous crop growth",
-      "es-mx": "",
+      "es-mx": "Crecimiento más vigoroso del cultivo",
       "pt-br": "",
       "sw": ""
     },
     "ci_range_body": {
       "en": "CI values typically range from 0 (bare soil or severely stressed vegetation) to 7+ (very healthy, dense vegetation). For sugarcane, values between 3-7 generally indicate good crop health, depending on the growth stage.",
-      "es-mx": "",
+      "es-mx": "Los valores de IC generalmente varían de 0 (suelo descubierto o vegetación severamente estresada) a 7+ (vegetación muy saludable y densa). Para la caña de azúcar, valores entre 3-7 generalmente indican buena salud del cultivo, dependiendo de la etapa de crecimiento.",
       "pt-br": "",
       "sw": ""
     },
     "data_structure_heading": {
       "en": "### Data File Structure and Columns",
-      "es-mx": "",
+      "es-mx": "### Estructura del Archivo de Datos y Columnas",
       "pt-br": "",
       "sw": ""
     },
     "data_structure_intro": {
       "en": "The data file is organized in rows, one row per agricultural field (polygon), and columns, providing field data, actual measurements, calculation results and descriptions. The data file can be directly integrated with existing farm management systems for further analysis. Each column is described hereunder:",
-      "es-mx": "",
+      "es-mx": "El archivo de datos está organizado en filas, una fila por campo agrícola (polígono), y columnas, proporcionando datos del campo, mediciones reales, resultados de cálculos y descripciones. El archivo de datos puede integrarse directamente con los sistemas de gestión agrícola existentes para análisis adicional. Cada columna se describe a continuación:",
       "pt-br": "",
       "sw": ""
     },
     "col_field_id_desc": {
       "en": "Unique identifier for a cane field combining field name and sub-field number. This can be the same as Field_Name but is also helpful in keeping track of cane fields should they change, split or merge.",
-      "es-mx": "",
+      "es-mx": "Identificador único para un campo de caña que combina el nombre del campo y el número de subcampo. Puede ser igual que Field_Name pero también es útil para dar seguimiento a los campos de caña en caso de cambios, divisiones o fusiones.",
       "pt-br": "",
       "sw": ""
     },
     "col_farm_section_desc": {
       "en": "Sub-area or section name",
-      "es-mx": "",
+      "es-mx": "Nombre de subárea o sección",
       "pt-br": "",
       "sw": ""
     },
     "col_field_name_desc": {
       "en": "Client name or label assigned to a cane field.",
-      "es-mx": "",
+      "es-mx": "Nombre o etiqueta del cliente asignado a un campo de caña.",
       "pt-br": "",
       "sw": ""
     },
     "col_acreage_desc": {
       "en": "Field size in acres",
-      "es-mx": "",
+      "es-mx": "Tamaño del campo en acres",
       "pt-br": "",
       "sw": ""
     },
     "col_status_trigger_desc": {
       "en": "Shows changes in crop status worth alerting. More detailed explanation of the possible alerts is written down under key concepts.",
-      "es-mx": "",
+      "es-mx": "Muestra cambios en el estado del cultivo que ameritan alerta. Una explicación más detallada de las alertas posibles se encuentra en conceptos clave.",
       "pt-br": "",
       "sw": ""
     },
     "col_last_harvest_desc": {
       "en": "Date of most recent harvest as per satellite detection algorithm / or manual entry",
-      "es-mx": "",
+      "es-mx": "Fecha de la cosecha más reciente según el algoritmo de detección satelital / o ingreso manual",
       "pt-br": "",
       "sw": ""
     },
     "col_age_week_desc": {
       "en": "Time elapsed since planting/harvest in weeks; used to predict expected growth phases. Reflects planting/harvest date.",
-      "es-mx": "",
+      "es-mx": "Tiempo transcurrido desde la siembra/cosecha en semanas; utilizado para predecir fases de crecimiento esperadas. Refleja la fecha de siembra/cosecha.",
       "pt-br": "",
       "sw": ""
     },
     "col_phase_desc": {
       "en": "Current growth phase (e.g., germination, tillering, stem elongation, grain fill, mature) inferred from crop age",
-      "es-mx": "",
+      "es-mx": "Fase de crecimiento actual (ej., germinación, macollamiento, elongación del tallo, llenado de grano, madurez) inferida de la edad del cultivo",
       "pt-br": "",
       "sw": ""
     },
     "col_germination_progress_desc": {
       "en": "Estimated percentage or stage of germination/emergence based on CI patterns and age. This goes for young fields (age < 4 months). Remains at 100% when finished.",
-      "es-mx": "",
+      "es-mx": "Porcentaje estimado o etapa de germinación/emergencia basado en patrones de IC y edad. Aplica para campos jóvenes (edad < 4 meses). Permanece en 100% al completarse.",
       "pt-br": "",
       "sw": ""
     },
     "col_mean_ci_desc": {
       "en": "Average Chlorophyll Index value across the field; higher values indicate healthier, greener vegetation. Calculated on a 7-day merged weekly image.",
-      "es-mx": "",
+      "es-mx": "Valor promedio del Índice de Clorofila en todo el campo; valores más altos indican vegetación más saludable y verde. Calculado sobre una imagen semanal fusionada de 7 días.",
       "pt-br": "",
       "sw": ""
     },
     "col_weekly_ci_change_desc": {
       "en": "Week-over-week change in Mean_CI; positive values indicate greening/growth, negative values indicate yellowing/decline",
-      "es-mx": "",
+      "es-mx": "Cambio semanal en Mean_CI; valores positivos indican reverdecimiento/crecimiento, valores negativos indican amarillamiento/declive",
       "pt-br": "",
       "sw": ""
     },
     "col_four_week_trend_desc": {
       "en": "Long term change in mean CI; smoothed trend (strong growth, growth, no growth, decline, strong decline)",
-      "es-mx": "",
+      "es-mx": "Cambio a largo plazo en el IC promedio; tendencia suavizada (crecimiento fuerte, crecimiento, sin crecimiento, declive, declive fuerte)",
       "pt-br": "",
       "sw": ""
     },
     "col_ci_range_desc": {
       "en": "Min-max Chlorophyll Index values within the field; wide ranges indicate spatial heterogeneity/patches. Derived from week mosaic.",
-      "es-mx": "",
+      "es-mx": "Valores mínimo-máximo del Índice de Clorofila dentro del campo; rangos amplios indican heterogeneidad espacial/parches. Derivado del mosaico semanal.",
       "pt-br": "",
       "sw": ""
     },
     "col_ci_percentiles_desc": {
       "en": "The CI-range without border effects",
-      "es-mx": "",
+      "es-mx": "El rango de IC sin efectos de borde",
       "pt-br": "",
       "sw": ""
     },
     "col_cv_desc": {
       "en": "Coefficient of variation of CI; measures field uniformity (lower = more uniform, >0.25 = poor uniformity). Derived from week mosaic. In percentages.",
-      "es-mx": "",
+      "es-mx": "Coeficiente de variación del IC; mide la uniformidad del campo (menor = más uniforme, >0.25 = uniformidad deficiente). Derivado del mosaico semanal. En porcentajes.",
       "pt-br": "",
       "sw": ""
     },
     "col_cv_trend_short_desc": {
       "en": "Trend of CV over two weeks. Indicating short-term heterogeneity.",
-      "es-mx": "",
+      "es-mx": "Tendencia del CV en dos semanas. Indica heterogeneidad a corto plazo.",
       "pt-br": "",
       "sw": ""
     },
     "col_cv_trend_long_desc": {
       "en": "Slope of 8-week trend line.",
-      "es-mx": "",
+      "es-mx": "Pendiente de la línea de tendencia de 8 semanas.",
       "pt-br": "",
       "sw": ""
     },
     "col_imminent_prob_desc": {
       "en": "Probability (0-1) that the field is ready for harvest based on LSTM harvest model predictions",
-      "es-mx": "",
+      "es-mx": "Probabilidad (0-1) de que el campo esté listo para cosecha basado en predicciones del modelo de cosecha LSTM",
       "pt-br": "",
       "sw": ""
     },
     "col_cloud_pct_desc": {
       "en": "Percentage of field visible in the satellite image (unobstructed by clouds); lower values indicate poor data quality",
-      "es-mx": "",
+      "es-mx": "Porcentaje del campo visible en la imagen satelital (sin obstrucción de nubes); valores más bajos indican calidad de datos deficiente",
       "pt-br": "",
       "sw": ""
     },
     "col_cloud_category_desc": {
       "en": "Classification of cloud cover level (e.g., clear, partial, heavy); indicates confidence in CI measurements",
-      "es-mx": "",
+      "es-mx": "Clasificación del nivel de cobertura de nubes (ej., despejado, parcial, denso); indica confianza en las mediciones de IC",
       "pt-br": "",
       "sw": ""
     },
     "key_concepts_heading": {
       "en": "# 3. Key Concepts",
-      "es-mx": "",
+      "es-mx": "# 3. Conceptos Clave",
       "pt-br": "",
       "sw": ""
     },
     "growth_phases_heading": {
       "en": "#### *Growth Phases (Age-Based)*",
-      "es-mx": "",
+      "es-mx": "#### *Fases de Crecimiento (Basadas en Edad)*",
       "pt-br": "",
       "sw": ""
     },
     "growth_phases_intro": {
       "en": "Each field is assigned to one of four growth phases based on age in weeks since planting:",
-      "es-mx": "",
+      "es-mx": "Cada campo se asigna a una de cuatro fases de crecimiento según la edad en semanas desde la siembra:",
       "pt-br": "",
       "sw": ""
     },
     "germination_age_range": {
       "en": "0-6 weeks",
-      "es-mx": "",
+      "es-mx": "0-6 semanas",
       "pt-br": "",
       "sw": ""
     },
     "tillering_age_range": {
       "en": "4-16 weeks",
-      "es-mx": "",
+      "es-mx": "4-16 semanas",
       "pt-br": "",
       "sw": ""
     },
     "grand_growth_age_range": {
       "en": "17-39 weeks",
-      "es-mx": "",
+      "es-mx": "17-39 semanas",
       "pt-br": "",
       "sw": ""
     },
     "maturation_age_range": {
       "en": "39+ weeks",
-      "es-mx": "",
+      "es-mx": "39+ semanas",
       "pt-br": "",
       "sw": ""
     },
     "germination_characteristics": {
       "en": "Crop emergence and early establishment; high variability expected",
-      "es-mx": "",
+      "es-mx": "Emergencia del cultivo y establecimiento temprano; alta variabilidad esperada",
       "pt-br": "",
       "sw": ""
     },
     "tillering_characteristics": {
       "en": "Shoot multiplication and plant establishment; rapid growth phase",
-      "es-mx": "",
+      "es-mx": "Multiplicación de brotes y establecimiento de la planta; fase de crecimiento rápido",
       "pt-br": "",
       "sw": ""
     },
     "grand_growth_characteristics": {
       "en": "Peak vegetative growth; maximum height and biomass accumulation",
-      "es-mx": "",
+      "es-mx": "Pico de crecimiento vegetativo; máxima altura y acumulación de biomasa",
       "pt-br": "",
       "sw": ""
     },
     "maturation_characteristics": {
       "en": "Ripening phase; sugar accumulation and preparation for harvest",
-      "es-mx": "",
+      "es-mx": "Fase de maduración; acumulación de azúcar y preparación para la cosecha",
       "pt-br": "",
       "sw": ""
     },
     "status_alert_heading": {
       "en": "#### *Status Alert*",
-      "es-mx": "",
+      "es-mx": "#### *Alerta de Estado*",
       "pt-br": "",
       "sw": ""
     },
     "status_alert_intro": {
       "en": "Status alerts indicate the current field condition based on CI and age-related patterns. Each field receives **one alert** reflecting its most relevant status:",
-      "es-mx": "",
+      "es-mx": "Las alertas de estado indican la condición actual del campo basada en patrones de IC y edad. Cada campo recibe **una alerta** que refleja su estado más relevante:",
       "pt-br": "",
       "sw": ""
     },
     "harvest_ready_condition": {
       "en": "Harvest model > 0.50 and crop is mature",
-      "es-mx": "",
+      "es-mx": "Modelo de cosecha > 0.50 y el cultivo está maduro",
       "pt-br": "",
       "sw": ""
     },
     "harvest_ready_phase_info": {
       "en": "Active from 52 weeks onwards",
-      "es-mx": "",
+      "es-mx": "Activo a partir de las 52 semanas",
       "pt-br": "",
       "sw": ""
     },
     "harvest_ready_msg": {
       "en": "Ready for harvest-check",
-      "es-mx": "",
+      "es-mx": "Listo para verificación de cosecha",
       "pt-br": "",
       "sw": ""
     },
     "harvested_bare_condition": {
       "en": "Field of 50 weeks or older either shows mean CI values lower than 1.5 (for a maximum of three weeks) OR drops from higher CI to lower than 1.5. Alert drops if CI rises and passes 1.5 again",
-      "es-mx": "",
+      "es-mx": "Campo de 50 semanas o más que muestra valores promedio de IC menores a 1.5 (por un máximo de tres semanas) O que cae de un IC mayor a menor de 1.5. La alerta se desactiva si el IC sube y supera 1.5 nuevamente",
       "pt-br": "",
       "sw": ""
     },
     "harvested_bare_phase_info": {
       "en": "Maturation (39+)",
-      "es-mx": "",
+      "es-mx": "Maduración (39+)",
       "pt-br": "",
       "sw": ""
     },
     "harvested_bare_msg": {
       "en": "Harvested or bare field",
-      "es-mx": "",
+      "es-mx": "Campo cosechado o descubierto",
       "pt-br": "",
       "sw": ""
     },
     "stress_condition": {
       "en": "Mean CI on field drops by 2+ points but field mean CI remains higher than 1.5",
-      "es-mx": "",
+      "es-mx": "El IC promedio del campo cae 2+ puntos pero el IC promedio del campo permanece mayor a 1.5",
       "pt-br": "",
       "sw": ""
     },
     "stress_phase_info": {
       "en": "Any",
-      "es-mx": "",
+      "es-mx": "Cualquiera",
       "pt-br": "",
       "sw": ""
     },
     "stress_msg": {
       "en": "Strong decline in crop health",
-      "es-mx": "",
+      "es-mx": "Declive fuerte en la salud del cultivo",
       "pt-br": "",
       "sw": ""
     },
     "harvest_date_heading": {
       "en": "#### *Harvest Date and Harvest Imminent*",
-      "es-mx": "",
+      "es-mx": "#### *Fecha de Cosecha y Cosecha Inminente*",
       "pt-br": "",
       "sw": ""
     },
     "harvest_date_body_1": {
       "en": "The SmartCane algorithm calculates the last harvest date and the probability of harvest approaching in the next 4 weeks. Two different algorithms are used.",
-      "es-mx": "",
+      "es-mx": "El algoritmo de SmartCane calcula la última fecha de cosecha y la probabilidad de que la cosecha se aproxime en las próximas 4 semanas. Se utilizan dos algoritmos diferentes.",
       "pt-br": "",
       "sw": ""
     },
     "harvest_date_body_2": {
       "en": "The **last harvest date** is a timeseries analysis of the CI levels of the past years, based on clean factory managed fields as data set for the machine learning, a reliability of over 90% has been reached. Smallholder managed fields of small size (0.3 acres) have specific side effects and field management characteristics, that influence the model results.",
-      "es-mx": "",
+      "es-mx": "La **última fecha de cosecha** es un análisis de series de tiempo de los niveles de IC de años anteriores, basado en campos limpios administrados por la fábrica como conjunto de datos para el aprendizaje automático, se ha alcanzado una confiabilidad superior al 90%. Los campos pequeños administrados por pequeños productores (0.3 acres) tienen efectos secundarios específicos y características de manejo que influyen en los resultados del modelo.",
       "pt-br": "",
       "sw": ""
     },
     "harvest_date_body_3": {
       "en": "**Imminent_probability** of harvest is a prediction algorithm, estimating the likelihood of a crop ready to be harvested in the near future. This prediction takes the CI-levels into consideration, building on the vegetative development of sugarcane in the last stage of Maturation, where all sucrose is pulled into the stalk, depleting the leaves from energy and productive function, reducing the levels of CI in the leaf tissue.",
-      "es-mx": "",
+      "es-mx": "**Imminent_probability** de cosecha es un algoritmo de predicción que estima la probabilidad de que un cultivo esté listo para ser cosechado en el futuro cercano. Esta predicción toma en cuenta los niveles de IC, basándose en el desarrollo vegetativo de la caña de azúcar en la última etapa de Maduración, donde toda la sacarosa se concentra en el tallo, agotando la energía y función productiva de las hojas, reduciendo los niveles de IC en el tejido foliar.",
       "pt-br": "",
       "sw": ""
     },
     "harvest_date_body_4": {
       "en": "Both algorithms are not always in sync, and can have contradictory results. Wider field characteristics analysis is suggested if such contradictory calculation results occur.",
-      "es-mx": "",
+      "es-mx": "Ambos algoritmos no siempre están sincronizados y pueden tener resultados contradictorios. Se sugiere un análisis más amplio de las características del campo si se presentan resultados de cálculo contradictorios.",
       "pt-br": "",
       "sw": ""
     }
@@ -544,127 +544,127 @@
   "status_translations": {
     "PHASE DISTRIBUTION": {
       "en": "Phase Distribution",
-      "es-mx": "",
+      "es-mx": "Distribución de Fases",
       "pt-br": "",
       "sw": ""
     },
     "OPERATIONAL ALERTS": {
       "en": "Operational Alerts",
-      "es-mx": "",
+      "es-mx": "Alertas Operativas",
       "pt-br": "",
       "sw": ""
     },
     "AREA CHANGE": {
       "en": "Area Change",
-      "es-mx": "",
+      "es-mx": "Cambio de Área",
       "pt-br": "",
       "sw": ""
     },
     "CLOUD INFLUENCE": {
       "en": "Cloud Influence",
-      "es-mx": "",
+      "es-mx": "Influencia de Nubes",
       "pt-br": "",
       "sw": ""
     },
     "TOTAL FARM": {
       "en": "Total Farm",
-      "es-mx": "",
+      "es-mx": "Total de la Finca",
       "pt-br": "",
       "sw": ""
     },
     "Germination": {
       "en": "Germination",
-      "es-mx": "",
+      "es-mx": "Germinación",
       "pt-br": "",
       "sw": ""
     },
     "Tillering": {
       "en": "Tillering",
-      "es-mx": "",
+      "es-mx": "Macollamiento",
       "pt-br": "",
       "sw": ""
     },
     "Grand Growth": {
       "en": "Grand Growth",
-      "es-mx": "",
+      "es-mx": "Gran Crecimiento",
       "pt-br": "",
       "sw": ""
     },
     "Maturation": {
       "en": "Maturation",
-      "es-mx": "",
+      "es-mx": "Maduración",
       "pt-br": "",
       "sw": ""
     },
     "Unknown Phase": {
       "en": "Unknown Phase",
-      "es-mx": "",
+      "es-mx": "Fase Desconocida",
       "pt-br": "",
       "sw": ""
     },
     "harvest_ready": {
       "en": "Ready for Harvest-Check",
-      "es-mx": "",
+      "es-mx": "Listo para Verificación de Cosecha",
       "pt-br": "",
       "sw": ""
     },
     "harvested_bare": {
       "en": "Harvested / Bare Field",
-      "es-mx": "",
+      "es-mx": "Cosechado / Campo Descubierto",
       "pt-br": "",
       "sw": ""
     },
     "stress_detected": {
       "en": "Stress Detected",
-      "es-mx": "",
+      "es-mx": "Estrés Detectado",
       "pt-br": "",
       "sw": ""
     },
     "germination_delayed": {
       "en": "Germination Delayed",
-      "es-mx": "",
+      "es-mx": "Germinación Retrasada",
       "pt-br": "",
       "sw": ""
     },
     "growth_on_track": {
       "en": "Growth on Track",
-      "es-mx": "",
+      "es-mx": "Crecimiento en Curso",
       "pt-br": "",
       "sw": ""
     },
     "No active triggers": {
       "en": "No Active Triggers",
-      "es-mx": "",
+      "es-mx": "Sin Alertas Activas",
       "pt-br": "",
       "sw": ""
     },
     "Improving": {
       "en": "Improving",
-      "es-mx": "",
+      "es-mx": "Mejorando",
       "pt-br": "",
       "sw": ""
     },
     "Stable": {
       "en": "Stable",
-      "es-mx": "",
+      "es-mx": "Estable",
       "pt-br": "",
       "sw": ""
     },
     "Declining": {
       "en": "Declining",
-      "es-mx": "",
+      "es-mx": "En Declive",
       "pt-br": "",
       "sw": ""
     },
     "Total Acreage": {
       "en": "Total Acreage",
-      "es-mx": "",
+      "es-mx": "Superficie Total",
       "pt-br": "",
       "sw": ""
     },
     "kpi_na": {
       "en": "KPI data not available for {project_dir} on this date.",
-      "es-mx": "",
+      "es-mx": "Datos de KPI no disponibles para {project_dir} en esta fecha.",
       "pt-br": "",
       "sw": ""
     }
@@ -672,159 +672,159 @@
   "figure_translations": {
     "kpi_col_category": {
       "en": "KPI Category",
-      "es-mx": "",
+      "es-mx": "Categoría de KPI",
       "pt-br": "",
       "sw": ""
     },
     "kpi_col_item": {
       "en": "Item",
-      "es-mx": "",
+      "es-mx": "Elemento",
       "pt-br": "",
       "sw": ""
     },
     "kpi_col_acreage": {
       "en": "Acreage",
-      "es-mx": "",
+      "es-mx": "Superficie",
       "pt-br": "",
       "sw": ""
     },
     "kpi_col_percent": {
       "en": "Percentage of Total Fields",
-      "es-mx": "",
+      "es-mx": "Porcentaje del Total de Campos",
       "pt-br": "",
       "sw": ""
     },
     "kpi_col_fields": {
       "en": "# Fields",
-      "es-mx": "",
+      "es-mx": "# Campos",
       "pt-br": "",
       "sw": ""
     },
     "hexbin_legend_acres": {
       "en": "Total Acres",
-      "es-mx": "",
+      "es-mx": "Acres Totales",
       "pt-br": "",
       "sw": ""
     },
     "hexbin_subtitle": {
       "en": "Acres of fields 'harvest ready within a month'",
-      "es-mx": "",
+      "es-mx": "Acres de campos listos para cosecha dentro de un mes",
       "pt-br": "",
       "sw": ""
     },
     "hexbin_not_ready": {
       "en": "Not harvest ready",
-      "es-mx": "",
+      "es-mx": "No listo para cosecha",
       "pt-br": "",
       "sw": ""
     },
     "metadata_caption": {
       "en": "Report Metadata",
-      "es-mx": "",
+      "es-mx": "Metadatos del Reporte",
       "pt-br": "",
       "sw": ""
     },
     "metric": {
       "en": "Metric",
-      "es-mx": "",
+      "es-mx": "Métrica",
       "pt-br": "",
       "sw": ""
     },
     "value": {
       "en": "Value",
-      "es-mx": "",
+      "es-mx": "Valor",
       "pt-br": "",
       "sw": ""
     },
     "report_gen": {
       "en": "Report Generated",
-      "es-mx": "",
+      "es-mx": "Reporte Generado",
       "pt-br": "",
       "sw": ""
     },
     "data_source": {
       "en": "Data Source",
-      "es-mx": "",
+      "es-mx": "Fuente de Datos",
       "pt-br": "",
       "sw": ""
     },
     "analysis_period": {
       "en": "Analysis Period",
-      "es-mx": "",
+      "es-mx": "Período de Análisis",
       "pt-br": "",
       "sw": ""
     },
     "tot_fields_number": {
       "en": "Total Fields [number]",
-      "es-mx": "",
+      "es-mx": "Total de Campos [número]",
       "pt-br": "",
       "sw": ""
     },
     "tot_area_acres": {
       "en": "Total Area [acres]",
-      "es-mx": "",
+      "es-mx": "Área Total [acres]",
       "pt-br": "",
       "sw": ""
     },
     "next_update": {
       "en": "Next Update",
-      "es-mx": "",
+      "es-mx": "Próxima Actualización",
       "pt-br": "",
       "sw": ""
     },
     "service_provided": {
       "en": "Service provided",
-      "es-mx": "",
+      "es-mx": "Servicio proporcionado",
       "pt-br": "",
       "sw": ""
     },
     "service_start": {
       "en": "Starting date service",
-      "es-mx": "",
+      "es-mx": "Fecha de inicio del servicio",
       "pt-br": "",
       "sw": ""
     },
     "next_wed": {
       "en": "Next Wednesday",
-      "es-mx": "",
+      "es-mx": "Próximo miércoles",
       "pt-br": "",
       "sw": ""
     },
     "week": {
       "en": "Week",
-      "es-mx": "",
+      "es-mx": "Semana",
       "pt-br": "",
       "sw": ""
     },
     "of": {
       "en": "of",
-      "es-mx": "",
+      "es-mx": "de",
       "pt-br": "",
       "sw": ""
     },
     "project": {
       "en": "Project",
-      "es-mx": "",
+      "es-mx": "Proyecto",
       "pt-br": "",
       "sw": ""
     },
     "cane_supply_service": {
       "en": "Cane Supply Office - Weekly",
-      "es-mx": "",
+      "es-mx": "Oficina de Abastecimiento de Caña - Semanal",
       "pt-br": "",
       "sw": ""
     },
     "unknown": {
       "en": "Unknown",
-      "es-mx": "",
+      "es-mx": "Desconocido",
       "pt-br": "",
       "sw": ""
     },
     "no_field_data": {
       "en": "No field-level KPI data available for this report period.",
-      "es-mx": "",
+      "es-mx": "No hay datos de KPI a nivel de campo disponibles para este período de reporte.",
       "pt-br": "",
       "sw": ""
     }
   }
-}
+}