# -*- coding: utf-8 -*-
# SKRIP 2.1: Buang awan dan bayangan memakai band SCL, per tanggal, dengan Raster Calculator
# Penulis: Badar Mubarok Yogaswara
# Syarat: Skrip 1.1 sudah dijalankan (hasil/ndvi/NDVI_*.tif ada). Piksel buruk menjadi NoData.
import os
import processing
from qgis.core import QgsRasterLayer, QgsProject
from osgeo import gdal
import numpy as np
import m3_umum as U

os.makedirs(os.path.join(U.HASIL, "ndvi_bersih"), exist_ok=True)
daftar = U.daftar_citra()
keluaran = []
for tgl, pr, ps in daftar:
    k = tgl.strftime("%Y%m%d")
    n_ndvi, n_scl = "NDVI_" + k, "SCL_" + k
    l_ndvi = QgsRasterLayer(os.path.join(U.HASIL, "ndvi", "NDVI_%s.tif" % k), n_ndvi)
    l_scl = QgsRasterLayer(ps, n_scl)
    if not l_ndvi.isValid() or not l_scl.isValid():
        raise SystemExit("Layer tidak terbaca (jalankan Skrip 1.1 dulu): " + n_ndvi)
    QgsProject.instance().addMapLayers([l_ndvi, l_scl], False)
    # piksel baik: SCL 4 (vegetasi), 5 (bukan vegetasi), 6 (air). Nilai baik dibagi 1; nilai buruk dibagi 0 -> NoData
    ekspresi = '"%s@1" / (("%s@1" = 4) OR ("%s@1" = 5) OR ("%s@1" = 6))' % (n_ndvi, n_scl, n_scl, n_scl)
    out = os.path.join(U.HASIL, "ndvi_bersih", "NDVI_bersih_%s.tif" % k)
    processing.run("native:rastercalc", {"LAYERS": [l_ndvi, l_scl], "EXPRESSION": ekspresi, "OUTPUT": out})
    keluaran.append(out)
    a = gdal.Open(out)
    arr = a.GetRasterBand(1).ReadAsArray()
    nd = a.GetRasterBand(1).GetNoDataValue()
    valid = np.isfinite(arr) if nd is None else (arr != nd) & np.isfinite(arr)
    print(tgl, "piksel valid: %5.1f%%" % (100 * valid.mean()))
print("selesai,", len(keluaran), "berkas")