From GEE to just Python

image = ee.Image("COPERNICUS/S2/...")

import rasterio
from pathlib import Path

# File path
image_path = Path("test_image.tif")

# Open image: Stacked GeoTIFF with multiple bands
with rasterio.open(image_path) as src:
    # Metadata
    print(src.crs)          # CRS
    print(src.bounds)       # Bounding box
    print(src.shape)        # (height, width)
    print(src.count)        # Number of bands
    print(src.transform)    # Affine transform

rgb = image.select(["B4", "B3", "B2"])
nir = image.select("B8")

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

# Open the image and read bands
with rasterio.open("test_image.tif") as src:
    all_bands = src.read()       # Read all bands at once
    band1 = src.read(1)          # Optional: read only first band (index inside .read())

# Assign Sentinel-2 bands
B2 = all_bands[1]  # Blue
B3 = all_bands[2]  # Green
B4 = all_bands[3]  # Red
B8 = all_bands[7]  # NIR

# Stack RGB for display
rgb = np.stack([B4, B3, B2], axis=2)  # shape: (height, width, 3)

# Optional: Display RGB image
plt.figure(figsize=(8, 8))
plt.imshow(rgb / np.max(rgb))  # normalize for display
plt.title("RGB Image (Red, Green, Blue)")
plt.axis("off")
plt.show()

# Optional: Display individual bands
bands = [B2, B3, B4, B8]
titles = ["Blue (B2)", "Green (B3)", "Red (B4)", "NIR (B8)"]

fig, axes = plt.subplots(1, 4, figsize=(20, 5))
for ax, band, title in zip(axes, bands, titles):
    ax.imshow(band, cmap="gray")
    ax.set_title(title)
    ax.axis("off")
plt.show()

ndvi = image.normalizedDifference(["B8", "B4"]).rename("NDVI")

ndvi = (B8 - B4) / (B8 + B4)

# Simple trick to avoid division by zero
ndvi = (B8 - B4) / (B8 + B4 + 1e-10)

nir = image.select("B8")
red = image.select("B4")
ndvi = nir.subtract(red).divide(nir.add(red))

nir = B8.astype(float)
red = B4.astype(float)

ndvi = (nir - red) / (nir + red)

# Display NDVI
plt.figure(figsize=(8, 8))
plt.imshow(ndvi, cmap="RdYlGn")  # Green = high NDVI, red = low NDVI
plt.colorbar(label="NDVI")
plt.title("NDVI Image")
plt.axis("off")
plt.show()

import rasterio
from rasterio.windows import Window
import matplotlib.pyplot as plt
from pathlib import Path

# Define window: col_off, row_off, width, height
window = Window(col_off=0, row_off=0, width=512, height=512)

# Read Windowed Subset
with rasterio.open(image_path) as src:
    subset = src.read(1, window=window)  # Read band 1 in the window

# Display Subset
plt.imshow(subset, cmap="gray")
plt.title("Windowed subset of Band 1")
plt.axis("off")
plt.show()

import rasterio                            # Read/write GeoTIFFs
import xarray as xr                        # Handle multi-band raster arrays
import numpy as np                         # Perform band math and array operations
import matplotlib.pyplot as plt            # Visualize results
from skimage.transform import resize       # Resize images
from pathlib import Path                   # Handle file paths

from rasterio.plot import show
from rasterio.warp import calculate_default_transform, reproject, Resampling

s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
        .filterBounds(aoi_region)
        .filterDate('2024-01-01', '2024-12-31')
        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', maxCloud))
        .median().clip(aoi_region)
        .select(['B2','B3','B4','B8','B11']))
    
# Landsat 8/9
l8 = (ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
        .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
        .filterBounds(aoi_region)
        .filterDate('2024-01-01', '2024-12-31')
        .filter(ee.Filter.lt('CLOUD_COVER', maxCloud))
        .median().clip(aoi_region)
        .multiply(0.0000275).subtract(0.2)
        .select(['SR_B2','SR_B3','SR_B4','SR_B5'])
        .rename(['B2','B3','B4','B8']))

# Load all data
data_folder = Path("data/S2_LS_test_data")

# Sentinel-2 bands - list of S2 files and read each band as array
s2_files = list(data_folder.glob("*S2*.tif"))                                 
s2_bands = [rasterio.open(f).read(1).astype(float) for f in s2_files]         

# Landsat 8 bands (scaled) - read and scale
l8_bands = [rasterio.open(f).read(1).astype(float)*0.0000275-0.2
            for f in data_folder.glob("*L8*.tif")]                            

# DEM - read DEM as array
dem = rasterio.open(data_folder / "DEM_elevation.tif").read(1).astype(float)  

# Combine S2 and L8
combined = s2.addBands(l8)

# Get reference size and georeference from first Sentinel-2 file
with rasterio.open(s2_files[0]) as ref:
    shape = ref.read(1).shape            # image dimensions (rows, cols)
    crs = ref.crs                        # coordinate reference system
    transform = ref.transform            # map location of top-left pixel
    bounds = ref.bounds

# Resize Landsat 8 and DEM to match Sentinel-2 
l8_resized = [resize(b, shape, preserve_range=True) for b in l8_bands]
dem_resized = resize(dem, shape, preserve_range=True)

# Combine all bands - 
all_bands = s2_bands + l8_resized + [dem_resized]          
stacked = np.stack(all_bands)                        

# Save as one georeferenced multiband TIFF
with rasterio.open("S2_LS_DEM_stack.tif", "w",
                   driver="GTiff",                      # file format
                   height=shape[0], width=shape[1],     # image size
                   count=stacked.shape[0],              # number of bands
                   dtype=stacked.dtype,                 # data type
                   crs=crs,                             # coordinate system
                   transform=transform) as dst:         # map location
    for i, band in enumerate(stacked, 1):               # write each band
        dst.write(band, i)

print("Saved: S2_LS_DEM_stack.tif") 
`

# Add DEM & Slope from existing ee resources
dem = ee.Image('USGS/SRTMGL1_003').clip(aoi_region)
slope = ee.Terrain.slope(dem).rename('slope')
combined = combined.addBands([dem.rename('DEM'), slope])

# Use DEM from stack (already resized to match Sentinel-2)
dem_band = stacked[-1]

# Reproject DEM to UTM (meters) for accurate slope
utm_crs = "EPSG:32651"
dem_utm = np.empty_like(dem_band)
utm_transform, _, _ = calculate_default_transform(crs, 
                                                  utm_crs, 
                                                  dem_band.shape[1], 
                                                  dem_band.shape[0], 
                                                  *bounds)

reproject(
    source=dem_band,
    destination=dem_utm,
    src_transform=transform,
    src_crs=crs,
    dst_transform=utm_transform,
    dst_crs=utm_crs,
    resampling=Resampling.bilinear,
    dst_shape=dem_band.shape 
)

# Compute slope in degrees
pixel_size = utm_transform.a
grad_y, grad_x = np.gradient(dem_utm, pixel_size, pixel_size)
slope_degrees_utm = np.arctan(np.sqrt(grad_x**2 + grad_y**2)) * 180/np.pi

# Calculate Indices
ndvi = combined.normalizedDifference(['B8','B4']).rename('NDVI')
ndwi = s2.normalizedDifference(['B3','B8']).rename('NDWI')
ndbi = s2.normalizedDifference(['B11','B8']).rename('NDBI')
    
# EVI calculation requires an expression mapping
evi = s2.expression(
    '2.5*((NIR-RED)/(NIR+6*RED-7.5*BLUE+1))',
    {'NIR': s2.select('B8'), 'RED': s2.select('B4'), 'BLUE': s2.select('B2')}
).rename('EVI')
    
combined = combined.addBands([ndvi, ndwi, ndbi, evi])

# Adjust to match the order of your stacked bands
band_names = [
    "S2_B2_Blue", "S2_B3_Green", "S2_B4_Red", "S2_B8_NIR", "S2_B11_SWIR",
    "L8_B2_Blue", "L8_B3_Green", "L8_B4_Red", "L8_B5_NIR", 
    "DEM_Elevation", "Slope"
]

# Compute Indices
BLUE = stacked_with_slope[0]
GREEN = stacked_with_slope[1]
RED = stacked_with_slope[2]
NIR = stacked_with_slope[3]
SWIR = stacked_with_slope[4]

NDVI = (NIR - RED) / (NIR + RED + 1e-10)
NDWI = (GREEN - NIR) / (GREEN + NIR + 1e-10)
NDBI = (SWIR - NIR) / (SWIR + NIR + 1e-10)
EVI = 2.5 * (NIR - RED) / (NIR + 6*RED - 7.5*BLUE + 1 + 1e-10)

# Add indices to stack
stacked_with_indices = np.concatenate([
    stacked,
    NDVI[np.newaxis, :, :],
    NDWI[np.newaxis, :, :],
    NDBI[np.newaxis, :, :],
    EVI[np.newaxis, :, :]
], axis=0)

# Add names to stack
index_names = ["NDVI", "NDWI", "NDBI", "EVI"]
band_names_with_indices = band_names + index_names

# Save the stacked array with indices as a multiband TIFF
with rasterio.open(
    "S2_LS_DEM_with_indices.tif", "w",
    driver="GTiff",
    height=stacked_with_indices.shape[1],   # rows
    width=stacked_with_indices.shape[2],    # cols
    count=stacked_with_indices.shape[0],    # total bands
    dtype=stacked_with_indices.dtype,
    crs=crs,
    transform=transform
) as dst:
    
    # Write each band
    for i, band in enumerate(stacked_with_indices, start=1):
        dst.write(band, i)
        dst.set_band_description(i, band_names_with_indices[i-1])

print("Saved multiband raster with indices to: S2_LS_DEM_with_indices.tif")