Bagian ini mencakup: - Perhitungan luas kawasan - Analisis desa - Analisis fungsi kawasan - Visualisasi peta interaktif
# Analisis Landscape Gunung Naning & HL Gunung Naning
# Visualisasi interaktif dengan leaflet + analisis luas & desa
# Input SHP:
# - landscape_naning.shp
# - status_landscape.shp
# - Desa_landscape_naning.shp
# - HLGN.shp
# - Desa_hlgn.shp
# ---- Paket ----
pkgs <- c("sf", "dplyr", "leaflet", "stringr", "htmltools", "tibble")
need <- pkgs[!vapply(pkgs, requireNamespace, FUN.VALUE = logical(1), quietly = TRUE)]
if (length(need) > 0) install.packages(need)
library(sf)
library(dplyr)
library(leaflet)
library(stringr)
library(htmltools)
library(tibble)
# ---- Lokasi data ----
base_dir <- "C:/Users/darma/Downloads/Landscape HLGN"
files <- list(
landscape = file.path(base_dir, "landscape_naning.shp"),
status = file.path(base_dir, "status_landscape.shp"),
desa_lc = file.path(base_dir, "Desa_landscape_naning.shp"),
hlgn = file.path(base_dir, "HLGN.shp"),
desa_hlgn = file.path(base_dir, "Desa_hlgn.shp")
)
# ---- Fungsi bantu ----
read_valid <- function(x) {
if (!file.exists(x)) stop("File tidak ditemukan: ", x)
out <- st_read(x, quiet = TRUE)
out <- st_make_valid(out)
out
}
name_clean <- function(x) {
x <- as.character(x)
str_to_upper(str_squish(x))
}
area_ha <- function(x) {
as.numeric(st_area(x)) / 10000
}
# ---- Baca data ----
landscape <- read_valid(files$landscape)
status <- read_valid(files$status)
desa_lc <- read_valid(files$desa_lc)
hlgn <- read_valid(files$hlgn)
desa_hlgn <- read_valid(files$desa_hlgn)
# ---- Samakan CRS untuk analisis luas ----
# Kalimantan Barat umumnya aman memakai UTM Zone 49S (EPSG:32749)
# Jika data Anda berada di zona lain, ganti jika diperlukan.
crs_area <- 32749
landscape_a <- st_transform(landscape, crs_area)
status_a <- st_transform(status, crs_area)
hlgn_a <- st_transform(hlgn, crs_area)
desa_lc_a <- st_transform(desa_lc, crs_area)
desa_hlgn_a <- st_transform(desa_hlgn, crs_area)
# ---- Luas Landscape & HL Gunung Naning ----
luas_landscape_ha <- sum(area_ha(landscape_a), na.rm = TRUE)
luas_hlgn_ha <- sum(area_ha(hlgn_a), na.rm = TRUE)
area_summary <- tibble(
objek = c("Landscape Gunung Naning", "HL Gunung Naning"),
luas_ha = c(luas_landscape_ha, luas_hlgn_ha)
)
print("=== LUAS AREA ===")## [1] "=== LUAS AREA ==="print(area_summary)## # A tibble: 2 × 2
## objek luas_ha
## <chr> <dbl>
## 1 Landscape Gunung Naning 363341.
## 2 HL Gunung Naning 231734.# ---- Jumlah desa & desa yang sama di dua kategori ----
# Ambil nama desa unik berdasarkan atribut NAMOBJ
desa_lc_u <- desa_lc %>%
st_drop_geometry() %>%
mutate(NAMOBJ_CLEAN = name_clean(NAMOBJ)) %>%
filter(!is.na(NAMOBJ_CLEAN), NAMOBJ_CLEAN != "") %>%
distinct(NAMOBJ_CLEAN, NAMOBJ)
desa_hlgn_u <- desa_hlgn %>%
st_drop_geometry() %>%
mutate(NAMOBJ_CLEAN = name_clean(NAMOBJ)) %>%
filter(!is.na(NAMOBJ_CLEAN), NAMOBJ_CLEAN != "") %>%
distinct(NAMOBJ_CLEAN, NAMOBJ)
desa_only_lc <- setdiff(desa_lc_u$NAMOBJ_CLEAN, desa_hlgn_u$NAMOBJ_CLEAN)
desa_only_hlgn <- setdiff(desa_hlgn_u$NAMOBJ_CLEAN, desa_lc_u$NAMOBJ_CLEAN)
desa_both <- intersect(desa_lc_u$NAMOBJ_CLEAN, desa_hlgn_u$NAMOBJ_CLEAN)
desa_both_table <- tibble(
NAMOBJ = desa_both
) %>%
left_join(desa_lc_u %>% select(NAMOBJ_CLEAN, NAMOBJ) %>% rename(NAMOBJ_LANDSCAPE = NAMOBJ),
by = c("NAMOBJ" = "NAMOBJ_CLEAN")) %>%
left_join(desa_hlgn_u %>% select(NAMOBJ_CLEAN, NAMOBJ) %>% rename(NAMOBJ_HLGN = NAMOBJ),
by = c("NAMOBJ" = "NAMOBJ_CLEAN"))
desa_summary <- tibble(
kategori = c("Desa dalam landscape", "Desa dalam HLGN", "Desa ditemukan di keduanya"),
jumlah = c(nrow(desa_lc_u), nrow(desa_hlgn_u), length(desa_both))
)
print("=== RINGKASAN DESA ===")## [1] "=== RINGKASAN DESA ==="print(desa_summary)## # A tibble: 3 × 2
## kategori jumlah
## <chr> <int>
## 1 Desa dalam landscape 45
## 2 Desa dalam HLGN 26
## 3 Desa ditemukan di keduanya 25print("=== DESA YANG ADA DI KEDUA KATEGORI ===")## [1] "=== DESA YANG ADA DI KEDUA KATEGORI ==="print(desa_both_table)## # A tibble: 25 × 3
## NAMOBJ NAMOBJ_LANDSCAPE NAMOBJ_HLGN
## <chr> <chr> <chr>
## 1 SINAR PEKAYAU Sinar Pekayau Sinar Pekayau
## 2 SUNGAI SEGAK Sungai Segak Sungai Segak
## 3 SUNSONG Sunsong Sunsong
## 4 BALAI LAGAS Balai Lagas Balai Lagas
## 5 BERNAYAU Bernayau Bernayau
## 6 GANJANG Ganjang Ganjang
## 7 HARAPAN JAYA Harapan Jaya Harapan Jaya
## 8 KARANG BETUNG Karang Betung Karang Betung
## 9 KELUAS HULU Keluas Hulu Keluas Hulu
## 10 KEMANTAN Kemantan Kemantan
## # ℹ 15 more rows# ---- 7. Analisis 3: Luas masing-masing status hutan/lahan berdasarkan Fungsi ----
if (!("Fungsi" %in% names(status_a))) {
stop("Kolom 'Fungsi' tidak ditemukan pada status_landscape.shp")
}
status_summary <- status_a %>%
mutate(luas_ha = area_ha(geometry)) %>%
st_drop_geometry() %>%
group_by(Fungsi) %>%
summarise(
luas_ha = sum(luas_ha, na.rm = TRUE),
jumlah_poligon = n(),
.groups = "drop"
) %>%
arrange(desc(luas_ha))
print("=== LUAS PER FUNGSI ===")## [1] "=== LUAS PER FUNGSI ==="print(status_summary)## # A tibble: 6 × 3
## Fungsi luas_ha jumlah_poligon
## <chr> <dbl> <int>
## 1 HL 208622. 10
## 2 HPT 134346. 5
## 3 HP 15279. 5
## 4 APL 4300. 1
## 5 HPK 774. 2
## 6 Air 27.3 1# ---- Simpan hasil ke CSV ----
write.csv(area_summary, file.path(base_dir, "ringkasan_luas_area.csv"), row.names = FALSE)
write.csv(desa_summary, file.path(base_dir, "ringkasan_desa.csv"), row.names = FALSE)
write.csv(desa_both_table, file.path(base_dir, "desa_yang_ada_di_keduanya.csv"), row.names = FALSE)
write.csv(status_summary, file.path(base_dir, "ringkasan_fungsi_status.csv"), row.names = FALSE)
# ===============================
# 1. Load paket
# ===============================
library(sf)
library(leaflet)
library(dplyr)
library(htmlwidgets)
# ===============================
# 2. Path data
# ===============================
base_dir <- "C:/Users/darma/Downloads/Landscape HLGN"
landscape <- st_read(file.path(base_dir, "landscape_naning.shp"), quiet = TRUE)
status <- st_read(file.path(base_dir, "status_landscape.shp"), quiet = TRUE)
desa_lc <- st_read(file.path(base_dir, "Desa_landscape_naning.shp"), quiet = TRUE)
hlgn <- st_read(file.path(base_dir, "HLGN.shp"), quiet = TRUE)
desa_hlgn <- st_read(file.path(base_dir, "Desa_hlgn.shp"), quiet = TRUE)
# ===============================
# 3. Validasi + CRS
# ===============================
landscape <- st_transform(st_make_valid(landscape), 4326)
status <- st_transform(st_make_valid(status), 4326)
desa_lc <- st_transform(st_make_valid(desa_lc), 4326)
hlgn <- st_transform(st_make_valid(hlgn), 4326)
desa_hlgn <- st_transform(st_make_valid(desa_hlgn), 4326)
# ===============================
# 4. Simplify (BIAR RINGAN)
# ===============================
landscape <- st_simplify(landscape, dTolerance = 0.005)
status <- st_simplify(status, dTolerance = 0.005)
desa_lc <- st_simplify(desa_lc, dTolerance = 0.005)
hlgn <- st_simplify(hlgn, dTolerance = 0.005)
desa_hlgn <- st_simplify(desa_hlgn, dTolerance = 0.005)
# ===============================
# 5. Warna
# ===============================
pal_status <- colorFactor("Set2", domain = status$Fungsi)
# ===============================
# 6. Titik tengah
# ===============================
center <- st_coordinates(st_centroid(st_union(landscape)))
# ===============================
# 7. GROUP (BIAR KONSISTEN)
# ===============================
grp_osm <- "🌍 OpenStreetMap"
grp_sat <- "🛰️ Satelit"
grp_lc_batas <- "🌳 LANDSCAPE | Batas"
grp_lc_status <- "🌳 LANDSCAPE | Status"
grp_lc_desa <- "🌳 LANDSCAPE | Desa"
grp_hl_batas <- "🏞️ HL | Batas"
grp_hl_desa <- "🏞️ HL | Desa"
# ===============================
# 8. Leaflet Map
# ===============================
m <- leaflet() %>%
# ===== BASEMAP =====
addProviderTiles(providers$OpenStreetMap, group = grp_osm) %>%
addProviderTiles(providers$Esri.WorldImagery, group = grp_sat) %>%
# ===== LANDSCAPE =====
addPolygons(
data = landscape,
color = "red",
weight = 2,
fill = FALSE,
group = grp_lc_batas
) %>%
addPolygons(
data = status,
fillColor = ~pal_status(Fungsi),
color = "white",
weight = 1,
fillOpacity = 0.5,
popup = ~paste("Fungsi:", Fungsi),
group = grp_lc_status
) %>%
addPolygons(
data = desa_lc,
fillColor = "blue",
color = "blue",
weight = 1,
fillOpacity = 0.2,
popup = ~NAMOBJ,
group = grp_lc_desa
) %>%
# ===== HL =====
addPolygons(
data = hlgn,
color = "green",
weight = 2,
fill = FALSE,
group = grp_hl_batas
) %>%
addPolygons(
data = desa_hlgn,
fillColor = "orange",
color = "orange",
weight = 1,
fillOpacity = 0.2,
popup = ~NAMOBJ,
group = grp_hl_desa
) %>%
# ===== LEGEND =====
addLegend(
position = "bottomright",
pal = pal_status,
values = status$Fungsi,
title = "Fungsi Hutan"
) %>%
# ===== CONTROL (Satelit di atas) =====
addLayersControl(
baseGroups = c(grp_sat, grp_osm), # ✅ urutan diperbaiki
overlayGroups = c(
grp_lc_batas,
grp_lc_status,
grp_lc_desa,
grp_hl_batas,
grp_hl_desa
),
options = layersControlOptions(collapsed = FALSE)
) %>%
# ===== DEFAULT: HL OFF =====
# ===== DEFAULT: hanya batas landscape aktif =====
hideGroup(c(
grp_lc_status,
grp_lc_desa,
grp_hl_batas,
grp_hl_desa
)) %>%
# ===== VIEW =====
setView(lng = center[1], lat = center[2], zoom = 9)
# ===============================
# 9. Tampilkan
# ===============================
m# ===============================
# 10. Export ke HTML
# ===============================
saveWidget(
m,
file = file.path(base_dir, "peta_landscape_hlgn_FINAL.html"),
selfcontained = TRUE
)