library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(tmap)
library(leaflet)
library(classInt)
library(RColorBrewer)
library(htmltools)
geojson_path <- "C:/Users/USER/OneDrive/Dokumen/Kulyah smt 5/Statistika Lingkungan/33.26_kelurahan.geojson"
kel <- st_read(geojson_path, quiet = TRUE)
print(kel)
## Simple feature collection with 285 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 109.4834 ymin: -7.244428 xmax: 109.7991 ymax: -6.840548
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
## First 10 features:
## kd_propinsi kd_dati2 kd_kecamatan kd_kelurahan nm_kelurahan
## 1 33 26 015 016 Curug
## 2 33 26 015 012 Karangjompo
## 3 33 26 008 007 Sinangohprendeng
## 4 33 26 008 019 Tanjungsari
## 5 33 26 001 010 Luragung
## 6 33 26 009 016 Ponolawen
## 7 33 26 008 005 Pringsurat
## 8 33 26 006 001 Pungangan
## 9 33 26 011 005 Randumuktiwaren
## 10 33 26 006 009 Randusari
## geometry
## 1 MULTIPOLYGON Z (((109.6439 ...
## 2 MULTIPOLYGON Z (((109.6444 ...
## 3 MULTIPOLYGON Z (((109.5915 ...
## 4 MULTIPOLYGON Z (((109.6026 ...
## 5 MULTIPOLYGON Z (((109.5032 ...
## 6 MULTIPOLYGON Z (((109.5209 ...
## 7 MULTIPOLYGON Z (((109.5916 ...
## 8 MULTIPOLYGON Z (((109.6718 ...
## 9 MULTIPOLYGON Z (((109.559 -...
## 10 MULTIPOLYGON Z (((109.7099 ...
cat(sprintf("Jumlah fitur: %d\n", nrow(kel)))
## Jumlah fitur: 285
# tampilkan nama kolom
attrs <- names(kel)
attrs
## [1] "kd_propinsi" "kd_dati2" "kd_kecamatan" "kd_kelurahan" "nm_kelurahan"
## [6] "geometry"
# cari kolom numerik (tanpa geometry)
num_cols <- kel %>% st_drop_geometry() %>% select(where(is.numeric)) %>% names()
cat("Kolom numerik yang tersedia:\n")
## Kolom numerik yang tersedia:
print(num_cols)
## character(0)
# jika tidak ada kolom numerik, buat kolom dummy (contoh: luas dalam m2)
if(length(num_cols) == 0){
kel$area_m2 <- as.numeric(st_area(kel))
num_cols <- "area_m2"
cat("Tidak ada kolom numerik: membuat kolom 'area_m2' berdasarkan luas geometris (m^2).\n")
}
## st_as_s2(): dropping Z and/or M coordinate
## Tidak ada kolom numerik: membuat kolom 'area_m2' berdasarkan luas geometris (m^2).
# Pilih kolom numerik terbaik untuk choropleth: gunakan yang memiliki variansi terbesar
if(length(num_cols) >= 1){
variances <- sapply(num_cols, function(cn) var(kel[[cn]], na.rm = TRUE))
target_col <- names(which.max(variances))
cat(sprintf("Kolom terpilih untuk peta tematik: %s\n", target_col))
} else {
stop("Tidak ada kolom numerik dan pembuatan kolom area gagal.")
}
## Kolom terpilih untuk peta tematik: area_m2
nclass <- 5
vals <- kel[[target_col]]
vals_numeric <- as.numeric(vals)
vals_numeric[is.na(vals_numeric)] <- median(vals_numeric, na.rm = TRUE)
ci <- classIntervals(vals_numeric, n = nclass, style = "quantile")
breaks <- as.numeric(ci$brks)
# Gunakan cut untuk label agar tidak error
labs <- levels(cut(vals_numeric, breaks = breaks, include.lowest = TRUE))
kel_fort <- st_transform(kel, 3857)
kel_fort$cat <- cut(as.numeric(kel[[target_col]]), breaks = breaks, include.lowest = TRUE, labels = labs)
ggplot() +
geom_sf(data = kel_fort, aes(fill = cat), color = "grey20", size = 0.2) +
labs(title = paste("Peta Tematik:", target_col), fill = target_col) +
theme_minimal() +
theme(axis.text = element_blank(), axis.ticks = element_blank())
tmap_mode("view")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
m_tmap <- tm_shape(kel) +
tm_polygons(col = target_col, style = "quantile", n = nclass, palette = "YlOrRd", title = target_col, alpha = 0.8) +
tm_layout(title = paste("Peta Choropleth -", target_col))
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
m_tmap
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
## "brewer.yl_or_rd"
## Multiple palettes called "yl_or_rd" found: "brewer.yl_or_rd", "matplotlib.yl_or_rd". The first one, "brewer.yl_or_rd", is returned.
library(sf)
library(leaflet)
# 1. Baca geojson
kel <- st_read("C:/Users/USER/OneDrive/Dokumen/Kulyah smt 5/Statistika Lingkungan/33.26_kelurahan.geojson")
## Reading layer `33.26_kelurahan' from data source
## `C:\Users\USER\OneDrive\Dokumen\Kulyah smt 5\Statistika Lingkungan\33.26_kelurahan.geojson'
## using driver `GeoJSON'
## Simple feature collection with 285 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 109.4834 ymin: -7.244428 xmax: 109.7991 ymax: -6.840548
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
# 2. Pastikan geometry valid + CRS sesuai
kel <- st_make_valid(kel)
kel <- st_transform(kel, crs = 4326)
# 3. Tentukan kolom target
target_col <- "area_m2"
vals_numeric <- kel[[target_col]]
# 4. Definisikan palet warna
pal_leaf <- colorBin(
palette = "YlOrRd",
domain = vals_numeric,
bins = 5,
na.color = "#808080"
)
# 5. Buat popup
popup_attrs <- names(kel)[1:min(4, ncol(kel))]
popup <- apply(
st_drop_geometry(kel)[, popup_attrs, drop = FALSE],
1,
function(r) paste(
paste0(popup_attrs, ": ", ifelse(is.na(r), "-", r)),
collapse = "<br/>"
)
)
# 6. Plot choropleth
leaflet(kel) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
fillColor = ~pal_leaf(vals_numeric),
weight = 0.5,
opacity = 1,
color = "#444444",
fillOpacity = 0.8,
highlight = highlightOptions(weight = 2, color = "#666", bringToFront = TRUE),
label = lapply(popup, HTML)
) %>%
addLegend(
pal = pal_leaf,
values = vals_numeric,
title = "Luas Area (m²)",
opacity = 1
)
# Hitung jumlah kelurahan per kecamatan
summary_table <- kel %>%
st_drop_geometry() %>%
group_by(kd_kecamatan) %>%
summarise(jumlah_kelurahan = n()) %>%
arrange(kd_kecamatan)
summary_table
## # A tibble: 19 × 2
## kd_kecamatan jumlah_kelurahan
## <chr> <int>
## 1 001 14
## 2 002 15
## 3 003 11
## 4 004 9
## 5 005 10
## 6 006 14
## 7 007 15
## 8 008 25
## 9 009 23
## 10 010 17
## 11 011 22
## 12 012 14
## 13 013 19
## 14 014 10
## 15 015 16
## 16 016 16
## 17 017 13
## 18 018 11
## 19 019 11