Skip to content

Instantly share code, notes, and snippets.

@dinizime
Created February 28, 2026 23:22
Show Gist options
  • Select an option

  • Save dinizime/ae6eac0863fdb89b3b6f9b7b3f3409d6 to your computer and use it in GitHub Desktop.

Select an option

Save dinizime/ae6eac0863fdb89b3b6f9b7b3f3409d6 to your computer and use it in GitHub Desktop.
Faz o stack de 4 imagens em lote, reprojeta para 3857 e pixel de 9.555 com reamostragem bilinear
#!/usr/bin/env python3
"""
Stack e reprojeção de imagens Sentinel-2 TCI para Super Resolution.
Replica o model_sem_corte do QGIS:
1. gdal_merge -separate (stack temporal como bandas)
2. gdalwarp para EPSG:3857, bilinear, pixel 9.555m
"""
from osgeo import gdal
from pathlib import Path
import re
import sys
import tempfile
import os
gdal.UseExceptions()
TARGET_SRS = "EPSG:3857"
PIXEL_SIZE = 9.555
RESAMPLE = "bilinear"
EXTENSIONS = (".tif", ".tiff", ".jp2")
def parse_sentinel(name):
"""Extrai tile e data de um filename Sentinel-2 TCI."""
m = re.search(r"(T\d{2}[A-Z]{3})_(\d{8})T\d{6}_TCI", name)
return (m.group(1), m.group(2)) if m else (None, None)
def process_folder(folder):
folder = Path(folder)
images = sorted(f for f in folder.iterdir() if f.suffix.lower() in EXTENSIONS)
if len(images) != 4:
print(f" SKIP {folder.name}: {len(images)} imagens (precisa 4)")
return False
parsed = []
for img in images:
tile, date = parse_sentinel(img.name)
if not tile:
print(f" ERRO nao consegui parsear: {img.name}")
return False
parsed.append((date, tile, img))
parsed.sort()
tiles = set(p[1] for p in parsed)
if len(tiles) != 1:
print(f" ERRO tiles misturados: {tiles}")
return False
tile = parsed[0][1]
date_min = parsed[0][0]
date_max = parsed[-1][0]
output_path = folder / f"{tile}_{date_min}_{date_max}.tif"
if output_path.exists():
print(f" SKIP {output_path.name} ja existe")
return True
print(f" {folder.name} -> {output_path.name}")
images_sorted = [p[2] for p in parsed]
with tempfile.TemporaryDirectory() as tmp:
# Expandir imagens multi-banda em VRTs single-band
# Replica gdal_merge -separate: img1_R, img1_G, img1_B, img2_R, ...
band_vrts = []
for img in images_sorted:
ds = gdal.Open(str(img))
nb = ds.RasterCount
ds = None
for b in range(1, nb + 1):
vrt = os.path.join(tmp, f"{img.stem}_b{b}.vrt")
gdal.Translate(vrt, str(img), bandList=[b], format="VRT")
band_vrts.append(vrt)
# Stack todas as bandas
stack_vrt = os.path.join(tmp, "stack.vrt")
ds = gdal.BuildVRT(
stack_vrt, band_vrts, options=gdal.BuildVRTOptions(separate=True)
)
ds.FlushCache()
ds = None
# Reprojetar para EPSG:3857 com pixel exato de 9.555m bilinear
ds = gdal.Warp(
str(output_path),
stack_vrt,
options=gdal.WarpOptions(
dstSRS=TARGET_SRS,
xRes=PIXEL_SIZE,
yRes=PIXEL_SIZE,
resampleAlg=RESAMPLE,
outputType=gdal.GDT_Byte,
format="GTiff",
),
)
ds.FlushCache()
ds = None
size_mb = output_path.stat().st_size / 1e6
print(f" OK {output_path.name} ({size_mb:.1f} MB)")
return True
def main():
if len(sys.argv) < 2:
print("Uso: python stack_sr.py <pasta_raiz>")
print(" Navega em cada subpasta, empilha 4 TCI e reprojeta para EPSG:3857")
sys.exit(1)
root = Path(sys.argv[1])
if not root.is_dir():
print(f"Erro: {root} nao e um diretorio")
sys.exit(1)
folders = sorted(d for d in root.iterdir() if d.is_dir())
if not folders:
print("Nenhuma subpasta encontrada")
sys.exit(1)
print(f"Encontradas {len(folders)} subpastas em {root}\n")
ok = err = 0
for f in folders:
if process_folder(f):
ok += 1
else:
err += 1
print(f"\nFinalizado: {ok} ok, {err} skip/erro")
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment