Modul ini disusun dengan merujuk pada beberapa referensi yang tertulis di bagian akhir dari modul ini. Secara umum, modul ini akan memberikan beberapa ilustrasi terkait dengan topik berikut:

  1. Pengenalan struktur data spasial
  2. Pengenalan R untuk impor data spasial
  3. Pengenalan R untuk manajemen data spasial
  4. Pengenalan R untuk visualisasi data spasial

Struktur Data Spasial

Perbedaan utama antara data spasial dan bukan spasial, adalah adanya atribut lokasi yang melekat pada data. Informasi tersebut dapat berupa informasi titik lokasi maupun area. 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
Gambar 1. Konstruksi jenis data vektor (Slingsby, 2023)
Gambar 1. Konstruksi jenis data vektor (Slingsby, 2023)

 

Data titik menyimpan informasi lokasi menurut coordinate reference system (CRS), data garis merupakan beberapa titik yang dihubungkan oleh garis sehingga membentuk informasi jarak/panjang, sedangkan data poligon berupa suatu bidang yang memiliki luasan area. Data poligon seringkali merepresentasikan suatu peubah yang bersifat diskret.

Gambar 2. Contoh visualisasi data poligon yang merepresentasikan jenis vegetasi di Cape Peninsula (Slingsby, 2023)
Gambar 2. Contoh visualisasi data poligon yang merepresentasikan jenis vegetasi di Cape Peninsula (Slingsby, 2023)

 

  1. Raster, dikenal juga sebagai data grid, bersifat pixel-based, serta umumnya merepresentasikan data yang bersifat kontinu, contohnya kualitas udara, temperatur, presipitasi, elevasi, dsb.
Gambar 3. Contoh visualisasi data raster yang merepresentasikan digital elevaton model (DEM) di Cape Peninsula (Slingsby, 2023)
Gambar 3. Contoh visualisasi data raster yang merepresentasikan digital elevaton model (DEM) di Cape Peninsula (Slingsby, 2023)

 

Paket R untuk Data Spasial

Secara umum, terdapat dua paket utama pada lingkungan R yang dapat digunakan pada data spasial, yaitu paket sp dan sf. Paket sp dikembangkan lebih dulu, sekitar tahun 2005, sedangkan paket sf dikembangkan sekitar tahun 2016. Oleh karenanya, paket sf dianggap lebih compatible dengan lingkungan tidymodels.

Paket sf konon memang dirancang untuk menggantikan sp, namun hingga saat ini masih cukup banyak metode pemodelan dan analisis yang bergantung pada paket sp. Oleh karenanya, beberapa sumber menyarankan agar kita tetap mempelajari keduanya, dan mungkin perlu membiasakan untuk mengonversi objek sf menjadi sp atau sebaliknya.

Paket lainnya yang akan digunakan adalah raster, untuk menangani data jenis raster/grid. Serupa dengan kondisi sebelumnya, terdapat paket lain yaitu terra yang memiliki fitur yang mirip dengan raster namun dianggap memiliki lebih banyak pengembangan dan pembaruan.

Selain dari paket yang telah disebutkan di atas, terdapat beberapa paket lain yang cukup banyak digunakan untuk memproses data spasial, masing-masing spesifik dengan kebutuhan analisis yang akan dilakukan. Untuk lebih lengkapnya, kita dapat melihatnya pada laman CRAN Task View: Analysis of Spatial Data, pada link berikut: https://cran.r-project.org/web/views/Spatial.html.

Operasi GIS pada data vektor

Impor dan Ekspor data

Membaca data spasial berjenis vektor dapat dilakukan menggunakan fungsi st_read(). Silahkan mengunduh data yang digunakan pada ilustrasi ini pada link berikut: klik disini

Data spasial dapat disimpan dalam beberapa format, seperti shp, json, tiff, dll. Namun jenis shp (shapefiles) termasuk yang paling banyak digunakan. Data dengan format shapefile akan memerlukan beberapa file dengan nama yang sama, di antaranya adalah file dengan format .shp, .shx, .dbf, .prj, dsb. Sehingga kita perlu menyimpan seluruh file tersebut dalam satu folder. Folder tersebut selanjutnya disebut sebagai dsn, sedangkan nama file disebut sebagai layer. Selanjutnya, impor data dapat dilakukan dengan fungsi seperti bawah ini.

<nama object> <-st_read(dsn="<nama folder>", layer="<nama file>")

Berikut contohnya:

library(sf)
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
veg<-st_read(dsn="cape_peninsula/veg", layer="Vegetation_Indigenous")
## Reading layer `Vegetation_Indigenous' from data source 
##   `D:\Dept.STK\Courses\Statistika Spasial (S2)\STA 1553 TA 23-24, genap\P1\cape_peninsula\veg' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1325 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -63972.95 ymin: -3803535 xmax: 430.8125 ymax: -3705149
## Projected CRS: WGS_1984_Transverse_Mercator

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

st_crs(veg)
## Coordinate Reference System:
##   User input: WGS_1984_Transverse_Mercator 
##   wkt:
## PROJCRS["WGS_1984_Transverse_Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["Hartebeesthoek94",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6148]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",19,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

Secara umum, terdapat tiga komponen utama dari informasi CRS di atas, yaitu:

  • BASEGEOGCRS, informasi CRS geografis (unprojected)

  • CONVERSION, informasi proyeksi, keluaran di atas memperlihatkan bahwa proyeksi menggunakan Tranverse Mercator Lo19 ("Longitude of natural origin",19), yang artinya bahwa data berpusat pada 19 derajat pada garis bujur.

  • CS, sumbu kartesius yang digunakan. Pada keluaran di atas, terlihat bahwa digunakan sumbu dengan arah utara dan timur dengan satuan meter.

Selanjutnya, kita juga dapat memerika kelas dari objek veg tersebut.

class(veg)
## [1] "sf"         "data.frame"

Karena objek tersebut mengandung dua kelas, maka kita dapat menerapkan fungsi yang compatible untuk keduanya pada objek tersebut.

head(veg)

Kolom geometry mengimpan informasi terkait dengan ukuran geometri pada suatu objek sf, contohnya titik, garis, poligon, yang bersesuaian dengan baris pada atribut data pada kolom-kolom yang berada di depannya.

Seandainya kita perlu menyimpan suatu objek sf ke dalam sebuah file, maka kita dapat menggunakan fungsi st_write(), seperti pada contoh berikut ini.

st_write(veg, "cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp", append=F)
## Deleting layer `Vegetation_Indigenous_duplicate' using driver `ESRI Shapefile'
## Writing layer `Vegetation_Indigenous_duplicate' to data source 
##   `cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp' using driver `ESRI Shapefile'
## Writing 1325 features with 5 fields and geometry type Polygon.

Jika kita ingin menimpa file dengan nama yang sama (jika sudah pernah disimpan), maka kita perlu mengganti argumen append=TRUE.

Visualisasi Sederhana

Cara paling mudah untuk menghasilkan visualisasi adalah dengan fungsi plot() seperti pada contoh berikut ini.

plot(veg)

Terlihat pada keluaran di atas bahwa terdapat 5 plot yang dihasilkan. Hal ini karena objek sf secara default akan menghasilkan plot untuk setiap kolom atribut pada data. Jadi, kita perlu memilih salah satu atribut, sesuai dengan kebutuhan, untuk diplotkan.

plot(veg[3])

Kita dapat pula membuat visualisasi menggunakan ggplot, seperti pada contoh di bawah ini.

library(tidyverse)

ggplot() +
  geom_sf(data=veg, aes(fill =`National_`))

Cropping Data

Seandainya kita ingin melakukan cropping, atau membuang sebagian pengamatan, sesuai dengan batas koordinat tertentu, maka kita bisa menggunakan fungsi st_crop() untuk melakukannya. Contohnya, misalnya kita menginginkan wilayah dengan batas sesuai dengan nilai yang disimpang pada objek ext berikut ini.

ext <- c(-66642.18, -3809853.29, -44412.18, -3750723.29) 
names(ext) <- c("xmin", "ymin", "xmax", "ymax") 

ext
##        xmin        ymin        xmax        ymax 
##   -66642.18 -3809853.29   -44412.18 -3750723.29

Selanjutnya kita akan crop sebagian data pada objek veg berdasarkan batas yang telah ditetapkan pada objek ext.

veg <- st_crop(veg, ext)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
ggplot() + geom_sf(data=veg, aes(fill = `National_`))

Memilih subset dari data berdasarkan atribut

Sebagai ilustrasi, misalnya kita mengingkan subset dari data menurut atribut/kolom/peubah National_ seperti di bawah ini.

#Make a vector of the veg types we want
split_veg <- c("Peninsula Granite Fynbos - North", 
               "Peninsula Granite Fynbos - South", 
               "Cape Flats Dune Strandveld - West Coast", 
               "Cape Flats Dune Strandveld - False Bay")

Selanjutnya subset dapat diambil menggunakan tidyverse seperti pada contoh dibawah ini, lalu dibuat plot dari subset tersebut.

#Using tidyverse piping to filter and plot
veg %>% 
  filter(National_ %in% split_veg) %>%
  ggplot() +
  geom_sf(aes(fill = `National_`))

Menggabungkan beberapa kategori dan memvisualisasikan hasil agregasi data

Misalnya kita ingin menggabungkan kategori Peninsula Granite Fynbos - North dan Peninsula Granite Fynbos - South menjadi satu yaitu Peninsula Granite Fynbos, sedangkan dua kategori lainnya menjadi Cape Flats Dune Strandveld.

veg$National_ <- str_replace_all(veg$National_,
                  c("Peninsula Granite Fynbos - North" = "Peninsula Granite Fynbos", 
                  "Peninsula Granite Fynbos - South" = "Peninsula Granite Fynbos", 
                  "Cape Flats Dune Strandveld - West Coast" = "Cape Flats Dune Strandveld", 
                  "Cape Flats Dune Strandveld - False Bay" = "Cape Flats Dune Strandveld"))
veg %>% group_by(National_) %>% 
  summarize() %>% 
  ggplot() + geom_sf(aes(fill = National_))

Ilustrasi menggunakan fitur iNaturalist (data titik)

Fitur iNaturalist dapat diakses langsung melalui R menggunakan package rinat.

library(rinat)

#Call the data directly from iNat
pc <- get_inat_obs(taxon_name = "Protea cynaroides",
                   bounds = c(-35, 18, -33.5, 18.5),
                   maxresults = 1000)

#View the first few rows of data
head(pc)
#Filter returned observations by a range of column attribute criteria
pc <- pc %>% filter(positional_accuracy<46 & 
                latitude<0 &
                !is.na(latitude) &
                captive_cultivated == "false" &
                quality_grade == "research")

class(pc)
## [1] "data.frame"

Seperti terlihat pada keluaran di atas, objek pc 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().

#Make the dataframe a spatial object of class = "sf"
pc <- st_as_sf(pc, coords = c("longitude", "latitude"), crs = 4326)

class(pc)
## [1] "sf"         "data.frame"
names(pc)
##  [1] "scientific_name"                  "datetime"                        
##  [3] "description"                      "place_guess"                     
##  [5] "tag_list"                         "common_name"                     
##  [7] "url"                              "image_url"                       
##  [9] "user_login"                       "id"                              
## [11] "species_guess"                    "iconic_taxon_name"               
## [13] "taxon_id"                         "num_identification_agreements"   
## [15] "num_identification_disagreements" "observed_on_string"              
## [17] "observed_on"                      "time_observed_at"                
## [19] "time_zone"                        "positional_accuracy"             
## [21] "public_positional_accuracy"       "geoprivacy"                      
## [23] "taxon_geoprivacy"                 "coordinates_obscured"            
## [25] "positioning_method"               "positioning_device"              
## [27] "user_id"                          "user_name"                       
## [29] "created_at"                       "updated_at"                      
## [31] "quality_grade"                    "license"                         
## [33] "sound_url"                        "oauth_application_id"            
## [35] "captive_cultivated"               "geometry"

Selanjutnya, kita dapat membuat visualisasi dari objek pc, yang telah memasukkan informasi spasial di dalamnya.

#Plot
ggplot() + geom_sf(data=pc)

Lebih jauh lagi, kita dapat menambahkan base map pada visualisasi di atas, untuk memudahkan interpretasi.

#install.packages(c("rosm", "ggspatial", "prettymapr"))
library(rosm)
library(ggspatial)

ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data=pc)

Membuat peta interaktif

#install.packages(c("leaflet", "htmltools"))
library(leaflet)
library(htmltools)
leaflet() %>%
  # Add default OpenStreetMap map tiles
  addTiles(group = "Default") %>%  
  # Add our points
  addCircleMarkers(data = pc,
                   group = "Protea cynaroides",
                   radius = 3, 
                   color = "green")

Operasi GIS pada data raster

Impor data

#install.packages("terra")

library(terra)
## terra 1.7.65
## 
## Attaching package: 'terra'
## The following object is masked from 'package:tidyr':
## 
##     extract
dem <- rast("cape_peninsula/CoCT_10m.tif")

class(dem)
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
dem
## class       : SpatRaster 
## dimensions  : 9902, 6518, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : -64180, 1000, -3804020, -3705000  (xmin, xmax, ymin, ymax)
## coord. ref. : GCS_WGS_1984 
## source      : CoCT_10m.tif 
## name        : 10m_BA 
## min value   :    -35 
## max value   :   1590

Untuk memastikan informasi tentang CRS dan proyeksi yang digunakan pada objek dem, kita dapat menggunakan fungsi crs().

crs(dem)
## [1] "PROJCRS[\"GCS_WGS_1984\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.25722356049,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"Transverse Mercator\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",19,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",1,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

Cropping

dem <- crop(dem, ext(c(-66642.18, -44412.18, -3809853.29, -3750723.29)))

dem
## class       : SpatRaster 
## dimensions  : 5330, 1977, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : -64180, -44410, -3804020, -3750720  (xmin, xmax, ymin, ymax)
## coord. ref. : GCS_WGS_1984 
## source(s)   : memory
## varname     : CoCT_10m 
## name        : 10m_BA 
## min value   :    -15 
## max value   :   1084

Aggregating / Resampling

Misalnya, kita ingin melakukan agregrasi data, yaitu dengan menggabungkan 9 pixel (3x3) menjadi 1 pixel. Kita dapat menggunakan fungi aggregate() seperti contoh di bawah ini.

dem30 <- aggregate(dem, fact = 3, fun = mean)
## 
|---------|---------|---------|---------|
=========================================
                                          
dem30
## class       : SpatRaster 
## dimensions  : 1777, 659, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -64180, -44410, -3804030, -3750720  (xmin, xmax, ymin, ymax)
## coord. ref. : GCS_WGS_1984 
## source(s)   : memory
## name        :   10m_BA 
## min value   :  -15.000 
## max value   : 1083.556

Basic Plotting

plot(dem30)

Atau bisa juga menggunakan ggplot seperti pada contoh di bawah ini.

#call tidyverse libraries and plot
library(tidyverse)

names(dem30)
## [1] "10m_BA"
dem30 %>% as.data.frame(xy = TRUE) %>%
  ggplot() +
  geom_raster(aes(x = x, y = y, fill = `10m_BA`)) + coord_fixed()

Berikut ini adalah plot dari data awal, sebelum dilakukan agregasi.

dem %>%
  as.data.frame(xy = TRUE) %>%
  ggplot() +
  geom_raster(aes(x = x, y = y, fill = `10m_BA`)) + coord_fixed()

Visualisasi beberapa gugus data dalam satu peta

Sebelum menggabungkan beberapa data, berikut ini disajikan contoh jika kita ingin mengubah data raster menjadi titik.

library(rinat)
library(sf)

#Call data for two species directly from iNat
pc <- get_inat_obs(taxon_name = "Protea cynaroides",
                   bounds = c(-35, 18, -33.5, 18.5),
                   maxresults = 1000)

ll <- get_inat_obs(taxon_name = "Leucadendron laureolum",
                   bounds = c(-35, 18, -33.5, 18.5),
                   maxresults = 1000)

#Combine the records into one dataframe
pc <- rbind(pc,ll)

#Filter returned observations by a range of attribute criteria
pc <- pc %>% filter(positional_accuracy<46 & 
                latitude<0 &
                !is.na(latitude) &
                captive_cultivated == "false" &
                quality_grade == "research")

#Make the dataframe a spatial object of class = "sf"
pc <- st_as_sf(pc, coords = c("longitude", "latitude"), crs = 4326)

#Set to the same projection as the elevation data
pc <- st_transform(pc, crs(dem30))

Sebelum menggabungkan beberapa data, berikut ini disajikan contoh jika kita ingin mengubah data raster menjadi vektor.

#Make the vegetation type a factor
veg$National_ <- as.factor(veg$National_)

#Rasterize
vegras <- rasterize(vect(veg), dem30, field = "National_")

#Plot
vegras %>% 
  as.data.frame(xy = TRUE) %>%
  ggplot() +
  geom_raster(aes(x = x, y = y, fill = National_)) + coord_fixed()

Selanjutnya, kita gabungkan beberapa objek yang telah kita simpan sebelumnya, ke dalam satu peta.

ggplot() +
  geom_raster(data = as.data.frame(vegras, xy = TRUE),
              aes(x = x, y = y, fill = National_)) +
  geom_contour(data = as.data.frame(dem30, xy = TRUE), 
               aes(x = x, y = y, z = `10m_BA`), breaks = seq(0, 1100, 100), colour = "black") +
  geom_sf(data=pc, colour = "white", size = 0.5)

Ada banyak hal yang dapat dicoba lebih lanjut tentang bagaimana membaca dan mengolah data spasial menggunakan R. Silahkan mempelajari hal-hal yang tidak tercakup pada modul ini pada referensi yang tercantum di bawah ini. Semoga bermanfaat!

Referensi

Brazil, N. (2023). Lab 3: Spatial Data in R. https://crd230.github.io/lab3.html (accessed at 23 January 2024)

Slingsby, J. (2023). A Minimal Introduction to GIS (in R). Retrieved from: https://www.ecologi.st/spatial-r/r-as-a-gis.html. (accessed at 24 January 2024)

https://tmieno2.github.io/R-as-GIS-for-Economists/before-you-start-1.html