ClearSKY provides near real‑time, daily–weekly cloud‑free mosaics compatible with Sentinel‑2. Your Nimbus composites are distributed as:

TOA (top‑of‑atmosphere) internally, and
BOA (bottom‑of‑atmosphere) for download in most sample packs

This recipe uses BOA GeoTIFFs (analysis‑ready). Sentinel‑2 BOA typically contains 12 bands — all bands except B10 (cirrus). That means:

B04 = Red, B08 = NIR (used for NDVI)
B03 = Green (used for NDWI)
B11 / B12 = SWIR (used for MNDWI / NBR etc.)

🧠 Tip: reflectance scaling can vary by source (e.g., 0–1 float vs 0–10000 int). The code below will auto‑scale to 0–1 when it detects values > 1.

What you'll need

macOS/Linux shell (or Windows PowerShell/WSL)
Python 3.9+ with: rasterio, numpy, pandas, matplotlib, shapely (optional), fiona (optional)
Access to ClearSKY sample datasets (non‑commercial)

1) Download a time series (daily–weekly BOA)

Replace the URL with a signed link from /samples after email unlock.

# Example: parcel time‑series pack (contains date‑stamped BOA mosaics)
curl -L -o parcel_timeseries.zip "REPLACE_WITH_SIGNED_URL"
unzip -q parcel_timeseries.zip -d timeseries

# Files should look like:
# timeseries/
#   20240410.tif
#   20240420.tif
#   20240430.tif
#   ...

2) (Optional) Provide a field polygon

If you want per‑field statistics, save a GeoJSON file named aoi.geojson with your polygon:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": { "name": "My Field" },
      "geometry": {
        "type": "Polygon",
        "coordinates": [[[ -50.11, -9.47 ], [ -50.10, -9.47 ], [ -50.10, -9.48 ], [ -50.11, -9.48 ], [ -50.11, -9.47 ]]]
      }
    }
  ]
}

CRS: EPSG:4326 (lon/lat) is fine — the script will reproject the AOI on the fly.

3) Run the NDVI time‑series extractor

Download the script below and run it. It computes NDVI, EVI and NDWI per date (masked to your AOI if provided), then writes a CSV and a PNG chart.

python calc_ndvi_timeseries.py --dir timeseries --aoi aoi.geojson
# or without AOI
python calc_ndvi_timeseries.py --dir timeseries

Outputs:

indices_timeseries.csv — date, ndvi_mean, evi_mean, ndwi_mean
ndvi_trend.png — simple NDVI trend chart

4) Band reference (BOA composites)

These composites have 12 bands (B10 omitted). Use 1‑based band indices in Rasterio:

Index	Sentinel‑2 band	Common name	Notes
1	B01	Coastal aerosol	60 m
2	B02	Blue	10 m
3	B03	Green	10 m
4	B04	Red	10 m
5	B05	Red‑edge 1	20 m
6	B06	Red‑edge 2	20 m
7	B07	Red‑edge 3	20 m
8	B08	NIR	10 m
9	B8A	NIR narrow	20 m
10	B09	Water vapor	60 m
11	B11	SWIR 1	20 m
12	B12	SWIR 2	20 m

B10 (cirrus) is not present in BOA delivery.

5) When you go live

Swap sample URLs for ordered AOIs.
Tune index thresholds per crop/region; aggregate by 7–10 day windows for smoother trends.
Use ClearSKY webhooks (or a cron) to re‑run when a new daily–weekly mosaic lands.

Full script (Python)

Save as calc_ndvi_timeseries.py:

#!/usr/bin/env python3
import argparse, glob, json, os, re
import numpy as np
import pandas as pd
import rasterio as rio
from rasterio.mask import mask as rio_mask
import matplotlib.pyplot as plt

def autoscale_reflectance(arr):
    # If the data looks like 0..10000 integers, scale to 0..1
    a = arr.astype("float32")
    finite = np.isfinite(a)
    if finite.any() and np.nanmax(a[finite]) > 2.0:
        a = a / 10000.0
    return a

def read_band(ds, idx):
    b = ds.read(idx).astype("float32")
    return autoscale_reflectance(b)

def compute_indices(bands):
    # bands dict must contain b2, b3, b4, b8 at minimum
    red = bands["b4"]; nir = bands["b8"]; blue = bands["b2"]; green = bands["b3"]

    # NDVI
    ndvi = (nir - red) / (nir + red + 1e-6)

    # EVI (L=1, C1=6, C2=7.5, G=2.5) using BOA
    evi = 2.5 * (nir - red) / (nir + 6.0*red - 7.5*blue + 1.0 + 1e-6)

    # NDWI (McFeeters; water index): (G - NIR)/(G + NIR)
    ndwi = (green - nir) / (green + nir + 1e-6)

    return ndvi, evi, ndwi

def read_aoi_geojson(path):
    if not path: return None
    with open(path, "r", encoding="utf-8") as f:
        gj = json.load(f)
    # Accept FeatureCollection or Feature
    if gj.get("type") == "FeatureCollection":
        geoms = [f["geometry"] for f in gj.get("features", []) if f.get("geometry")]
    elif gj.get("type") == "Feature":
        geoms = [gj["geometry"]]
    else:
        geoms = [gj]
    return geoms

def list_tifs(dirpath):
    tif_paths = sorted(glob.glob(os.path.join(dirpath, "*.tif")) + glob.glob(os.path.join(dirpath, "*.tiff")))
    # Expect filenames like YYYY-MM-DD_*.tif
    def date_from_name(p):
        m = re.search(r"(20\d{2}-\d{2}-\d{2})", os.path.basename(p))
        return m.group(1) if m else None
    items = [(date_from_name(p), p) for p in tif_paths]
    items = [(d, p) for (d, p) in items if d is not None]
    items.sort(key=lambda x: x[0])
    return items

def masked_mean(arr, aoi_mask=None):
    a = arr
    if aoi_mask is not None:
        a = np.where(aoi_mask, a, np.nan)
    return float(np.nanmean(a))

def main():
    ap = argparse.ArgumentParser(description="Compute NDVI/EVI/NDWI time series from BOA composites")
    ap.add_argument("--dir", required=True, help="Folder with *_BOA.tif files (one per date)")
    ap.add_argument("--aoi", help="Optional AOI GeoJSON (lon/lat)")
    ap.add_argument("--out_csv", default="indices_timeseries.csv")
    ap.add_argument("--out_png", default="ndvi_trend.png")
    args = ap.parse_args()

    items = list_tifs(args.dir)
    if not items:
        raise SystemExit("No GeoTIFFs found in --dir")

    geoms = read_aoi_geojson(args.aoi)

    rows = []
    for date, path in items:
        with rio.open(path) as ds:
            # If AOI present, reproject AOI and build mask using band 4 shape
            aoi_mask = None
            if geoms:
                try:
                    data_mask, _ = rio_mask(ds, geoms, crop=False, filled=False)
                    aoi_mask = data_mask[0].astype(bool)
                except Exception:
                    aoi_mask = None  # proceed without

            # BOA band indices (1-based): B2=2, B3=3, B4=4, B8=8
            bands = {
                "b2": read_band(ds, 2),
                "b3": read_band(ds, 3),
                "b4": read_band(ds, 4),
                "b8": read_band(ds, 8),
            }
            ndvi, evi, ndwi = compute_indices(bands)

            rows.append({
                "date": date,
                "ndvi_mean": masked_mean(ndvi, aoi_mask),
                "evi_mean": masked_mean(evi, aoi_mask),
                "ndwi_mean": masked_mean(ndwi, aoi_mask),
            })

    df = pd.DataFrame(rows).sort_values("date")
    df.to_csv(args.out_csv, index=False)

    # Plot NDVI trend
    plt.figure(figsize=(9, 4.5))
    plt.plot(pd.to_datetime(df["date"]), df["ndvi_mean"], marker="o")
    plt.title("NDVI Trend")
    plt.xlabel("Date")
    plt.ylabel("NDVI (mean)")
    plt.grid(True, linestyle="--", alpha=0.4)
    plt.tight_layout()
    plt.savefig(args.out_png, dpi=120)
    print(f"Wrote {args.out_csv} and {args.out_png}")

if __name__ == "__main__":
    main()

Notes

Cadence: These mosaics are generated daily–weekly, so your chart will update often.
Near real‑time: New BOA composites appear shortly after overpass; schedule the script accordingly.
Bands: BOA has 12 bands (B10 missing). We used B4/B8 for NDVI, B2/B4/B8 for EVI, B3/B8 for NDWI.

Happy building!