What you’ll build

You’ll create a Python workflow to analyze crop trends over time using ClearSKY’s cloud-free Sentinel-2 mosaics. This includes loading weekly (or daily) mosaics, masking out NoData, extracting NDVI for a defined field or plot, and visualizing crop health changes through time. The workflow is robust for operational agriculture monitoring, leveraging cloud- and shadow-free, BoA-equivalent data for clean temporal analytics.

The timeseries approach is invaluable for:

• Early detection of crop stress or abnormal change
• Monitoring seasonal progression
• Building up data for later yield modeling or anomaly detection

Requirements

• Python 3.8+
• rasterio
• xarray and rioxarray
• numpy
• matplotlib

Install all requirements:

pip install rasterio xarray rioxarray numpy matplotlib

• Obtain ClearSKY sample data by requesting the Agriculture sample pack at /developers.
  You may substitute clean, cloud-free original Sentinel-2 L2A imagery, ensuring correct band order (B10/cirrus absent).

Load sample data

You’ll need a stack of ClearSKY mosaics, covering the same AOI (tile) over several dates. Each file is a Sentinel-2 L2A-like GeoTIFF with:

• Resolution: 10 m across all bands
• No clouds and no cloud shadows – quality-masked in fusion
• Band order: Sentinel-2 L2A (B10/cirrus absent)
• NoData: only outside AOI, value -32768
• Bands: int16 reflectance, scale ÷10000 to [0,1]

Place each GeoTIFF in a folder, naming them e.g.,

• clearsky_AOI_20240401.tif
• clearsky_AOI_20240408.tif
• clearsky_AOI_20240415.tif
• ...

Prepare a field boundary (optional)

Clip to a single field by providing a shapefile, GeoJSON, or even simple pixel coordinates.

Core logic

This recipe loads the ClearSKY stacks, extracts NDVI for each date, masks NoData, and plots NDVI evolution for the field.

Helper functions

Use these helpers for proper scaling and masking:

import numpy as np
NODATA = -32768

def read_band_reflectance(ds, idx):
    a = ds.read(idx).astype("float32")
    nodata = ds.nodata if ds.nodata is not None else NODATA
    a = np.where(a == nodata, np.nan, a)
    # Sentinel-2 reflectance is scaled by 10000; convert to [0,1]
    if np.nanmax(a) > 1.0:
        a = a / 10000.0
    return a

def image_date_from_tags(ds):
    v = ds.tags().get("IMAGE_DATE")
    if v and len(v) == 8:
        return f"{v[0:4]}-{v[4:6]}-{v[6:8]}"
    return None

Reading and analyzing the time stack

Assuming all GeoTIFFs are in ./clearsky_timeseries/:

import glob
import rasterio
import numpy as np
import matplotlib.pyplot as plt

folder = './clearsky_timeseries/'
filelist = sorted(glob.glob(folder + '*.tif'))
print(f"Found {len(filelist)} mosaics.")

ndvi_list = []
dates = []

for fname in filelist:
    with rasterio.open(fname) as ds:
        print(f"Processing {fname} ...")
        # Sentinel-2 L2A order; NIR=B08 (idx=8), RED=B04 (idx=4)
        nir = read_band_reflectance(ds, 8)
        red = read_band_reflectance(ds, 4)
        # NDVI = (NIR - RED) / (NIR + RED), avoid zero divide
        ndvi = np.where((nir + red) > 0, (nir - red) / (nir + red), np.nan)
        # Mask field region if desired (e.g., use a binary mask array)
        ndvi_list.append(ndvi)
        # Extract date from ClearSKY tag
        date_str = image_date_from_tags(ds)
        dates.append(date_str)

(Optional) Add a spatial mask for your field

This can be a binary numpy array with shape matching your mosaic, where True marks the area-of-interest; e.g.,

# For demonstration: select a center 50x50 window
field_mask = np.zeros_like(ndvi, dtype=bool)
ysize, xsize = ndvi.shape
field_mask[ysize//2-25:ysize//2+25, xsize//2-25:xsize//2+25] = True

To apply to each NDVI:

# Apply mask after NDVI computation
mean_ndvi_timeseries = [np.nanmean(nd[ field_mask ]) for nd in ndvi_list]

If working with entire AOI, just take np.nanmean(ndvi) across each image.

Visualize/Export

Plot NDVI trend

plt.figure(figsize=(10,5))
plt.plot(dates, mean_ndvi_timeseries, marker='o', linewidth=2)
plt.xlabel('Date')
plt.ylabel('Mean NDVI')
plt.title('Crop NDVI Timeseries (ClearSKY, cloud-free)')
plt.xticks(rotation=45)
plt.grid(True)
plt.tight_layout()
plt.show()

Show NDVI spatial maps for select dates

import matplotlib.cm as cm
for i, ndvi in enumerate(ndvi_list[::len(ndvi_list)//3]):  # Pick 3 dates
    plt.figure()
    plt.imshow(ndvi, cmap=cm.YlGn, vmin=0, vmax=1)
    plt.title(f'NDVI map ({dates[::len(ndvi_list)//3][i]})')
    plt.colorbar(label='NDVI')
    plt.tight_layout()
plt.show()

Troubleshooting

• If you see all-NaN NDVI values, ensure you’re masking only for NoData (-32768)—not valid zero reflectance.
• When band values seem too large or small, recheck that you’ve properly scaled: int16 reflectance ÷10000 → [0,1].
• NDVI outside expected [-1,+1]? Clip manually if necessary; if you use ClearSKY’s precomputed NDVI from API, scale by ÷32767 → [-1,1].
• If arrays don't match, check that all images use the same AOI/extent and projection—ClearSKY sample packs are pre-mosaicked/gridded.
• Daily–weekly cadence means there may be small date gaps.

Next steps

• Expand to multiple fields and plot multiple timeseries per field or crop type.
• Compute rolling averages or detect anomalies (e.g., with scipy.signal.find_peaks).
• Export trend results as CSV for downstream use.
• Integrate with other indices: NDWI, EVI, etc. (see docs.clearsky.vision for details).
• For API access and automation, see api.clearsky.vision.

Ready to scale? To ingest Agriculture sample packs or subscribe to production mosaics, see /developers.