Praktikum 3 : Eksplorasi Data Areal

Author

Rosita Ria Rusesta - Statistika dan Sains Data

Review Data Areal/Lattice

Data areal (lattice data) merupakan bentuk data spasial yang bersifat diskrit karena merepresentasikan informasi dalam batas wilayah tertentu. Data ini muncul ketika suatu fenomena atau jumlah kejadian digabungkan ke dalam area administratif, misalnya jumlah kasus penyakit dalam epidemiologi spasial atau jumlah kelahiran pada suatu kabupaten. Dengan pendekatan ini, peneliti dapat mengidentifikasi pola geografis, memahami distribusi fenomena, dan menelusuri faktor risiko yang mungkin berhubungan dengan kondisi lingkungan maupun sosial.

Dalam praktik analisis, data areal umumnya tersedia dalam bentuk objek vektor dengan format shapefile. Format ini menjadi standar paling umum untuk menyimpan data spasial karena mampu memuat informasi geometri (seperti poligon batas wilayah) sekaligus atribut nonspasial yang menyertainya. Di R, data spasial tersebut dapat dikelola menggunakan paket sf yang dirancang khusus untuk memudahkan manipulasi dan analisis data berbasis spasial. Sementara itu, data nonspasial biasanya disimpan dalam format umum seperti file CSV atau .xlsx. Dengan kombinasi kedua jenis data ini, analisis spasial dapat dilakukan secara lebih komprehensif, baik dari sisi geografis maupun dari sisi atribut yang melekat pada setiap area.

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(sf)
library(sp)
library(mapview)
library(spData)
library(ggplot2)
library(spatstat)
library(tidyverse)
library(spData)
library(readxl)
list.files(system.file("shapes", package = "spData"))
 [1] "auckland.gpkg"          "boston_tracts.gpkg"     "columbus.gpkg"         
 [4] "cycle_hire_osm.geojson" "cycle_hire.geojson"     "eire.gpkg"             
 [7] "NY8_bna_utm18.gpkg"     "NY8_utm18.gpkg"         "sids.gpkg"             
[10] "wheat.gpkg"             "world.gpkg"            
d2 <- st_read(system.file("shapes/columbus.gpkg",
                         package = "spData"), quiet = TRUE)
st_crs(d2) <- NA
names(d2)
 [1] "AREA"       "PERIMETER"  "COLUMBUS_"  "COLUMBUS_I" "POLYID"    
 [6] "NEIG"       "HOVAL"      "INC"        "CRIME"      "OPEN"      
[11] "PLUMB"      "DISCBD"     "X"          "Y"          "NSA"       
[16] "NSB"        "EW"         "CP"         "THOUS"      "NEIGNO"    
[21] "geom"      
print(d2)
Simple feature collection with 49 features and 20 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 5.874907 ymin: 10.78863 xmax: 11.28742 ymax: 14.74245
CRS:           NA
First 10 features:
       AREA PERIMETER COLUMBUS_ COLUMBUS_I POLYID NEIG  HOVAL    INC     CRIME
1  0.309441  2.440629         2          5      1    5 80.467 19.531 15.725980
2  0.259329  2.236939         3          1      2    1 44.567 21.232 18.801754
3  0.192468  2.187547         4          6      3    6 26.350 15.956 30.626781
4  0.083841  1.427635         5          2      4    2 33.200  4.477 32.387760
5  0.488888  2.997133         6          7      5    7 23.225 11.252 50.731510
6  0.283079  2.335634         7          8      6    8 28.750 16.029 26.066658
7  0.257084  2.554577         8          4      7    4 75.000  8.438  0.178269
8  0.204954  2.139524         9          3      8    3 37.125 11.337 38.425858
9  0.500755  3.169707        10         18      9   18 52.600 17.586 30.515917
10 0.246689  2.087235        11         10     10   10 96.400 13.598 34.000835
       OPEN    PLUMB DISCBD     X     Y NSA NSB EW CP THOUS NEIGNO
1  2.850747 0.217155   5.03 38.80 44.07   1   1  1  0  1000   1005
2  5.296720 0.320581   4.27 35.62 42.38   1   1  0  0  1000   1001
3  4.534649 0.374404   3.89 39.82 41.18   1   1  1  0  1000   1006
4  0.394427 1.186944   3.70 36.50 40.52   1   1  0  0  1000   1002
5  0.405664 0.624596   2.83 40.01 38.00   1   1  1  0  1000   1007
6  0.563075 0.254130   3.78 43.75 39.28   1   1  1  0  1000   1008
7  0.000000 2.402402   2.74 33.36 38.41   1   1  0  0  1000   1004
8  3.483478 2.739726   2.89 36.71 38.71   1   1  0  0  1000   1003
9  0.527488 0.890736   3.17 43.44 35.92   1   1  1  0  1000   1018
10 1.548348 0.557724   4.33 47.61 36.42   1   1  1  0  1000   1010
                             geom
1  POLYGON ((8.624129 14.23698...
2  POLYGON ((8.25279 14.23694,...
3  POLYGON ((8.653305 14.00809...
4  POLYGON ((8.459499 13.82035...
5  POLYGON ((8.685274 13.63952...
6  POLYGON ((9.401384 13.5504,...
7  POLYGON ((8.037741 13.60752...
8  POLYGON ((8.247527 13.58651...
9  POLYGON ((9.333297 13.27242...
10 POLYGON ((10.08251 13.03377...
ggplot(d2) + geom_sf(aes(fill = CRIME)) + theme_minimal()

Membuat objek sf

Objek sf (simple feature) di R adalah struktur data yang digunakan untuk menyimpan data spasial vektor. Objek ini menggabungkan data atribut (seperti data frame) dengan geometri spasial.

  1. Data Atribut:

    Mirip dengan data frame biasa, bagian ini berisi kolom-kolom yang menyimpan informasi non-spasial (atribut) tentang setiap fitur. Misalnya, nama, ID, populasi, dll.

  2. Kolom Geometri:

    Kolom khusus yang menyimpan informasi geometri dari setiap fitur. Kolom ini biasanya dinamai geometry atau geom, tetapi bisa diberi nama lain. Geometri ini bisa berupa titik, garis, poligon, multipoint, atau multipoligon.

  3. Sistem Referensi Koordinat (CRS):

    Informasi tentang sistem referensi koordinat yang digunakan untuk geometri. CRS menentukan bagaimana koordinat dipetakan ke permukaan bumi.

  4. Atribut Metadata:

    Informasi tambahan tentang objek sf, seperti nama kolom geometri, CRS, dan tipe geometri.

Kita dapat menggunakan st_sf() fungsi tersebut untuk membuat objek sf dengan menyediakan dua elemen, yaitu data frame dengan atribut setiap fitur dan kolom daftar geometri fitur sederhana (simple feature geometri list-column/sfc) yang berisi geometri fitur sederhana (simple feature geomteries/sfg).

# Membuat data frame dengan atribut
df <- data.frame(id = 1:3,
                 name = c("A", "B", "C"))

# Membuat geometri (titik)
p1 <- st_point(c(7.35, 52.42))
p2 <- st_point(c(7.22, 52.18))
p3 <- st_point(c(7.44, 52.19))
geom <- st_sfc(list(p1, p2, p3), crs = 'EPSG:4326')

# Menggabungkan data frame dan geometri menjadi objek sf
sf_obj <- st_sf(df, geometry = geom)

# Menampilkan objek sf
print(sf_obj)
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.22 ymin: 52.18 xmax: 7.44 ymax: 52.42
Geodetic CRS:  WGS 84
  id name           geometry
1  1    A POINT (7.35 52.42)
2  2    B POINT (7.22 52.18)
3  3    C POINT (7.44 52.19)

  • POINT Satu titik. Contoh: lokasi sensor atau titik GPS.

  • MULTIPOINT Sekumpulan titik dalam satu fitur. Contoh: sebaran pohon dalam satu plot.

  • LINESTRING Garis tunggal yang terdiri atas segmen-segmen. Contoh: satu ruas jalan atau sungai.

  • MULTILINESTRING Beberapa garis dalam satu fitur. Contoh: jaringan jalan yang dimiliki satu entitas.

  • POLYGON Satu poligon dengan cangkang luar (boleh ada lubang/holes). Contoh: batas satu kelurahan.

  • MULTIPOLYGON Sekumpulan poligon dalam satu fitur. Contoh: kabupaten yang terdiri atas beberapa pulau terpisah.

  • GEOMETRYCOLLECTION Kumpulan berbagai tipe geometri (titik, garis, poligon) dalam satu fitur.

Implementasi Data Area : Data Jawa Barat

Data yang digunakan dalam praktik ini merupakan data sekunder berupa data kependudukan yang diperoleh dari Badan Pusat Statistik (BPS) dan peta shapefile berisi peta Provinsi Jawa Barat. Berdasarkan data tersebut, tahap awal yaitu proses menggabungkan data bentuk geometri wilayah (peta shp Jawa Barat) dengan data statistik atribut dalam format xlsx. Dataset mencakup informasi spasial dan tabular tingkat kabupaten/kota di Provinsi Jawa Barat untuk tahun 2015 dan 2016. Praktik ini melibatkan variabel utama, yaitu %Kemiskinan tahun 2015. Data dapat didownload pada link berikut https://bit.ly/data_praktikum3

setwd("/Users/rositariarusesta/Documents/Regresi Spasial")
shp_jabar <- st_read("petaJabar2/Jabar2.shp")
Reading layer `Jabar2' from data source 
  `/Users/rositariarusesta/Documents/Regresi Spasial/PetaJabar2/Jabar2.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 26 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
CRS:           NA
data_jabar <- read_excel("data_jabar.xlsx")
data_jabar$KABKOTNO <- as.character(data_jabar$KABKOTNO)

# join
data_spasial <- shp_jabar %>%
  left_join(data_jabar, by = "KABKOTNO") %>%
  na.omit()
#Cek semua nama kolom
names(data_spasial)
 [1] "PROVNO.x"   "KABKOTNO"   "KODE2010.x" "PROVINSI.x" "KABKOT.x"  
 [6] "SUMBER"     "IDSP2010.x" "PROVNO.y"   "KODE2010.y" "PROVINSI.y"
[11] "KABKOT.y"   "IDSP2010.y" "Long"       "Lat"        "p.miskin15"
[16] "p.miskin16" "j.miskin15" "j.miskin16" "AHH2015"    "AHH2016"   
[21] "EYS2015"    "EYS2016"    "MYS2015"    "MYS2016"    "EXP2015"   
[26] "EXP2016"    "APM.SD15"   "APM.SMP15"  "APM.SMA15"  "APM.PT15"  
[31] "APK.SD15"   "APK.SMP15"  "APK.SMA15"  "APK.PT15"   "APS.USIA15"
[36] "APS.USIA2"  "APS.USIA3"  "APS.USIA4"  "geometry"  
#Memilih data yang mengandung variabel saja
data_filter<-data_spasial[,15:37]
#Visualisasi data Jabar
plot(data_filter)

library(ggplot2)

plot_kemiskinan <- ggplot(data = data_spasial) +  
  geom_sf(aes(fill = p.miskin15), color = "white", linewidth = 0.2) +  
  # Menambahkan label nama kabupaten/kota:
  geom_sf_text(aes(label = KABKOT.x), size = 2, color = "black", check_overlap = TRUE) +
  scale_fill_viridis_c(option = "magma", name = "Persentase (%)") + 
  theme_minimal() +  
  labs(
    title = "%Kemiskinan Jawa Barat",       
    subtitle = "Tahun 2015",       
    x = "Longitude", 
    y = "Latitude"
  ) +  
  theme(legend.position = "bottom")

plot_kemiskinan

Explorasi Data Area

Import Dara Area

Kita dapat membaca objek sf dengan fungsi st_read(). Misalnya, di sini kita membaca shapefile ncyang berisi kabupaten-kabupaten di North Carolina, AS, beserta nama kabupaten, jumlah kelahiran, dan jumlah kematian bayi mendadak pada tahun 1974 dan 1979.

imp <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(imp, quiet = TRUE)
class(nc)
[1] "sf"         "data.frame"

Selanjutnya kita dapat mengidentifikasi CRS yang digunakan pada data tersebut dengan fungsi st_crs().

st_crs(nc)
Coordinate Reference System:
  User input: NAD27 
  wkt:
GEOGCRS["NAD27",
    DATUM["North American Datum 1927",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4267]]
plot(nc)

Kita dapat membuat subset set fitur dengan menggunakan notasi tanda kurung siku dan menggunakan argumen drop untuk menghilangkan geometri.

plot(nc[9])

Menghapus Beberapa Poligon

map <- nc[-which(nc$FIPS %in% c("37125", "37051")), ]
ggplot(nc) + geom_sf(aes(fill = SID79))

Menghitung Jumlah Titik dalam Poligon

Kita dapat menggunakan st_intersects() untuk menghitung jumlah titik dalam poligon suatu objek sfdan st_sample() menghasilkan titik sampel acak. Objek yang dikembalikan adalah daftar dengan id fitur yang berpotongan di setiap poligon.

# Peta (objek sf)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

# Titik di atas peta (kolom daftar geometri fitur sederhana sfc)
points <- st_sample(map, size = 100)

# Peta titik di dalam poligon
ggplot() + geom_sf(data = map) + geom_sf(data = points)

# Interseksi (argumen pertama peta, kemudian titik)
inter <- st_intersects(map, points)

# Menambahkan jumlah titik ke setiap poligon
map$count <- lengths(inter)

# Peta jumlah titik di dalam poligon
ggplot(map) + geom_sf(aes(fill = count))

Peta Choropleth

Peta Choropleth merupakan peta visualisasi utama untuk data areal dengan memberikan gradasi warna pada setiap wilayah berdasarkan nilai peubahnya. Penggunaan nilai relatif sangat disarakankan (seperti rate, proporsi(%), atau kepadatan) karena nilai mutlak (jumlah mentah) sering kali menyesatkan jika tidak memperhitungkan ukuran populasi atau skala wilayah.

Skema Kelas dan Warna

Proses pengelompokkan data kontinu ke dalam kategori (class interval) untuk menentukan pewarnaan peta. Metode Klasifikasi : (1) Interval sama (Equal Interval), (2) Kuantil (Quantile), (3) Simpangan Baku (Standard Deviation), dan (4) Natural Breaks (Jenks)

Implementasi dalam R

Instalasi dan pemanggilan package

# Instalasi (cukup dilakukan sekali)
# install.packages(c("sf", "classInt", "tmap"))

# Memanggil package
library(sf)
library(classInt)
library(tmap)

Kita gunakan data wilayah bawaan sf, yaitu North Carolina pada contoh sebelumnya. Misalnya kita ingin memetakan jumlah penduduk per county. Variabel yang dapat digunakan misalnya:

#jumlah kelahiran pada tahun 1974.
nc$BIR74
  [1]  1091   487  3188   508  1421  1452   286   420   968  1612  1035  4449
 [13]  1671  1556  2180  3608  1638  3146  1323   484   751   781  1269  1399
 [25] 11858 16184  4672  1324  3164  7970  4021   671  3657  3609   770  1549
 [37] 14484   765  4139  1207  1333  5509  3573   990   248  1946  4456  1646
 [49]  3702  4606  5094  5754  7515  3999  2110   521  2692   675   870  2252
 [61]  2992  6638  3776  4866  2216  1143  2648 21588  4099  1258  2356  2574
 [73]   415  3589  1173  9014   533   797  3025   542  1027 20366   578  3915
 [85]  1570  1494   338  2483  2756   284  5868  2255 11158  7889  2414  1782
 [97]  1228  3350  5526  2181

Melihat Distribusi

summary(nc$BIR74)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    248    1077    2180    3300    3936   21588

Tentukan jumlah kelas

Gunakan aturan Sturges:

k=1+3,322 log⁡10(n)

dengan:

  • k = jumlah kelas

  • n = jumlah observasi/wilayah

n <- length(nc$BIR74)

k <- 1 + 3.322 * log10(n)
k <- round(k)
k
[1] 8

Metode Klasifikasi

1) Equal Interval

Equal Interval (Interval Sama) digunakan untuk membagi rentang nilai atribut menjadi subrentang dengan ukuran yang sama. Equal Interval paling sesuai digunakan untuk data dengan rentang nilai yang sudah umum atau mudah dipahami, seperti persentase dan suhu. Metode ini menekankan besarnya suatu nilai atribut dibandingkan dengan nilai lainnya.

Contoh:

Jika data berkisar dari 0–300 dan ingin dibuat 3 kelas:

Lebar interval=300 / 3 =100

sehingga:

Kelas Rentang
1 0–100
2 100–200
3 200–300 
# Membuat klasifikasi Equal Interval
ei <- classIntervals(
  nc$BIR74,
  n = k,
  style = "equal"
)

ei
style: equal
  one of 14,887,031,544 possible partitions of this variable into 8 classes
   [248,2915.5)   [2915.5,5583)   [5583,8250.5)  [8250.5,10918) [10918,13585.5) 
             61              26               6               1               2 
[13585.5,16253) [16253,18920.5) [18920.5,21588] 
              2               0               2 
ei$brks
[1]   248.0  2915.5  5583.0  8250.5 10918.0 13585.5 16253.0 18920.5 21588.0

Membuat Variabel Kelas

nc$kelas_ei <- cut(
  nc$BIR74,
  breaks = ei$brks,
  include.lowest = TRUE
)

head(nc[, c("BIR74", "kelas_ei")])
Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
  BIR74            kelas_ei                       geometry
1  1091      [248,2.92e+03] MULTIPOLYGON (((-81.47276 3...
2   487      [248,2.92e+03] MULTIPOLYGON (((-81.23989 3...
3  3188 (2.92e+03,5.58e+03] MULTIPOLYGON (((-80.45634 3...
4   508      [248,2.92e+03] MULTIPOLYGON (((-76.00897 3...
5  1421      [248,2.92e+03] MULTIPOLYGON (((-77.21767 3...
6  1452      [248,2.92e+03] MULTIPOLYGON (((-76.74506 3...

Membuat Peta

tm_shape(nc) +
  tm_polygons(
    "kelas_ei",
    title = "Equal Interval"
  )

2) Quantile

Pada klasifikasi Quantile (Kuantil), setiap kelas berisi jumlah objek/wilayah (features) yang sama. Metode Quantile memberikan jumlah nilai data yang sama pada setiap kelas. Dengan demikian, tidak ada kelas yang kosong atau kelas yang memiliki jumlah nilai terlalu sedikit maupun terlalu banyak.

Namun, karena objek/wilayah dikelompokkan dalam jumlah yang sama pada setiap kelas, peta yang dihasilkan dapat terkadang menyesatkan. Objek yang memiliki nilai sangat mirip dapat ditempatkan pada kelas yang berbeda, sementara objek yang memiliki nilai sangat berbeda dapat dimasukkan ke dalam kelas yang sama. Distorsi tersebut dapat dikurangi dengan meningkatkan jumlah kelas.

Misalnya terdapat 100 kecamatan dan dibuat 5 kelas:

Kelas Jumlah kecamatan
Sangat rendah 20
Rendah 20
Sedang 20
Tinggi 20
Sangat tinggi 20

Jadi, yang dibuat sama adalah jumlah observasinya, bukan rentang nilainya.

q <- classIntervals(
  nc$BIR74,
  n = k,
  style = "quantile"
)

q$brks
[1]   248.000   612.875  1077.000  1457.250  2180.500  3020.875  3936.000
[8]  5668.500 21588.000
nc$kelas_quantile <- cut(
  nc$BIR74,
  breaks = q$brks,
  include.lowest = TRUE
)
tm_shape(nc) +
  tm_polygons(
    "kelas_quantile",
    title = "Quantile"
  )

3) Standard Deviation

Standard Deviation adalah metode klasifikasi yang mengelompokkan data berdasarkan seberapa jauh nilai suatu wilayah dari nilai rata-rata. Batas kelas ditentukan berdasarkan kelipatan tertentu dari simpangan baku, misalnya ±1 SD, ±½ SD, ±⅓ SD, atau ±¼ SD.

sd_class <- classIntervals(
  nc$BIR74,
  n = k,
  style = "sd"
)

sd_class$brks
 [1]  -548.5451  1375.5374  3299.6200  5223.7026  7147.7851  9071.8677
 [7] 10995.9503 12920.0328 14844.1154 16768.1979 18692.2805 20616.3631
[13] 22540.4456
nc$kelas_sd <- cut(
  nc$BIR74,
  breaks = sd_class$brks,
  include.lowest = TRUE
)
tm_shape(nc) +
  tm_polygons(
    "kelas_sd",
    title = "Standard Deviation"
  )

4) Natural Breaks (Jenks)

Natural Breaks (Jenks) adalah metode klasifikasi yang mencari kelompok alami dalam data. Nilai yang mirip akan dikelompokkan dalam kelas yang sama, sedangkan batas kelas ditempatkan pada bagian yang memiliki perbedaan nilai yang relatif besar. Tujuannya adalah meminimalkan perbedaan di dalam kelas dan memaksimalkan perbedaan antar kelas.

Catatan : Natural Breaks kurang tepat untuk membandingkan peta antarperiode jika batas kelas dihitung ulang untuk setiap peta, karena kelasnya mengikuti distribusi masing-masing data. Jadi, misalnya untuk membandingkan peta 2015, 2020, dan 2025, penggunaan Jenks perlu diperhatikan agar perubahan warna tidak hanya mencerminkan perubahan batas kelas

jenks <- classIntervals(
  nc$BIR74,
  n = k,
  style = "jenks"
)

jenks
style: jenks
  one of 14,887,031,544 possible partitions of this variable into 8 classes
   [248,1421]   (1421,2756]   (2756,4139]   (4139,5868]   (5868,9014] 
           37            24            18            10             5 
 (9014,11858] (11858,16184] (16184,21588] 
            2             2             2 
jenks <- classIntervals(
  nc$BIR74,
  n = 5,
  style = "jenks"
)

jenks
style: jenks
  one of 3,764,376 possible partitions of this variable into 5 classes
   [248,2356]   (2356,5094]   (5094,9014]  (9014,16184] (16184,21588] 
           55            30             9             4             2 
jenks$brks
[1]   248  2356  5094  9014 16184 21588
nc$kelas_jenks <- cut(
  nc$BIR74,
  breaks = jenks$brks,
  include.lowest = TRUE
)
tm_shape(nc) +
  tm_polygons(
    "kelas_jenks",
    title = "Natural Breaks (Jenks)"
  )

Gabungan Semua Pendekatan Klasifikasi

library(sf)
library(tmap)

# 2. Definisikan gaya klaisfikasi
style_list <- c("equal", "quantile", "jenks", "sd")

# 3. Buat peta dalam loop
maps <- lapply(style_list, function(s) {
  tm_shape(nc) +
    tm_polygons(
      fill = "BIR74",  # Menyebutkan nama kolom variabel
      fill.scale = tm_scale_intervals(
        style = s,
        values = "YlOrRd",
        n = k
      ),
      fill.legend = tm_legend(title = s)
    ) +
    tm_title(paste("style =", s), size = 0.9)
})

# 4. Tampilkan 4 peta sekaligus menggunakan do.call (Aman dari out of bounds)
do.call(tmap_arrange, c(maps, list(ncol = 2)))

Peta Cartogram

Kartogram adalah peta tematik yang mengubah ukuran atau bentuk wilayah secara proporsional terhadap nilai variabel tertentu. Pada peta biasa, ukuran wilayah mengikuti kondisi geografis sebenarnya. Pada kartogram, ukuran wilayah didistorsi agar lebih mencerminkan besarnya variabel yang dianalisis.

Misalnya kita ingin memetakan jumlah penduduk:

Kabupaten dengan jumlah penduduk besar → ukuran wilayah diperbesar
Kabupaten dengan jumlah penduduk kecil → ukuran wilayah diperkecil.

Sehingga kartogram menjawab pertanyaan:

“Di mana nilai variabel paling besar atau paling kecil?”

Bayangkan kita membuat peta Indonesia berdasarkan jumlah penduduk.

Pada peta geografis biasa:

  • Kalimantan memiliki wilayah yang sangat luas.

  • Jawa memiliki wilayah yang jauh lebih kecil.

Tetapi jumlah penduduk Jawa jauh lebih besar.

Jika menggunakan luas geografis, kita bisa mendapatkan kesan bahwa Kalimantan “lebih besar” secara demografis. Kartogram mencoba mengubah perspektif tersebut.

library(cartogram)
ggplot(nc) +
  geom_sf(aes(fill = BIR74)) +
  theme_minimal() +
  labs(
    title = "Jumlah Kelahiran Tahun 1974",
    fill = "Jumlah"
  )

Pada peta ini, ukuran wilayah tidak dipengaruhi oleh jumlah kelahiran. Ukuran mengikuti kondisi geografis sebenarnya.

Karena algoritma cartogram melakukan perhitungan geometris seperti jarak, centroid, dan perubahan ukuran/bentuk wilayah sehinggamembutuhkan sistem koordinat terproyeksi. Kalau masih menggunakan longitude/latitude (derajat), jarak antar-koordinat tidak memiliki satuan panjang yang seragam seperti meter.

#Cek CRS
st_crs(nc)
Coordinate Reference System:
  User input: NAD27 
  wkt:
GEOGCRS["NAD27",
    DATUM["North American Datum 1927",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4267]]
#projected_crs untuk data nc yang awalnya masih berbentung geometri(long, lat)
nc_proj <- st_transform(nc, 32119)
st_crs(nc_proj)
Coordinate Reference System:
  User input: EPSG:32119 
  wkt:
PROJCRS["NAD83 / North Carolina",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["SPCS83 North Carolina zone (meter)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",33.75,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-79,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",36.1666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",34.3333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",609601.22,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["United States (USA) - North Carolina - counties of Alamance; Alexander; Alleghany; Anson; Ashe; Avery; Beaufort; Bertie; Bladen; Brunswick; Buncombe; Burke; Cabarrus; Caldwell; Camden; Carteret; Caswell; Catawba; Chatham; Cherokee; Chowan; Clay; Cleveland; Columbus; Craven; Cumberland; Currituck; Dare; Davidson; Davie; Duplin; Durham; Edgecombe; Forsyth; Franklin; Gaston; Gates; Graham; Granville; Greene; Guilford; Halifax; Harnett; Haywood; Henderson; Hertford; Hoke; Hyde; Iredell; Jackson; Johnston; Jones; Lee; Lenoir; Lincoln; Macon; Madison; Martin; McDowell; Mecklenburg; Mitchell; Montgomery; Moore; Nash; New Hanover; Northampton; Onslow; Orange; Pamlico; Pasquotank; Pender; Perquimans; Person; Pitt; Polk; Randolph; Richmond; Robeson; Rockingham; Rowan; Rutherford; Sampson; Scotland; Stanly; Stokes; Surry; Swain; Transylvania; Tyrrell; Union; Vance; Wake; Warren; Washington; Watauga; Wayne; Wilkes; Wilson; Yadkin; Yancey."],
        BBOX[33.83,-84.33,36.59,-75.38]],
    ID["EPSG",32119]]
nc_cartogram <- cartogram_cont(
  nc_proj,
  weight = "BIR74"
)

Implementasi Peta Cartogram di R

library(cartogram)

nc_cartogram <- cartogram_cont(
  nc_proj,
  weight = "BIR74"
)
ggplot(nc_cartogram) +
  geom_sf(aes(fill = BIR74)) +
  theme_void() +
  labs(
    title = "Contiguous Cartogram",
    subtitle = "Ukuran wilayah berdasarkan jumlah kelahiran 1974",
    fill = "Jumlah"
  )

Kartogram digunakan ketika fokus utama adalah menunjukkan besarnya suatu variabel, bukan mempertahankan bentuk geografis wilayah.

Dan langsung beri peringatan:

⚠️ Keterbatasan: semakin besar distorsi geometris, semakin sulit pembaca mengenali lokasi dan hubungan spasial antarwilayah.

Menambahkan label pada peta

#Label pada top10
top10 <- nc_cartogram |>
  dplyr::arrange(desc(BIR74)) |>
  dplyr::slice_head(n = 10)

ggplot(nc_cartogram) +
  geom_sf(aes(fill = BIR74)) +
  geom_sf_text(
    data = top10,
    aes(label = NAME),
    size = 3
  ) +
  theme_void() +
  labs(
    title = "Contiguous Cartogram",
    subtitle = "10 wilayah dengan jumlah kelahiran tertinggi",
    fill = "Jumlah"
  )

Isu Data Areal

1. MAUP (Modifiable Areal Unit Problem) MAUP terjadi ketika hasil statistik atau analisis spasial berubah-ubah akibat perubahan batas administrasi/wilayah pengamatan. MAUP terdiri dari dua efek utama:

  • Efek Skala (Scale Effect): Hasil analisis berubah jika unit wilayah digabungkan (aggregated) ke tingkat yang lebih besar

  • Efek Zonasi (Zoning Effect): Hasil analisis berubah jika bentuk atau susunan batas wilayah diubah, meskipun agregasi/jumlah wilayahnya sama.

2. MIDP (Modifiable Individual Data Problem) Mirip dengan MAUP, namun berfokus pada titik data level individu (point data) yang dikelompokkan ke dalam area statistik buatan.

3. Bias Ekologis (Ecological Fallacy / Bias) Kesalahan logika saat mengambil kesimpulan mengenai individu berdasarkan hasil analisis dari data agregat (level kelompok/area). Contoh: Kabupaten A memiliki rerata pendapatan tinggi, sehingga dapat disimpulkan seluruh warga di Kabupaten A pasti kaya (ini adalah bias ekologis).

Ilustrasi MAUP

a) Efek Skala

Contoh kasus yaitu mengagregasikan data spasial dari level Desa (skala kecil) ke level Kecamatan (skala besar) menggunakan paket sfdan dplyr. Hal ini untuk mengilustrasikan bagaimana nilai rata-rata/korelasi berubah akibat Efek Skala (MAUP).

library(patchwork) 

#Membuat Data Spasial Buatan (Grid 4x4 = 16 Desa)
set.seed(123)
# Buat bounding box dasar
bbox_dasar <- st_bbox(c(xmin = 0, ymin = 0, xmax = 4, ymax = 4), crs = 4326)
# Buat grid berdasarkan bbox dasar
grid <- st_make_grid(
  x = bbox_dasar,
  cellsize = c(1, 1),
  n = c(4, 4)
) %>% st_sf()

# Menambahkan atribut data simulasi level Desa
data_desa <- grid %>%
  mutate(
    desa_id = paste0("Desa_", 1:n()),
    # Mengelompokkan 16 Desa menjadi 4 Kecamatan (2x2 grid)
    kecamatan_id = paste0("Kecamatan_", c(1, 1, 2, 2, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 4, 4)),
    pendapatan = round(runif(n(), min = 20, max = 100), 2),
    kemiskinan = round(100 - pendapatan + rnorm(n(), mean = 0, sd = 5), 2)
  )
# Agregasi Data ke Level Kecamatan (Efek Skala)
data_kecamatan <- data_desa %>%
  group_by(kecamatan_id) %>%
  summarise(
    pendapatan_mean = mean(pendapatan),
    kemiskinan_mean = mean(kemiskinan),
    .groups = "drop"
  )
# Bandingkan Statistik & Korelasi (Level Desa vs. Kecamatan)
#==================================================
#EFEK SKALA MAUP: HASIL ANALISIS
#==================================================
kor_desa <- cor(data_desa$pendapatan, data_desa$kemiskinan)
kor_desa
[1] -0.9799158
kor_kec <- cor(data_kecamatan$pendapatan_mean, data_kecamatan$kemiskinan_mean)
kor_kec
[1] -0.847398
# Visualisasi Perubahan Skala Wilayah
plot_desa <- ggplot(data_desa) +
  geom_sf(aes(fill = pendapatan), color = "white") +
  geom_sf_text(aes(label = desa_id), size = 2.5) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Skala Asli: Desa (16 Unit)", fill = "Pendapatan") +
  theme_minimal()

plot_kecamatan <- ggplot(data_kecamatan) +
  geom_sf(aes(fill = pendapatan_mean), color = "white") +
  geom_sf_text(aes(label = kecamatan_id), size = 3, fontface = "bold") +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Skala Agregat: Kecamatan (4 Unit)", fill = "Rerata") +
  theme_minimal()

plot_desa + plot_kecamatan

b) Efek Zonasi

# 2. Membuat Data Spasial Buatan (Grid 4x4 = 16 Desa)
set.seed(123)

bbox_dasar <- st_bbox(c(xmin = 0, ymin = 0, xmax = 4, ymax = 4), crs = 4326)

grid <- st_make_grid(
  x = bbox_dasar,
  cellsize = c(1, 1),
  n = c(4, 4)
) %>% st_sf()

# Menambahkan atribut data simulasi & 2 Skema Zonasi Berbeda (Sama-sama 4 Unit)
data_desa <- grid %>%
  mutate(
    desa_id = paste0("D", 1:n()),
    
    # ZONASI A: Skema Grid 2x2 (4 Kecamatan)
    zona_grid = paste0("Zona_A_", c(1,1,2,2, 1,1,2,2, 3,3,4,4, 3,3,4,4)),
    
    # ZONASI B: Skema Strip Horizontal 4x1 (4 Kecamatan)
    zona_horiz = paste0("Zona_B_", c(1,1,1,1, 2,2,2,2, 3,3,3,3, 4,4,4,4)),
    
    # Data Atribut Simulasi
    pendapatan = round(runif(n(), min = 20, max = 100), 2),
    kemiskinan = round(100 - pendapatan + rnorm(n(), mean = 0, sd = 5), 2)
  )
head(data_desa)
Simple feature collection with 6 features and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 4 ymax: 2
Geodetic CRS:  WGS 84
                        geometry desa_id zona_grid zona_horiz pendapatan
1 POLYGON ((0 0, 1 0, 1 1, 0 ...      D1  Zona_A_1   Zona_B_1      43.01
2 POLYGON ((1 0, 2 0, 2 1, 1 ...      D2  Zona_A_1   Zona_B_1      83.06
3 POLYGON ((2 0, 3 0, 3 1, 2 ...      D3  Zona_A_2   Zona_B_1      52.72
4 POLYGON ((3 0, 4 0, 4 1, 3 ...      D4  Zona_A_2   Zona_B_1      90.64
5 POLYGON ((0 1, 1 1, 1 2, 0 ...      D5  Zona_A_1   Zona_B_2      95.24
6 POLYGON ((1 1, 2 1, 2 2, 1 ...      D6  Zona_A_1   Zona_B_2      23.64
  kemiskinan
1      53.56
2      14.71
3      53.40
4      11.16
5       6.76
6      76.91
# =============================================================================
# 3. AGREGASI DATA BERDASARKAN DUA SKEMA ZONASI
# =============================================================================
# Agregasi Zonasi A (Grid 2x2)
data_zona_A <- data_desa %>%
  group_by(zona_grid) %>%
  summarise(
    pendapatan = mean(pendapatan),
    kemiskinan = mean(kemiskinan),
    .groups = "drop"
  )

# Agregasi Zonasi B (Strip Horizontal 4x1)
data_zona_B <- data_desa %>%
  group_by(zona_horiz) %>%
  summarise(
    pendapatan = mean(pendapatan),
    kemiskinan = mean(kemiskinan),
    .groups = "drop"
  )
# =============================================================================
# 4. HASIL PERBANDINGAN STATISTIK (EFEK ZONASI)
# ==============================================================================
r_desa   <- cor(data_desa$pendapatan, data_desa$kemiskinan)
r_desa
[1] -0.9799158
r_zona_A <- cor(data_zona_A$pendapatan, data_zona_A$kemiskinan)
r_zona_A
[1] -0.847398
r_zona_B <- cor(data_zona_B$pendapatan, data_zona_B$kemiskinan)
r_zona_B
[1] 0.2049043
# =============================================================================
# 5. VISUALISASI PERBANDINGAN ZONASI
# =============================================================================

plot_A <- ggplot(data_zona_A) +
  geom_sf(aes(fill = pendapatan), color = "white") +
  geom_sf_text(aes(label = zona_grid), size = 3, fontface = "bold") +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Zonasi A: Grid (2x2)", subtitle = "4 Unit Wilayah", fill = "Rerata") +
  theme_minimal()

plot_B <- ggplot(data_zona_B) +
  geom_sf(aes(fill = pendapatan), color = "white") +
  geom_sf_text(aes(label = zona_horiz), size = 3, fontface = "bold") +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Zonasi B: Horizontal (4x1)", subtitle = "4 Unit Wilayah", fill = "Rerata") +
  theme_minimal()

# Tampilkan visualisasi berdampingan
plot_A + plot_B

Ilustrasi MIDP

Misalkan dalam menganalisis korelasi antara Tingkat Curah Hujan (mm) dan Jumlah Kasus Demam Berdarah (DBD) di suatu kota.

  • Agregasi Harian: Hubungan tidak terlihat atau korelasi sangat rendah (r≈0.05) karena ada time lag (masa inkubasi nyamuk butuh waktu beberapa hari/minggu setelah hujan).

  • Agregasi Mingguan: Hubungan mulai terlihat (r≈0.45).

  • Agregasi Bulanan: Korelasi menjadi sangat kuat (r≈0.82) karena variasi harian (noise) sudah terratakan.

Poin Inti: Fenomena fenomena spasial-temporal yang sama dapat menghasilkan kesimpulan statistik yang bertolak belakang hanya karena perubahan skala interval waktu pengamatan.

library(lubridate)
set.seed(42)

# Buat Data Harian selama 1 Tahun (365 Hari)
data_harian <- data.frame(
  tanggal = seq(as.Date("2024-01-01"), as.Date("2024-12-31"), by = "day")
) %>%
  mutate(
    bulan = floor_date(tanggal, "month"),
    # Simulasi Curah Hujan Harian
    curah_hujan = round(rpois(n(), lambda = 20) + sin(1:n()/30)*10, 1),
    # Kasus DBD dipengaruhi hujan tetapi ada noise harian yang tinggi
    kasus_dbd = round(curah_hujan * 0.3 + rnorm(n(), mean = 10, sd = 8), 0)
  )
head(data_harian)
     tanggal      bulan curah_hujan kasus_dbd
1 2024-01-01 2024-01-01        26.3        15
2 2024-01-02 2024-01-01        17.7        16
3 2024-01-03 2024-01-01        21.0        14
4 2024-01-04 2024-01-01        16.3        28
5 2024-01-05 2024-01-01        20.7        22
6 2024-01-06 2024-01-01        28.0        42
# Agregasi ke Level Bulanan
data_bulanan <- data_harian %>%
  group_by(bulan) %>%
  summarise(
    hujan_mean = mean(curah_hujan),
    dbd_mean   = mean(kasus_dbd),
    .groups    = "drop"
  )
head(data_bulanan)
# A tibble: 6 × 3
  bulan      hujan_mean dbd_mean
  <date>          <dbl>    <dbl>
1 2024-01-01      26.0      16.6
2 2024-02-01      30.5      18.7
3 2024-03-01      25.3      19  
4 2024-04-01      14.9      13.8
5 2024-05-01       9.89     13.0
6 2024-06-01      13.2      15.7
# 3. Hitung Korelasi
r_harian  <- cor(data_harian$curah_hujan, data_harian$kasus_dbd)
r_harian
[1] 0.3054636
r_bulanan <- cor(data_bulanan$hujan_mean, data_bulanan$dbd_mean)
r_bulanan
[1] 0.8511978
# 4. Visualisasi Scatter Plot
p_harian <- ggplot(data_harian, aes(x = curah_hujan, y = kasus_dbd)) +
  geom_point(alpha = 0.4, color = "darkblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Skala Harian", subtitle = paste0("r = ", round(r_harian, 2))) +
  theme_minimal()

p_bulanan <- ggplot(data_bulanan, aes(x = hujan_mean, y = dbd_mean)) +
  geom_point(size = 3, color = "firebrick") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Skala Bulanan (Agregat)", subtitle = paste0("r = ", round(r_bulanan, 2))) +
  theme_minimal()

p_harian + p_bulanan