Praktikum 1 : Pengenalan R untuk Analisis Statistika Spasial

Author

Rosita Ria Rusesta - Statistika dan Sains Data, IPB

Apa itu Data Spasial?

Perbedaan utama antara data spasial dan bukan spasial, adalah adanya atribut lokasi yang melekat pada data. Data spasial menjawab pertanyaan “Where is it? (dimana lokasi/posisinya)” dan atribut non spasialnya menjawab “what is it?” atau “how much/what condition?” (Apa objeknya, berapa jumlahnya, atau bagaimana sifatnya?).

Menurut Tobler pada tahun 1979 mengungkapkan di dalam Hukum Geografi pertama, bahwa “Segala sesuatu saling berhubungan satu dengan yang lainnya, tetapi sesuatu yang dekat lebih mempunyai pengaruh daripada sesuatu yang jauh”.

Jenis Data Spasial

Menurut Cressie & Moores (2021), suatu pemodelan proses spasial dibedakan menjadi tiga, yaitu:

  1. Point process (titik)

  2. Lattice process (data bersifat diskret)

  3. Geostatistical process (data bersifat kontinu)

Sedangkan, berdasarkan penyimpanannya, data spasial (lebih spesifiknya, data GIS) secara umum dibedakan menjadi dua jenis, yaitu:

  1. Vektor, terbagi menjadi tiga yaitu:

    • data titik

    • data garis (line / arc)

    • data poligon

  2. Data Raster

Paket R untuk Data Spasial

Dalam analisis data spasial di R, terdapat berbagai paket yang memiliki fungsi berbeda sesuai jenis data yang digunakan, baik itu data vektor, raster, maupun titik. Untuk mempermudah pengolahan, visualisasi, hingga analisis statistik spasial, berikut adalah beberapa paket utama yang perlu diaktifkan dengan perintah library(), sebagaimana ditunjukkan pada potongan kode berikut.

library(sp) # mengolah data peta (titik, garis, poligon) (versi lama)
library(sf) # mengolah data spasial vektor modern
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(mapview) # membuat peta interaktif
library(spData) # berisi contoh data spasial
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
library(spatstat) # untuk data berupa titik ()
Loading required package: spatstat.data
Loading required package: spatstat.univar
spatstat.univar 3.1-6
Loading required package: spatstat.geom
spatstat.geom 3.7-0
Loading required package: spatstat.random
spatstat.random 3.4-4
Loading required package: spatstat.explore
Loading required package: nlme
spatstat.explore 3.7-0
Loading required package: spatstat.model
Loading required package: rpart
spatstat.model 3.6-1
Loading required package: spatstat.linnet
spatstat.linnet 3.4-1

spatstat 3.5-1 
For an introduction to spatstat, type 'beginner' 
library(dplyr) # manipulasi dan pengolahan data

Attaching package: 'dplyr'
The following object is masked from 'package:nlme':

    collapse
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(ggplot2) # membuat grafik
library(rnaturalearth) # menyediakan data peta vektor(base map)
library(terra) # untuk data raster
terra 1.9.27

Attaching package: 'terra'
The following objects are masked from 'package:spatstat.geom':

    area, delaunay, is.empty, rescale, rotate, shift, where.max,
    where.min
library(raster) # versi lama mengolah data raster

Attaching package: 'raster'
The following object is masked from 'package:dplyr':

    select
The following object is masked from 'package:nlme':

    getData
library(rinat) # untuk data iNaturalist

1. Data Titik (point patterns)

Data titik (point pattern) digunakan untuk merepresentasikan fenomena spasial berdasarkan posisi lokasinya. Contoh fenomena yang dapat direpresentasikan dengan pola titik antara lain titik kebakaran hutan, lokasi pohon, atau titik kemunculan satwa.

Ilustrasi 1 : clmfires

Contoh dari data titik adalah kejadian kebakaran di Castilla-La Mancha, Spanyol, antara tahun 1998 dan 2007 yang terdapat dalam data clmfires pada paket spatstat (Baddeley, Turner, dan Rubak 2022). Data clmfires merupakan pola titik yang berisi informasi dari setiap kejadian kebakaran.

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
plot(clmfires, use.marks = FALSE, pch = ".",cols = "red")

Ilustrasi 2 : Data iNaturalist (Orangutan Kalimantan : Pongo pygmaeus)

Fitur iNaturalist dapat diakses langsung melalui R menggunakan package rinat.

# 1. Mengambil data dari iNaturalist
orangutan <- get_inat_obs(taxon_name = "Pongo pygmaeus",
                          bounds = c(-4, 108, 7, 119),
                          quality = "research",
                          maxresults = 1000)
head(orangutan,3)
       scientific_name                  datetime
1 Pongo pygmaeus morio 2026-08-05 06:48:03 +0800
2 Pongo pygmaeus morio 2026-08-16 14:49:00 +0800
3 Pongo pygmaeus morio 2026-08-16 13:07:00 +0800
                                                 description
1 Two found climbing around trees next to the Kinabatangan. 
2                                                           
3                                                           
          place_guess latitude longitude tag_list                 common_name
1           Sabah, MY 5.410289  118.3076          Northeast Bornean Orangutan
2 Sandakan, MY-SA, MY 5.923935  117.9758          Northeast Bornean Orangutan
3 Sandakan, MY-SA, MY 5.871293  117.9246          Northeast Bornean Orangutan
                                                 url
1 https://www.inaturalist.org/observations/392522275
2 https://www.inaturalist.org/observations/392462205
3 https://www.inaturalist.org/observations/392462161
                                                                   image_url
1 https://inaturalist-open-data.s3.amazonaws.com/photos/719359508/medium.jpg
2 https://inaturalist-open-data.s3.amazonaws.com/photos/719227982/medium.jpg
3 https://inaturalist-open-data.s3.amazonaws.com/photos/719227146/medium.jpg
   user_login        id               species_guess iconic_taxon_name taxon_id
1 rimbawan_37 392522275 Northeast Bornean Orangutan          Mammalia   147760
2 ericcameron 392462205 Northeast Bornean Orangutan          Mammalia   147760
3 ericcameron 392462161 Northeast Bornean Orangutan          Mammalia   147760
  num_identification_agreements num_identification_disagreements
1                             1                                0
2                             3                                0
3                             3                                0
   observed_on_string observed_on        time_observed_at    time_zone
1 2026-08-05 06:48:03  2026-08-05 2026-08-04 22:48:03 UTC Asia/Kuching
2  2026/08/16 2:49 PM  2026-08-16 2026-08-16 06:49:00 UTC Asia/Kuching
3  2026/08/16 1:07 PM  2026-08-16 2026-08-16 05:07:00 UTC Asia/Kuching
  positional_accuracy public_positional_accuracy geoprivacy taxon_geoprivacy
1               44249                      44249   obscured         obscured
2                  16                      31368                    obscured
3                  24                      31368                    obscured
  coordinates_obscured positioning_method positioning_device user_id
1                 true                gps                gps  736477
2                 true                                       5349844
3                 true                                       5349844
     user_name              created_at              updated_at quality_grade
1   Jeremy Loh 2026-08-18 23:58:07 UTC 2026-08-19 01:00:52 UTC      research
2 Eric Cameron 2026-08-18 19:43:42 UTC 2026-08-19 01:02:36 UTC      research
3 Eric Cameron 2026-08-18 19:43:36 UTC 2026-08-19 01:02:41 UTC      research
   license sound_url oauth_application_id captive_cultivated
1 CC-BY-NC      <NA>                    2              false
2 CC-BY-NC      <NA>                   NA              false
3 CC-BY-NC      <NA>                   NA              false
# Filter / Hapus baris yang memiliki nilai NA pada koordinat
orangutan_clean <- orangutan %>%
  filter(!is.na(latitude) & !is.na(longitude)) %>%
  filter(geoprivacy == "" | is.na(geoprivacy) | geoprivacy == "open") %>%
  filter(positional_accuracy <= 1000 | is.na(positional_accuracy))
class(orangutan_clean)
[1] "data.frame"

Seperti terlihat pada keluaran di atas, objek orangutan_clean merupakan kelas data frame, dimana terdapat kolom koordinat bujur dan lintang di dalamnya. Namun, R belum mengenali objek tersebut sebagai objek spasial. Konversi data frame menjadi objek spasial dapat dilakukan menggunakan fungsi st_as_sf().

# Konversi ke objek spasial (sf)
orangutan_sf <- st_as_sf(orangutan_clean, 
                         coords = c("longitude", "latitude"), 
                         crs = 4326)
class(orangutan_sf)
[1] "sf"         "data.frame"

Visualisasi Data obyek orangutan_sf

world <- ne_countries(scale = "medium", returnclass = "sf")

ggplot() +
  geom_sf(data = world, fill = "#e0e0e0", color = "white") +
  geom_sf(data = orangutan_sf, color = "#d95f02", alpha = 0.7, size = 2) +
  coord_sf(xlim = c(108, 119), ylim = c(-4, 7), expand = FALSE) +
  theme_minimal() +
  labs(
    title = "Persebaran Observasi Orangutan Kalimantan",
    subtitle = "Data iNaturalist",
    x = "Bujur (Longitude)",
    y = "Lintang (Latitude)"
  )

2. Data Area/Lattice

Data area biasanya muncul ketika banyaknya kasus / kejadian dari suatu peubah yang diamati tersebut diagregasi pada suatu area. Misalnya, dalam epidemiologi spasial, lokasi individu dengan penyakit tertentu sering digabungkan dalam area administratif. Data ini dapat dianalisis untuk memahami pola geografis dan mengidentifikasi faktor risiko penyakit, dengan mempertimbangkan konfigurasi lingkungan dan faktor lain yang diketahui memengaruhi risiko penyakit (Moraga, 2018) 

Berikut ini contoh data area yang merepresentasikan jumlah kematian bayi mendadak di setiap kabupaten di North Carolina, AS, pada tahun 1974 dari paket sf (Pebesma, 2022).

d <- st_read(system.file("shape/nc.shp", package = "sf"),
             quiet = TRUE)
mapview(d, zcol = "SID74")

Berbeda dengan pola titik yang merepresentasikan posisi individu, pada data areal analisis spasial dilakukan pada unit wilayah dengan batas yang jelas.

d <- st_read(system.file("shape/nc.shp", package = "sf"),
             quiet = TRUE)
mapview(d, zcol = "SID74")  # SID74 (jumlah kasus kanker)

3. Data Geostatistik (Data Kontinu)

Dataset properties dalam paket R spData adalah contoh data spasial berkelas sf berisi 1.000 titik data lokasi apartemen di Kota Athena, Yunani, pada tahun 2017. Dataset ini menyajikan variabel prpsqm (Price Per Square Meter) sebagai variabel respons utama yang merupakan data geostatistik, serta beberapa atribut pendukung lainnya. Variabel prpsqm (Price Per Square Meter) dikategorikan sebagai data geostatistik karena harga per meter persegi tersebut mewakili fenomena kontinu yang pada dasarnya ada di seluruh permukaan wilayah Kota Athena, namun nilainya diamati melalui sampel di 1.000 titik koordinat spesifik.

mapview(properties, zcol = "prpsqm")

Data Raster

Fungsi rast() dapat digunakan untuk membuat dan membaca data raster. Fungsi ini writeRaster() memungkinkan kita untuk menulis data raster. Misalnya, di sini kita menggunakan rast() untuk membaca data raster yang mewakili ketinggian/elevasi di Luksemburg dari paket terra, dan menetapkannya ke sebuah objek r kelas SpatRaster.

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

Coordinate Reference System

Sistem referensi koordinat (CRS) data spasial menentukan asal dan satuan pengukuran koordinat spasial. CRS penting untuk manipulasi, analisis, dan visualisasi data spasial, dan memungkinkan untuk menangani beberapa data dengan mengubahnya menjadi CRS umum. Lokasi di Bumi dapat dirujuk menggunakan CRS yang tidak diproyeksikan (juga disebut geografis) atau yang diproyeksikan. CRS yang tidak diproyeksikan atau geografis menggunakan nilai lintang dan bujur untuk mewakili lokasi pada permukaan elipsoid tiga dimensi Bumi. CRS yang diproyeksikan menggunakan koordinat Cartesian untuk merujuk lokasi pada representasi dua dimensi Bumi.

Kode EPSG

CRS yang paling umum dapat ditentukan dengan memberikan kode EPSG (European Petroleum Survey Group). Proyeksi spasial umum dapat ditemukan di https://spatialreference.org/ref/.

Rincian proyeksi tertentu dapat diperiksa menggunakan st_crs() fungsi paket sf. Misalnya, kode EPSG 4326 mengacu pada proyeksi garis bujur/lintang WGS84.

st_crs("EPSG:4326")$Name
[1] "WGS 84"
st_crs("EPSG:4326")$epsg
[1] 4326

Transformasi CRS dengan sf dan terra

Fungsi sf::st_crs() dan terra::crs() memungkinkan kita untuk mendapatkan CRS data spasial. Fungsi-fungsi ini juga memungkinkan kita untuk menetapkan CRS ke data spasial dengan menggunakan st_crs(x) <- value jika x adalah sfobjek, dan crs(r) <- value jika r adalah raster.

Perhatikan bahwa menetapkan CRS tidak mengubah data, itu hanya mengubah label CRS. Kita bisa menggunakansf::st_transform() dan terra::project() untuk mengubah data sfatau raster, masing-masing, ke CRS baru.

Untuk data sf , kita dapat membaca dan mendapatkan CRS, lalu mengubah data tersebut menjadi CRS baru sebagai berikut:

pathshp <- system.file("shape/nc.shp", package = "sf")
map <- st_read(pathshp, quiet = TRUE)

# Get CRS
st_crs(map)
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]]
# Transform CRS
map2 <- st_transform(map, crs = "EPSG:4326")

Kita dapat menggunakan terra untuk membaca dan mendapatkan CRS, dan mengubah data menjadi CRS raster baru sebagai berikut:

pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)

# Get CRS
# crs(r)
# Transform CRS
r2 <- terra::project(r, "EPSG:2169")
# Get CRS
# crs(r2)
plot(r2)