In [1]:
import planetary_computer
import pystac_client
import stackstac
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings('ignore')
print("All imports OK")
All imports OK
In [2]:
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
# AOI: Jalisco agricultural valleys (Tequila / Amatitán region)
bbox = [-104.0, 20.5, -103.5, 20.9]
search = catalog.search(
collections=["sentinel-2-l2a"],
bbox=bbox,
datetime="2024-03-01/2024-05-31",
query={"eo:cloud_cover": {"lt": 10}},
)
items = list(search.items())
print(f"Scenes found: {len(items)}")
for item in items[:5]:
print(f" {item.id} — cloud cover: {item.properties['eo:cloud_cover']}%")
Scenes found: 50 S2B_MSIL2A_20240531T171859_R012_T13QFD_20240531T223148 — cloud cover: 2.33111% S2B_MSIL2A_20240531T171859_R012_T13QFC_20240531T222356 — cloud cover: 1.222408% S2B_MSIL2A_20240531T171859_R012_T13QEC_20240531T222341 — cloud cover: 0.830179% S2A_MSIL2A_20240526T171901_R012_T13QFD_20240527T012538 — cloud cover: 0.002644% S2A_MSIL2A_20240526T171901_R012_T13QFC_20240527T011215 — cloud cover: 3.91525%
In [4]:
# Best scene: 0.002% cloud cover
best_item = items[3]
print(f"Using: {best_item.id}")
print(f"Cloud cover: {best_item.properties['eo:cloud_cover']}%")
print(f"Date: {best_item.properties['datetime']}")
# Stack the bands we need:
# B04 = Red, B08 = NIR, B03 = Green, B05 = Red Edge, B11 = SWIR
stack = stackstac.stack(
[best_item],
assets=["B04", "B08", "B03", "B05", "B11"],
epsg=32614,
resolution=10,
bounds_latlon=bbox,
)
print(f"\nStack shape: {stack.shape}")
print(f"Bands: {stack.coords['band'].values}")
Using: S2A_MSIL2A_20240526T171901_R012_T13QFD_20240527T012538 Cloud cover: 0.002644% Date: 2024-05-26T17:19:01.024000Z Stack shape: (1, 5, 4592, 5351) Bands: ['B04' 'B08' 'B03' 'B05' 'B11']
In [5]:
import xarray as xr
# Load into memory (this may take 30-60 seconds)
print("Loading bands into memory...")
data = stack.compute()
# Extract individual bands
B04 = data.sel(band="B04").values[0].astype(float) # Red
B08 = data.sel(band="B08").values[0].astype(float) # NIR
B03 = data.sel(band="B03").values[0].astype(float) # Green
B05 = data.sel(band="B05").values[0].astype(float) # Red Edge
B11 = data.sel(band="B11").values[0].astype(float) # SWIR
print(f"Bands loaded. Shape: {B04.shape}")
print(f"B04 (Red) min/max: {B04.min():.0f} / {B04.max():.0f}")
print(f"B08 (NIR) min/max: {B08.min():.0f} / {B08.max():.0f}")
Loading bands into memory... Bands loaded. Shape: (4592, 5351) B04 (Red) min/max: nan / nan B08 (NIR) min/max: nan / nan
In [6]:
print(f"B04 dtype: {B04.dtype}")
print(f"Total pixels: {B04.size}")
print(f"NaN count: {np.isnan(B04).sum()}")
print(f"Zero count: {(B04 == 0).sum()}")
# Check raw stack fill value
print(f"\nStack fill value: {stack.attrs.get('_FillValue', 'not set')}")
print(f"Stack nodata: {stack.attrs.get('nodata', 'not set')}")
# Check a raw sample before compute
raw_sample = stack.isel(band=0, x=slice(2000,2010), y=slice(2000,2010)).values
print(f"\nRaw sample values:\n{raw_sample}")
B04 dtype: float64 Total pixels: 24571792 NaN count: 12675138 Zero count: 0 Stack fill value: not set Stack nodata: not set Raw sample values: [[[2732. 2758. 2656. 2754. 2770. 2556. 2422. 2372. 2448. 2570.] [2750. 2770. 2750. 2804. 2832. 2574. 2466. 2374. 2406. 2544.] [2786. 2810. 2790. 2832. 2782. 2508. 2458. 2378. 2412. 2516.] [2804. 2828. 2788. 2816. 2698. 2478. 2438. 2380. 2418. 2526.] [2824. 2810. 2824. 2852. 2660. 2462. 2416. 2422. 2494. 2556.] [2914. 2896. 2844. 2810. 2566. 2450. 2402. 2532. 2642. 2676.] [2948. 2874. 2826. 2686. 2472. 2468. 2410. 2522. 2580. 2626.] [2984. 2860. 2860. 2808. 2606. 2478. 2398. 2434. 2460. 2494.] [2956. 2924. 2992. 3018. 2780. 2512. 2432. 2436. 2466. 2472.] [2884. 2942. 3032. 3038. 2858. 2612. 2532. 2494. 2526. 2494.]]]
In [7]:
# Mask NaN pixels consistently
valid = (~np.isnan(B04)) & (~np.isnan(B08)) & (~np.isnan(B03)) & (~np.isnan(B05)) & (~np.isnan(B11))
# Initialize output arrays with NaN
NDVI = np.full(B04.shape, np.nan)
NDWI = np.full(B04.shape, np.nan)
NDRE = np.full(B04.shape, np.nan)
SAVI = np.full(B04.shape, np.nan)
L = 0.5 # SAVI soil correction factor
# Calculate only on valid pixels
NDVI[valid] = (B08[valid] - B04[valid]) / (B08[valid] + B04[valid])
NDWI[valid] = (B03[valid] - B08[valid]) / (B03[valid] + B08[valid])
NDRE[valid] = (B08[valid] - B05[valid]) / (B08[valid] + B05[valid])
SAVI[valid] = ((B08[valid] - B04[valid]) / (B08[valid] + B04[valid] + L)) * (1 + L)
print("Indices calculated.")
print(f"NDVI — mean: {np.nanmean(NDVI):.3f} | min: {np.nanmin(NDVI):.3f} | max: {np.nanmax(NDVI):.3f}")
print(f"NDWI — mean: {np.nanmean(NDWI):.3f} | min: {np.nanmin(NDWI):.3f} | max: {np.nanmax(NDWI):.3f}")
print(f"NDRE — mean: {np.nanmean(NDRE):.3f} | min: {np.nanmin(NDRE):.3f} | max: {np.nanmax(NDRE):.3f}")
print(f"SAVI — mean: {np.nanmean(SAVI):.3f} | min: {np.nanmin(SAVI):.3f} | max: {np.nanmax(SAVI):.3f}")
Indices calculated. NDVI — mean: 0.186 | min: -0.585 | max: 0.837 NDWI — mean: -0.254 | min: -0.658 | max: 0.545 NDRE — mean: 0.117 | min: -0.756 | max: 0.761 SAVI — mean: 0.279 | min: -0.877 | max: 1.255
In [10]:
# Better RGB: NIR=R, Red=G, Green=B → classic false color for agriculture
rgb_r = B08[::step, ::step] # NIR → Red channel
rgb_g = B04[::step, ::step] # Red → Green channel
rgb_b = B03[::step, ::step] # Green → Blue channel
def normalize(arr):
a = np.nanpercentile(arr, 2)
b = np.nanpercentile(arr, 98)
return np.clip((arr - a) / (b - a), 0, 1)
rgb = np.dstack([normalize(rgb_r), normalize(rgb_g), normalize(rgb_b)])
# Black out nodata
nan_mask = np.isnan(NDVI[::step, ::step])
rgb[nan_mask] = 0
plt.style.use('dark_background')
fig = plt.figure(figsize=(20, 14))
fig.patch.set_facecolor('#0d0d0d')
gs = GridSpec(2, 3, figure=fig, hspace=0.35, wspace=0.25)
ax0 = fig.add_subplot(gs[0, 0])
ax1 = fig.add_subplot(gs[0, 1])
ax2 = fig.add_subplot(gs[0, 2])
ax3 = fig.add_subplot(gs[1, 0])
ax4 = fig.add_subplot(gs[1, 1])
for ax in [ax0, ax1, ax2, ax3, ax4]:
ax.set_facecolor('#0d0d0d')
# False color composite (NIR/Red/Green)
ax0.imshow(rgb, interpolation='nearest')
ax0.set_title('False Color Composite\n(R=NIR G=Red B=Green)', color='white', fontsize=11, pad=10)
ax0.axis('off')
# NDVI — tighter range to show contrast
im1 = ax1.imshow(ndvi_plot, cmap='RdYlGn', vmin=-0.1, vmax=0.6, interpolation='nearest')
ax1.set_title('NDVI — Vegetation Health\n(NIR−Red / NIR+Red)', color='white', fontsize=11, pad=10)
ax1.axis('off')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04).ax.yaxis.set_tick_params(color='white', labelcolor='white')
# NDWI
im2 = ax2.imshow(ndwi_plot, cmap='RdBu', vmin=-0.5, vmax=0.3, interpolation='nearest')
ax2.set_title('NDWI — Water Stress\n(Green−NIR / Green+NIR)', color='white', fontsize=11, pad=10)
ax2.axis('off')
plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04).ax.yaxis.set_tick_params(color='white', labelcolor='white')
# NDRE
im3 = ax3.imshow(ndre_plot, cmap='PiYG', vmin=-0.3, vmax=0.6, interpolation='nearest')
ax3.set_title('NDRE — Chlorophyll (Red Edge)\n(NIR−RedEdge / NIR+RedEdge)', color='white', fontsize=11, pad=10)
ax3.axis('off')
plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04).ax.yaxis.set_tick_params(color='white', labelcolor='white')
# SAVI
im4 = ax4.imshow(savi_plot, cmap='YlGn', vmin=-0.2, vmax=0.8, interpolation='nearest')
ax4.set_title('SAVI — Soil-Adjusted Vegetation\n(L=0.5)', color='white', fontsize=11, pad=10)
ax4.axis('off')
plt.colorbar(im4, ax=ax4, fraction=0.046, pad=0.04).ax.yaxis.set_tick_params(color='white', labelcolor='white')
fig.suptitle(
'Sprint 04 — Multispectral Agriculture Analysis\nSentinel-2 L2A · Jalisco, Mexico · 2024-05-26 · 10m/px · ~245 km²',
color='white', fontsize=14, fontweight='bold', y=0.98
)
plt.savefig('../data/sprint04_multispectral_jalisco.png', dpi=150, bbox_inches='tight', facecolor='#0d0d0d')
plt.show()
print("Saved.")
Saved.
In [11]:
print("=" * 55)
print("SPRINT 04 — MULTISPECTRAL AGRICULTURE ANALYSIS")
print("=" * 55)
print(f"Scene: {best_item.id}")
print(f"Date: 2024-05-26")
print(f"Region: Jalisco, Mexico (Tequila / Amatitán)")
print(f"Source: Sentinel-2 L2A · Microsoft Planetary Computer")
print(f"CRS: EPSG:32614 (UTM Zone 14N)")
print(f"Resolution: 10m/px")
print(f"Area: ~245 km²")
print(f"Cloud Cover: 0.002%")
print("-" * 55)
print(f"NDVI mean: {np.nanmean(NDVI):.3f} | max: {np.nanmax(NDVI):.3f}")
print(f"NDWI mean: {np.nanmean(NDWI):.3f} | max: {np.nanmax(NDWI):.3f}")
print(f"NDRE mean: {np.nanmean(NDRE):.3f} | max: {np.nanmax(NDRE):.3f}")
print(f"SAVI mean: {np.nanmean(SAVI):.3f} | max: {np.nanmax(SAVI):.3f}")
print("-" * 55)
print(f"Healthy vegetation (NDVI > 0.4): "
f"{(NDVI[valid] > 0.4).sum() / valid.sum() * 100:.1f}% of valid pixels")
print(f"Water stress (NDWI > 0): "
f"{(NDWI[valid] > 0.0).sum() / valid.sum() * 100:.1f}% of valid pixels")
print("=" * 55)
======================================================= SPRINT 04 — MULTISPECTRAL AGRICULTURE ANALYSIS ======================================================= Scene: S2A_MSIL2A_20240526T171901_R012_T13QFD_20240527T012538 Date: 2024-05-26 Region: Jalisco, Mexico (Tequila / Amatitán) Source: Sentinel-2 L2A · Microsoft Planetary Computer CRS: EPSG:32614 (UTM Zone 14N) Resolution: 10m/px Area: ~245 km² Cloud Cover: 0.002% ------------------------------------------------------- NDVI mean: 0.186 | max: 0.837 NDWI mean: -0.254 | max: 0.545 NDRE mean: 0.117 | max: 0.761 SAVI mean: 0.279 | max: 1.255 ------------------------------------------------------- Healthy vegetation (NDVI > 0.4): 5.2% of valid pixels Water stress (NDWI > 0): 0.5% of valid pixels =======================================================
In [ ]: