# -*- coding: utf-8 -*-
# I2 Bab 6 (latihan terpadu): lahan terbuka di lereng curam, dari peta kelas (1 m) dan kemiringan (2 m).
# Jalankan setelah skrip bab 2 dan bab 5. Penulis: Badar Mubarok Yogaswara
import os
import numpy as np
import processing
from osgeo import gdal
from qgis.core import QgsRasterLayer

HASIL = r"C:/temp/hasil_i2"
kelas = QgsRasterLayer(os.path.join(HASIL, "terbimbing_maks_lik.tif"), "peta_kelas")
persen = QgsRasterLayer(os.path.join(HASIL, "slope_persen.tif"), "slope_persen")

# Dua raster berbeda ukuran piksel (1 m dan 2 m). Raster calculator memakai ukuran dan luas yang Anda tentukan,
# di sini mengikuti peta kelas: piksel 1 m, area 300 m x 300 m. Nilai kemiringan diambil dari piksel 2 m yang menutupinya.
AMBANG = 25                       # batas kemiringan (persen)
KELAS_UJI = 4                     # kelas yang dicari (4 = Lahan Terbuka)
NAMA_KELAS = {1: "hutan", 2: "kebun", 3: "sawah", 4: "lahan terbuka"}[KELAS_UJI]
luas_ha = kelas.rasterUnitsPerPixelX() * kelas.rasterUnitsPerPixelY() / 10000.0     # luas satu piksel (ha)
keluar = os.path.join(HASIL, "terbuka_curam.tif")
processing.run("native:rastercalc", {
    "EXPRESSION": '("peta_kelas@1" = %d) AND ("slope_persen@1" > %d)' % (KELAS_UJI, AMBANG), "LAYERS": [kelas, persen],
    "CELL_SIZE": kelas.rasterUnitsPerPixelX(), "EXTENT": kelas.extent(), "CRS": kelas.crs(), "OUTPUT": keluar})
arr = gdal.Open(keluar).ReadAsArray()
print("Piksel %s di lereng > %d persen:" % (NAMA_KELAS, AMBANG), int(arr.sum()), "= %.3f ha" % (arr.sum() * luas_ha))
terbuka = gdal.Open(os.path.join(HASIL, "terbimbing_maks_lik.tif")).ReadAsArray() == KELAS_UJI
print("Seluruh %s pada peta: %d piksel = %.3f ha; bagian di lereng curam: %.1f%%" % (
    NAMA_KELAS, terbuka.sum(), terbuka.sum() * luas_ha, 100.0 * arr.sum() / terbuka.sum()))