spasial tugas 1
nabila
2025-10-18
DATA PRODUKTIVITAS PADI DI JAWA TIMUR
#loading package
library(readxl)
library(sp)## Warning: package 'sp' was built under R version 4.3.3library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(spdep)## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`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# import data
getwd()## [1] "C:/Users/Nabila/OneDrive/Dokumen/RStudio/SDS SMT 5"data <- read_excel('Produktivitas Padi Menurut Kabupaten_Kota di Provinsi Jawa_Timur, 2023.xlsx')
str(data)## tibble [38 × 2] (S3: tbl_df/tbl/data.frame)
## $ Kabupaten/Kota : chr [1:38] "Bangkalan" "Banyuwangi" "Blitar" "Bojonegoro" ...
## $ produktivitas_padi: num [1:38] 49 60.6 66.1 53.2 49.2 ...#import shp
indo_desa <- st_read("D:/STATISTIKA RECAP/SEMESTER 5/Sains Data Spasial/BATAS WILAYAH KELURAHAN-DESA 10K/Batas_Wilayah_KelurahanDesa_10K_AR.shp")## Reading layer `Batas_Wilayah_KelurahanDesa_10K_AR' from data source
## `D:\STATISTIKA RECAP\SEMESTER 5\Sains Data Spasial\BATAS WILAYAH KELURAHAN-DESA 10K\Batas_Wilayah_KelurahanDesa_10K_AR.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 83518 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range: zmin: 0 zmax: 1876.818
## Geodetic CRS: WGS 84# Subset langsung untuk jawa barat
jatim <- indo_desa %>% filter(WADMPR == "Jawa Timur")# Plot hasilnya
plot(st_geometry(jatim), col = "#8bab81", main = "Wilayah Jawa Barat (batas desa)")
library(lwgeom)## Warning: package 'lwgeom' was built under R version 4.3.3## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.11.2, PROJ 9.3.1##
## Attaching package: 'lwgeom'## The following objects are masked from 'package:sf':
##
## st_minimum_bounding_circle, st_perimeterjatim <- jatim %>%
st_make_valid() # perbaiki geometri yang tidak valid
sum(!st_is_valid(jatim)) # kalau hasilnya 0, semua geometri valid ✅## [1] 0# menggabungkan seluruh desa menjadi polygon kab/kota (Union/dissolve)
jatim_kab <- jatim %>%
group_by(WADMKK) %>%
summarise(geometry = st_union(geometry))
plot(st_geometry(jatim_kab),
col = "#8bab81", main = "Batas Kabupaten/Kota Wilayah Jawa Timur")
nrow(jatim_kab)## [1] 39jatim_kab$WADMKK## [1] "Bangkalan" "Banyuwangi" "Blitar" "Bojonegoro"
## [5] "Bondowoso" "Gresik" "Jember" "Jombang"
## [9] "Kediri" "Kota Batu" "Kota Blitar" "Kota Kediri"
## [13] "Kota Madiun" "Kota Malang" "Kota Mojokerto" "Kota Pasuruan"
## [17] "Kota Probolinggo" "Kota Surabaya" "Lamongan" "Lumajang"
## [21] "Madiun" "Magetan" "Malang" "Mojokerto"
## [25] "Nganjuk" "Ngawi" "Pacitan" "Pamekasan"
## [29] "Pasuruan" "Ponorogo" "Probolinggo" "Sampang"
## [33] "Sidoarjo" "Situbondo" "Sumenep" "Trenggalek"
## [37] "Tuban" "Tulungagung" NAjatim_kab[is.na(jatim_kab$WADMKK), ]## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 112.6637 ymin: -7.578155 xmax: 112.8824 ymax: -7.191753
## Geodetic CRS: WGS 84
## # A tibble: 1 × 2
## WADMKK geometry
## <chr> <MULTIPOLYGON [°]>
## 1 <NA> (((112.8791 -7.574963, 112.8789 -7.574806, 112.8788 -7.574723, 112.878…jatim_kab <- jatim_kab[!is.na(jatim_kab$WADMKK), ]
unique(jatim_kab$WADMKK)## [1] "Bangkalan" "Banyuwangi" "Blitar" "Bojonegoro"
## [5] "Bondowoso" "Gresik" "Jember" "Jombang"
## [9] "Kediri" "Kota Batu" "Kota Blitar" "Kota Kediri"
## [13] "Kota Madiun" "Kota Malang" "Kota Mojokerto" "Kota Pasuruan"
## [17] "Kota Probolinggo" "Kota Surabaya" "Lamongan" "Lumajang"
## [21] "Madiun" "Magetan" "Malang" "Mojokerto"
## [25] "Nganjuk" "Ngawi" "Pacitan" "Pamekasan"
## [29] "Pasuruan" "Ponorogo" "Probolinggo" "Sampang"
## [33] "Sidoarjo" "Situbondo" "Sumenep" "Trenggalek"
## [37] "Tuban" "Tulungagung"table(jatim_kab$WADMKK)##
## Bangkalan Banyuwangi Blitar Bojonegoro
## 1 1 1 1
## Bondowoso Gresik Jember Jombang
## 1 1 1 1
## Kediri Kota Batu Kota Blitar Kota Kediri
## 1 1 1 1
## Kota Madiun Kota Malang Kota Mojokerto Kota Pasuruan
## 1 1 1 1
## Kota Probolinggo Kota Surabaya Lamongan Lumajang
## 1 1 1 1
## Madiun Magetan Malang Mojokerto
## 1 1 1 1
## Nganjuk Ngawi Pacitan Pamekasan
## 1 1 1 1
## Pasuruan Ponorogo Probolinggo Sampang
## 1 1 1 1
## Sidoarjo Situbondo Sumenep Trenggalek
## 1 1 1 1
## Tuban Tulungagung
## 1 1head(data)## # A tibble: 6 × 2
## `Kabupaten/Kota` produktivitas_padi
## <chr> <dbl>
## 1 Bangkalan 49
## 2 Banyuwangi 60.6
## 3 Blitar 66.1
## 4 Bojonegoro 53.2
## 5 Bondowoso 49.2
## 6 Gresik 63.0nrow(data)## [1] 38input data ke dalam peta
padi <- data[1:38,2]
padi.matrix = as.matrix(padi)
data.gab <- data.frame(jatim_kab, padi.matrix)
jatim_sf <- merge(jatim_kab,data.gab)AUTOKORELASI SPASIAL GLOBAL
# ubah ke objek spasial
jatim_sp <- as(jatim_sf, "Spatial")
spplot(jatim_sp, colnames(padi.matrix),scales=list(draw=F),colorkey=list(space="bottom"))
# memakai tmap
library(tmap)
library(cols4all)
library(sf)
padi_df <- data.frame(
Kabupaten_Kota = data$`Kabupaten/Kota`,
Jumlah = data$produktivitas_padi
)# Lalu gabungkan
jatim_kab <- jatim_kab %>%
left_join(padi_df, by = c("WADMKK" = "Kabupaten_Kota"))Plot peta
tm_shape(jatim_kab) +
tm_polygons(
col = "Jumlah", # variabel yang ingin dipetakan
palette = "YlGn", # palet warna dari RColorBrewer atau cols4all
title = "Produktivitas Padi" # judul legenda
) +
tm_text("WADMKK", size = 0.5, col = "black", fontface = "bold") +
tm_layout(
title = "Peta Sebaran Produktivitas Padi di Jawa Timur",
frame = FALSE
)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<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()`: 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 = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.
Uji autokorelasi spasial global
library(spdep)
nb <- poly2nb(jatim_sp, queen = TRUE)## Warning in poly2nb(jatim_sp, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style = "W")
moran_result <- moran.test(jatim_sp@data$produktivitas_padi, lw)# --- Ambil hasil penting
I_value <- moran_result$estimate["Moran I statistic"]
E_I <- moran_result$estimate["Expectation"]
p_value <- moran_result$p.value# --- Kesimpulan berdasarkan teori
cat("=== UJI AUTOKORELASI SPASIAL (MORAN’S I) ===\n")## === UJI AUTOKORELASI SPASIAL (MORAN’S I) ===cat("H0 : I = E[I] (tidak ada autokorelasi spasial)\n")## H0 : I = E[I] (tidak ada autokorelasi spasial)cat("H1 : I ≠E[I] (ada autokorelasi spasial)\n\n")## H1 : I ≠E[I] (ada autokorelasi spasial)cat("Nilai Moran’s I :", round(I_value, 4), "\n")## Nilai Moran’s I : 0.3523cat("Nilai E[I] :", round(E_I, 4), "\n")## Nilai E[I] : -0.027cat("p-value :", round(p_value, 4), "\n\n")## p-value : 0.0018# --- Interpretasi sesuai teori
if (p_value < 0.05) {
if (I_value > E_I) {
cat("Kesimpulan: TOLAK H0.\n")
cat("Terdapat autokorelasi spasial POSITIF (klaster),\n",
"karena wilayah yang berdekatan cenderung memiliki nilai serupa.\n")
} else {
cat("Kesimpulan: TOLAK H0.\n")
cat("Terdapat autokorelasi spasial NEGATIF (disperse),\n",
"karena wilayah yang berdekatan cenderung memiliki nilai yang berbeda.\n")
}
} else {
cat("Kesimpulan: GAGAL menolak H0.\n")
cat("Tidak terdapat autokorelasi spasial yang signifikan.\n",
"Nilai antar wilayah bersifat acak (random spatial pattern).\n")
}## Kesimpulan: TOLAK H0.
## Terdapat autokorelasi spasial POSITIF (klaster),
## karena wilayah yang berdekatan cenderung memiliki nilai serupa.auotkorelasi spasial lokal
# 0. load
library(sf)
library(tmap)
library(cols4all) # untuk palet
library(dplyr)# 1. Periksa struktur data dan nama kolom
print("NAMA KOLOM:")## [1] "NAMA KOLOM:"print(names(jatim_sf))## [1] "WADMKK" "produktivitas_padi" "geometry"print("Tipe objek (harus 'sf'):")## [1] "Tipe objek (harus 'sf'):"print(class(jatim_sf))## [1] "sf" "data.frame"print("Ringkasan kolom produktivitas (ganti nama jika perlu):")## [1] "Ringkasan kolom produktivitas (ganti nama jika perlu):"if("produktivitas_padi" %in% names(jatim_sf)) {
print(summary(jatim_sf$produktivitas_padi))
print(str(jatim_sf$produktivitas_padi))
} else {
message("Kolom 'produktivitas_padi' tidak ditemukan. Cek nama kolom di output 'names(jatim_sf)'.")
}## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 48.05 53.30 57.27 57.52 60.98 71.63
## num [1:38] 49 60.6 66.1 53.2 49.2 ...
## NULL# 2. Jika kolom masih bernama 'Jumlah', alias lain, kita buat kolom konsisten 'Jumlah'
# (ini tidak merusak; hanya menyalin jika memang ada)
if("Jumlah" %in% names(jatim_sf) & !"produktivitas_padi" %in% names(jatim_sf)) {
jatim_sf <- jatim_sf %>% mutate(produktivitas_padi = as.numeric(Jumlah))
} else if("produktivitas_padi" %in% names(jatim_sf)) {
jatim_sf <- jatim_sf %>% mutate(produktivitas_padi = as.numeric(produktivitas_padi))
} else {
stop("Tidak menemukan kolom 'Jumlah' atau 'produktivitas_padi'. Sesuaikan nama kolom dan jalankan ulang.")
}# 3. Cek apakah ada nilai non-numeric terselubung (NAs setelah coercion)
na_count <- sum(is.na(jatim_sf$produktivitas_padi))
message("Jumlah NA di kolom produktivitas_padi setelah coercion: ", na_count)## Jumlah NA di kolom produktivitas_padi setelah coercion: 0if(na_count > 0) {
message("Contoh nilai asli (variabel sebelum as.numeric) mungkin berisi karakter. Periksa dengan:")
print(head(jatim_sf %>% dplyr::select(starts_with("produktivitas")) , 20))
}# 4. Pastikan geom valid dan bukan GEOMETRYCOLLECTION aneh
if(!inherits(jatim_sf, "sf")) stop("jatim_sf bukan objek sf.")
if(any(!sf::st_is_valid(jatim_sf))) {
message("Beberapa geometri tidak valid — memperbaiki dengan st_make_valid()")
jatim_sf <- sf::st_make_valid(jatim_sf)
}# 5. Sekarang plot dengan sintaks tmap v4 (fill, fill.scale, fill.legend)
tm_shape(jatim_sf) +
tm_polygons(
fill = "produktivitas_padi", # harus persis nama kolom di data
fill.scale = tm_scale_continuous(
values = cols4all::c4a("brewer.yl_or_rd", n = 9) # palet dari cols4all
),
fill.legend = tm_legend(title = "Produktivitas Padi"),
col = "black", # warna batas
lwd = 0.5 # ketebalan batas
) +
tm_title("Peta Sebaran Produktivitas Padi di Jawa Timur (LISA-ready)") +
tm_layout(legend.outside = TRUE)## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
# Buat matriks bobot spasial (queen contiguity)
nb <- poly2nb(jatim_sf) # neighbors## Warning in poly2nb(jatim_sf): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw <- nb2listw(nb, style = "W") # weight list## Moran's Plot
mp <- moran.plot(as.vector(scale(jatim_sf$produktivitas_padi)), lw)
# Hitung Local Moran's I
loc_moran <- localmoran(jatim_sf$produktivitas_padi, lw, alternative = "greater")
head(loc_moran)## Ii E.Ii Var.Ii Z.Ii Pr(z > E(Ii))
## 1 2.9299557 -0.07124278 2.51435540 1.8926983 0.029199005
## 2 -0.5841307 -0.00926288 0.10978505 -1.7349876 0.958628503
## 3 1.7454531 -0.07230312 0.58411479 2.3784090 0.008693765
## 4 -0.2904635 -0.01830934 0.09802545 -0.8692518 0.807645294
## 5 0.3470295 -0.06793670 0.55142304 0.5588179 0.288143007
## 6 0.0229234 -0.03004145 0.25375179 0.1051437 0.458130899# Tambahkan hasil ke data sf
jatim_sf$lmI <- loc_moran[, "Ii"]
jatim_sf$lmZ <- loc_moran[, "Z.Ii"]
jatim_sf$lmp <- loc_moran[, "Pr(z > E(Ii))"]# Mode tampilan peta
tmap_mode("plot")## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`# Peta variabel asli (Sanitasi)
p1 <- tm_shape(jatim_sf) +
tm_polygons(
fill = "produktivitas_padi",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.yl_or_rd" # palet bawaan
),
fill.legend = tm_legend(title = "Sebaran Produktivitas Padi di Jawa Spasial")
) +
tm_layout(legend.outside = TRUE)
p1
# Peta Local Moran's I
p2 <- tm_shape(jatim_sf) +
tm_polygons(
fill = "lmI",
fill.scale = tm_scale_intervals(
style = "quantile",
values = "brewer.blues",
midpoint = NA
),
fill.legend = tm_legend(title = "Local Moran's I")
) +
tm_layout(legend.outside = TRUE)
p2
# Peta Z-score
p3 <- tm_shape(jatim_sf) +
tm_polygons(
fill = "lmZ",
fill.scale = tm_scale_intervals(
breaks = c(-Inf, 1.65, Inf),
labels = c("Tidak Signifikan", "Signifikan"),
values = c("grey80", "orange"),
midpoint = NA
),
fill.legend = tm_legend(title = "Z-score")
) +
tm_layout(legend.outside = TRUE)
p3
# Peta p-value
p4 <- tm_shape(jatim_sf) +
tm_polygons(
fill = "lmp",
fill.scale = tm_scale_intervals(
breaks = c(-Inf, 0.05, Inf),
labels = c("Signifikan", "Tidak signifikan"),
values = c("purple", "grey80")
),
fill.legend = tm_legend(title = "p-value")
) +
tm_layout(legend.outside = TRUE)
p4
# Gabungkan semua peta
tmap_arrange(p1, p2, p3, p4)## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
# cluster for quadrant
jatim_sf$quadrant <- NA
# high-high
jatim_sf[(mp$x >= 0 & mp$wx >= 0) & (jatim_sf$lmp <= 0.05), "quadrant"]<- 1
# low-low
jatim_sf[(mp$x <= 0 & mp$wx <= 0) & (jatim_sf$lmp <= 0.05), "quadrant"]<- 2
# high-low
jatim_sf[(mp$x >= 0 & mp$wx <= 0) & (jatim_sf$lmp <= 0.05), "quadrant"]<- 3
# low-high
jatim_sf[(mp$x <= 0 & mp$wx >= 0) & (jatim_sf$lmp <= 0.05), "quadrant"]<- 4
# non-significant
jatim_sf[(jatim_sf$lmp > 0.05), "quadrant"] <- 5
# Peta cluster quadrant
p_quadrant <- tm_shape(jatim_sf) +
tm_polygons(
fill = "quadrant",
fill.scale = tm_scale_intervals(
breaks = c(1, 2, 3, 4, 5, 6),
labels = c("High-High", "Low-Low", "High-Low", "Low-High", "Non-significant"),
values = c("red", "blue", "lightpink", "skyblue2", "white")
),
fill.legend = tm_legend(title = "")
) +
tm_borders(fill_alpha = 0.5) +
tm_text("WADMKK", size = 0.5, col = "black", fontface = "bold") +
tm_title("Clusters") +
tm_layout(frame = FALSE, legend.outside = TRUE)
p_quadrant