# [SKRIP 1.1: Memuat COG & Eksekusi NDVI Virtual di QGIS 4.x]
# Penulis: Badar Mubarok Yogaswara
# Syarat: Pastikan layer batas area (misal: 'Batas_KPH') aktif di QGIS.

import processing
from qgis.core import QgsProject, QgsRasterLayer

print("🔗 Memulai koneksi ke Cloud-Optimized GeoTIFF (COG)...")

# 1. Tentukan URL COG dari STAC API (Contoh URL dummy Sentinel-2)
url_b4 = "/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/2026/T48M/B04.tif"
url_b8 = "/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/2026/T48M/B08.tif"

# 2. Muat Layer secara Virtual
layer_b4 = QgsRasterLayer(url_b4, "B04_Red_Cloud", "gdal")
layer_b8 = QgsRasterLayer(url_b8, "B08_NIR_Cloud", "gdal")

if not layer_b4.isValid() or not layer_b8.isValid():
    print("❌ Gagal memuat citra dari Cloud. Cek koneksi internet.")
else:
    QgsProject.instance().addMapLayers([layer_b4, layer_b8])
    print("✅ Citra Red dan NIR berhasil dimuat.")

    # 3. Ambil layer Vektor untuk batas pemotongan
    batas_kph = QgsProject.instance().mapLayersByName('Batas_KPH')[0]

    # 4. Potong (Clip) Citra langsung dari Cloud menggunakan Algoritma Native GDAL
    print("✂️ Memotong citra berdasarkan batas KPH...")

    # Clip Band 8 (NIR)
    clip_b8 = processing.run("gdal:cliprasterbymasklayer", {
        'INPUT': layer_b8,
        'MASK': batas_kph,
        'CROP_TO_CUTLINE': True,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    })['OUTPUT']

    # Clip Band 4 (Red)
    clip_b4 = processing.run("gdal:cliprasterbymasklayer", {
        'INPUT': layer_b4,
        'MASK': batas_kph,
        'CROP_TO_CUTLINE': True,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    })['OUTPUT']

    # 5. Kalkulasi NDVI Menggunakan Raster Calculator Native
    # Rumus: (NIR - RED) / (NIR + RED)
    print("📊 Menghitung NDVI...")

    ndvi_result = processing.run("qgis:rastercalculator", {
        'EXPRESSION': '("B08@1" - "B04@1") / ("B08@1" + "B04@1")',
        'LAYERS': [clip_b8, clip_b4],
        'CELLSIZE': 10,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    })

    QgsProject.instance().addMapLayer(QgsRasterLayer(ndvi_result['OUTPUT'], "NDVI_KPH", "gdal"))
    print("🚀 Selesai! Layer NDVI_KPH telah ditambahkan ke kanvas.")