Data Spasial

Author

Alfa Nugraha

Pengantar

Fenomena spasial1 secara umum dapat dikatakan sebagai objek diskret dengan batas wilayah yang jelas atau fenomena kontinu yang dapat diamati lokasinya namun tidak memiliki batasan ruang.

  • 1 spatial (adj.) 1840 (spacial is from 1838), “occupying space, characterized by space,” from Latin spatium + adjectival suffix -al (1); formed in English as an adjective to space (n.), that goes with temporal. The meaning “of or relating to space” is from 1857.

    • Objek spasial diskret dapat berupa sungai, jalan, negara, kota, atau daerah penelitian. Fenomena yang bersifat kontinu atau dikenal dengan istilah “spatial fields”, misalnya elevasi, temperatur, dan kualitas udara.

    Representasi objek spasial biasanya berbentuk (data) vector. Gambaran dari data ini berupa objek yang memiliki bentuk (shape), koordinat (geometry), dan peubah/variabel/atribut. Contoh dari gugus data vector yaitu peta batas negara di dunia (geometri) yang menyimpan nama dan jumlah populasi di dalamnya. Data spasial kontinu biasanya direprentasikan dengan struktur data raster.

    Representasi Sederhana Data Spasial

    Data spasial, secara sederhana, dapat dibangun dengan beberapa kombinasi tipe data dasar di R. Sebagai contoh dibuat peta dimana lokasi (diwakili oleh longitude dan latitude)2 dari 10 stasiun cuaca dengan simbol A-J beserta presipitasi tahunan.

  • 2 Longitude dan latitude adalah sistem koordinat yang umum digunakan secara global. Longitude atau garis bujur adalah garis yang membentang dari utara ke selatan bumi. Latitude atau garis lintang adalah garis mendatar dari khatulistiwa (0) ke kutub selatan (90), atau khatulistiwa (0) ke kutub utara (-90)

    • name <- LETTERS[1:10]

longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5,
               -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7,
              38.9, 36.2, 39, 41.6, 36.9)
stations <- cbind(longitude, latitude)

# Generated rainfall data
set.seed(0)
precip <- round((runif(length(latitude))*10)^3)

    Peta adalah plot data geospasial yang memiliki label dan objek grafis seperti skala dan legenda. Data spasial sendiri tidak sepenuhnya berupa peta

    Berikut ini titik lokasi stasiun cuaca dimana besarnya diameter titik menggambarkan proporsi dari curah hujan.

    psize <- 1 + precip/500

plot(stations, cex=psize, pch=20, col='red', main='Curah hujan (mm/tahun)')
# add names to plot
text(stations, name, pos=4)
# add a legend
breaks <- c(100, 250, 500, 1000)
legend.psize <- 1+breaks/500
legend("topright", legend=breaks, pch=20, pt.cex=legend.psize, col='red', bg='gray')

    Figure 1: 10 stasiun cuaca dengan nilai curah hujan per lokasi

    Data di atas dibangkitkan dengan “longitude, latitude” namun terkadang pada sebagian besar peta, garis lintang (Utara/Selatan) digunakan untuk sumbu vertikal dan garis bujur (Timur/Barat) untuk sumbu horizontal. Hal ini penting untuk dicatat, karena ini adalah sumber kesalahan yang sangat umum!

    Plot di tersebut dapat ditambahkan beberapa titik, garis, dan poligon:

    lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
x <- cbind(lon, lat)
plot(stations, main='Curah hujan (mm/tahun)')
polygon(x, col='blue', border='light blue')
lines(stations, lwd=3, col='red')
points(x, cex=2, pch=20)
points(stations, cex=psize, pch=20, col='red', main='Curah hujan (mm/tahun)')

    Figure 2: Data vector

    Data spasial juga dapat dibentuk dengan sebuah matriks atau array berordo tinggi. Pendekatan yang sangat praktis dapat diimplementasikan pada tipe data titik. Kode berikut sebagai contoh dalam membuat gugus data spasial titik dengan atribut yang dapat dikombinasikan dengan geometri beserta atributnya dalam suatu data frame.

    wst <- data.frame(longitude, latitude, name, precip)
kable(wst)
    longitude latitude name precip
    -116.7 45.3 A 721
    -120.4 42.6 B 19
    -116.7 38.9 C 52
    -113.5 42.1 D 188
    -115.5 35.7 E 749
    -120.8 38.9 F 8
    -119.5 36.2 G 725
    -113.7 39.0 H 843
    -113.7 41.6 I 289
    -110.7 36.9 J 249

    wst hanya sebuah data frame dan R tidak otomatis mengenal makna khusus pada dua kolom pertama. Oleh karena itu, diperlukan adanya sistem (referensi) koordinat sebagai salah satu aspek penting yang digunakan untuk memproyeksi data spasial. Untuk memfasilitasi hal tersebut beberapa R package telah tersedia untuk mengoperasikan data spasial seperti terra dan sf yang menggantikan package lama yaitu raster dan sp.

    Tipe Data Spasial

    Data Vector

    Tipe utama dari data vector adalah titik, garis, dan poligon. Ketiga tipe tersebut memiliki koordinat berpasangan \((x, y)\).

    • Titik adalah tipe yang paling sederhana. Setiap titik memiliki satu pasang koordinat dan \(n\) jumlah atribut. Misalnya titik lokasi dimana tikus ditangkap. Atribut yang disimpan yaitu waktu ditangkap, penangkap, spesies, dan informasi tentang habitatnya. Kombinasi beberapa titik atau struktur multi-point dapat dimungkinkan dengan atribut tunggal di dalamnya seperti semua kedai kopi di suatu kota yang dapat dianggap sebagai suatu geometri tunggal.

    • Garis sedikit lebih kompleks mengacu ke suatu garis tunggal atau gugus garis yang saling terhubung (polyline). Misalnya sungai dengan atribut penyusun dari banyak garis. Perbedaan garis dengan titik, garis adalah sekumpulan koordinat yang terurut (node) atau serupa dengan struktur multi-titik.

    • Poligon merupakan sekumpulan garis yang saling menutup pada setiap pasang ujungnya. Setiap poligon memiliki atribut, sama seperti tipe lainnya kompilasi dari banyak poligon dapat juga dihitung sebagai geometri tunggal. Contohnya Indonesia terdiri dari banyak pulau. Masing-masing pulau diwakili oleh poligon tunggal tapi bisa juga mewakili satu negara dengan bentuk multi-poligon tunggal.

    Membuat data vector

    terra membuat sejumlah object-oriented class yang dimulai dengan kode Spat. SpatVectoradalah objek kelas yang digunakan untuk jenis data vector sebagai representasi geometri dan atribut yang menggambarkannya.

    Titik
    longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)

    Konversi koordinat di atas dapat dilakukan dengan perintah vect

    pts <- vect(lonlat)
class(pts)
## [1] "SpatVector"
## attr(,"package")
## [1] "terra"

    Deskripsi variabel pts sebagai berikut

    pts
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 10, 0  (geometries, attributes)
##  extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
##  coord. ref. :
    geom(pts)
##       geom part      x    y hole
##  [1,]    1    1 -116.7 45.3    0
##  [2,]    2    1 -120.4 42.6    0
##  [3,]    3    1 -116.7 38.9    0
##  [4,]    4    1 -113.5 42.1    0
##  [5,]    5    1 -115.5 35.7    0
##  [6,]    6    1 -120.8 38.9    0
##  [7,]    7    1 -119.5 36.2    0
##  [8,]    8    1 -113.7 39.0    0
##  [9,]    9    1 -113.7 41.6    0
## [10,]   10    1 -110.7 36.9    0

    Dapat dilihat bahwa koordinat di atas digunakan sebagai extent, atau disebut spatial extent. Terdapat juga informasi coord. ref. atau CRS (coordinate reference system). Cara menambahkan CRS sebagai berikut

    crdref <- "+proj=longlat +datum=WGS84"
pts <- vect(lonlat, crs=crdref)
pts
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 10, 0  (geometries, attributes)
##  extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs

crs(pts)
## [1] "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"

    Atribut dapat ditambahkan pada SpatVector pada nomor baris yang sama sesuai geometrinya.

    # Generate df nilai curah hujan random dengan kuantitas yang sama sebanyak titik
precipvalue <- runif(nrow(lonlat), min=0, max=100)
df <- data.frame(ID=1:nrow(lonlat), precip=precipvalue)

# Combine dengan df
ptv <- vect(lonlat, atts=df, crs=crdref)
ptv
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 10, 2  (geometries, attributes)
##  extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       :    ID precip
##  type        : <int>  <num>
##  values      :     1  6.179
##                    2   20.6
##                    3  17.66
    Garis dan poligon

    Membuat data vector titik sedikit lebih kompleks dibanding garis atau poligon

    lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
lonlat
##      id part    lon  lat
## [1,]  1    1 -116.8 41.3
## [2,]  1    1 -114.2 42.9
## [3,]  1    1 -112.9 42.4
## [4,]  1    1 -111.9 39.8
## [5,]  1    1 -114.2 37.6
## [6,]  1    1 -115.4 38.3
## [7,]  1    1 -117.7 37.6

# membuat garis
lns <- vect(lonlat, type="lines", crs=crdref)
lns
##  class       : SpatVector 
##  geometry    : lines 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs
    # membuat poligon
pols <- vect(lonlat, type="polygons", crs=crdref)
pols
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs

    Jadi, dalam membuat SpatVector diperlukan koordinat (geometri), CRS, dan atribut. Data vector di atas dapat divisualisasikan dalam bentuk peta.

    plot(pols, las=1)
plot(pols, border='blue', col='yellow', lwd=3, add=TRUE)
points(pts, col='red', pch=20, cex=3)

    Figure 3: Garis dan poligon

    Memuat data vector

    Terdapat banyak sekali format data vector. Shapefile adalah format yang paling banyak digunakan dan sangat penting untuk mengetahui sejumlah file pendukungnya setidaknya 3-4 file utama dengan nama yang sama tapi ekstensi yang berbeda.

    Sebagai contoh, shapefile x harus memiliki file pada direktori yang sama dengan ektensi: x.shp, x.shx, x.dbf, biasanya dilengkapi dengan x.prj.

    Berikut ini adalah daftar file yang biasa digunakan untuk mendefinisikan shapefile

    ekstensi file konten
    .dbf Attribute information
    .shp Feature geometry
    .shx Feature geometry index
    .aih Attribute index
    .ain Attribute index
    .prj Coordinate system information
    .sbn Spatial index file
    .sbx Spatial index file 
    depok_vect <- vect("kota_depok_2020.shp")
depok_vect
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 17  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  source      : kota_depok_2020.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area    ADM3_EN ADM3_PCODE ADM3_REF ADM3ALT1EN
##  type        :      <num>      <num>      <chr>      <chr>    <chr>      <chr>
##  values      :     0.1951   0.001207       Beji  ID3276050                    
##                    0.2806   0.001583 Bojongsari  ID3276011                    
##                    0.2442   0.001271   Cilodong  ID3276031                    
##  ADM3ALT2EN    ADM2_EN ADM2_PCODE    ADM1_EN (and 7 more)
##       <chr>      <chr>      <chr>      <chr>             
##             Kota Depok     ID3276 Jawa Barat             
##             Kota Depok     ID3276 Jawa Barat             
##             Kota Depok     ID3276 Jawa Barat

    Visualisasi data vector

    Visualisasi data vector yaitu berupa peta dengan perintah sebagai berikut

    # pilih satu atribut representasi dari poligon
plot(depok_vect, "ADM3_EN")

    Figure 4: Kecamatan Kota Depok

    Dengan menggunakan dua atribut, kustomisasi pada plot juga dapat dilakukan dengan sintaks berikut

    par(mfrow=c(1,2))
m <- c(3.1, 3.1, 2.1, 2.1)
plot(depok_vect, "ADM3_EN", col=rainbow(25), mar=m, plg=list(x="topright", title="Kecamatan", bty = "o"), pax=list(las=1))
plot(depok_vect, "ADM2_EN", col=rainbow(25), mar=m, plg=list(x="topright", title="Kota Depok", bty = "o", cex=.75), pax=list(las=1))

    d <- aggregate(depok_vect, "ADM3_EN")
plot(depok_vect, col="light blue", lty=2, border="red", lwd=2)
lines(d, lwd=5)
lines(d, col="white", lwd=1)
text(depok_vect, "ADM3_EN", cex=.8, halo=TRUE)

    Figure 5: Kecamatan Kota Depok dengan label

    Menyimpan data vector

    outfile <- "shp_test.shp"
writeVector(depok_vect, outfile, overwrite=TRUE)

    Data Raster

    Jenis data raster umumnya digunakan untuk menggambarkan fenomena kontinu secara spasial seperti elevasi. Sebuah raster dibentuk dari sebuah grid yang bebentuk persegi, dikenal dengan cell atau dalam konteks penginderaan jauh (remote sensing), piksel yang mengandung suatu nilai tertentu atau lebih (bahkan NA).

    Nilai cell pada raster normalnya adalah rata-rata (atau mayoritas) luasan dari suatu wilayah. Namun, dalam beberapa studi kasus nilai tersebut adalah nilai center suatu cell. Cell tersebut dapat berkumpul menjadi suatu atribut.

    Berbeda dengan data vector, geometri pada data raster tidak secara eksplisit disimpan sebagai koordinat tetapi digunakan untuk mengetahui spasial extent atau jumlah baris dan kolom sebagai pembagian atau batas areanya. Dari informasi tersebut, ukuran cell pada raster dapat dikomputasikan.

    Membuat data raster

    Banyak sekali fungsi atau geoprocessing tools untuk analisis data raster pada terra. Contohnya fungsi untuk membaca sepotong data raster dari file atau mengubahnya nomor baris cell menjadi suatu koordinat atau sebaliknya. Package ini juga dapat melakukan operasi aljabar dan berbagai fungsi manipulasi data raster.

    SpatRaster digunakan sebagai gambaran data raster multi-layer (multi-variable) yang menyimpan sejumlah parameter dasar menggambarkan geometrinya termasuk CRS. SpatRaster juga dapat menyimpan informasi file raster yang dimuat atau menyimpannya ke dalam memory.

    r <- rast(ncol=10, nrow=10, xmin=-150, xmax=-80, ymin=20, ymax=60)
r
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84

    Variabel r di atas hanya menyimpan geometrinya saja. Untuk menambahkan nilai pada masing-masing cellnya dapat dilakukan dengan cara berikut

    values(r) <- runif(ncell(r))
r
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        :      lyr.1 
## min value   : 0.01307758 
## max value   :  0.9926841

    Atau dapat diisi dengan nomor cell

    values(r) <- 1:ncell(r)
r
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     1 
## max value   :   100

    Plot objek tersebut

    plot(r)

# tambah poligon dan titik 
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)

pts <- vect(lonlat)
pols <- vect(lonlat, type="polygons", crs="+proj=longlat +datum=WGS84")
points(pts, col="red", pch=20, cex=3)
lines(pols, col="blue", lwd=2)

    Multi-layer objek pada data raster dapat dilakukan dengan perintah c()

    r2 <- r * r
r3  <- sqrt(r)
s <- c(r, r2, r3)
s
## class       : SpatRaster 
## dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## sources     : memory  
##               memory  
##               memory  
## names       : lyr.1, lyr.1, lyr.1 
## min values  :     1,     1,     1 
## max values  :   100, 10000,    10

plot(s)

    Memuat data raster

    depok_rast <- rast("kota_depok.tif")
depok_rast
## class       : SpatRaster 
## dimensions  : 100, 75, 1  (nrow, ncol, nlyr)
## resolution  : 0.002694673, 0.001468689  (x, y)
## extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : kota_depok.tif 
## name        : ADM3_EN 
## min value   :    Beji 
## max value   :   Tapos

    Visualisasi data raster

    Visualisasi data raster yaitu berupa peta juga dengan perintah sebagai berikut

    plot(depok_rast)

    Plot di atas adalah plot data raster dengan label kategori sebagai legendanya. Ilustrasi lain yang lebih mewakili penggunaan data raster misalnya adalah elevasi

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

    SpatRaster dapat dikombinasikan dengan data vector (titik, garis, poligon).

    plot(elev)
lines(p, lwd=2)
set.seed(12)
xy <- spatSample(elev, 20, "random", na.rm=TRUE, xy=TRUE)
points(xy, pch=20, col="red", cex=2)

    Kemudian dapat menggunakan legenda dengan interval nilai elevasi tertentu

    m <- c(3.1, 3.1, 1.1, 1.1)
plot(elev, type="interval", plg=list(x="topright"), mar=m)

    Menyimpan data raster

    
# konversi data vector menjadi raster
depok_ext <- rast(depok, ncols=75, nrows=100)
depok_rast <- rasterize(depok, depok_ext, "ADM3_EN")

# save as raster
writeRaster(depok_rast, "kota_depok.tif", overwrite=TRUE)

    Coordinate Reference Systems

    Koordinat angular

    Bumi memiliki bentuk seperti bola yang tidak beraturan (geoid). Dalam menentukan posisi atau titik suatu wilayah, sistem yang sudah sangat umum digunakan adalah sistem koordinat dengan ukuran derajat longitude dan latitude.

    Pada dasarnya koordinat tersebut tidak sepenuhnya dapat dijadikan ukuran muka bumi tetapi perlu pendugaan model yang disebut datum. Datum adalah serangkaian paramater yang menentukan posisi asal, skala, dan orientasi suatu sistem koordinat.

    Datum yang lebih kompleks dapat memperkirakan berbagai variasi bentuk bumi (ellipsoid) berdasarkan proyeksi yang berasosiasi dengan lokasinya. Datum yang paling umum, banyak digunakan, dan cocok secara global adalah WGS84 (World Geodesic System 1984). Datum yang fit secara lokal dapat digunakan pada negara atau wilayah bumi bagian tertentu dengan aproksimasi error yang kecil.

    Proyeksi

    Pertanyaan utama pada analisis spasial dan kartografi adalah bagaimana melakukan transformasi sistem angular 3D ke dalam sistem planar 2D. Beberapa tipe planar mengacu pada projections. Proyeksi mengacu kepada perhitungan matematis yang dilakukan untuk meratakan data menjadi bentuk ruang 2D. Sistem referensi koordinat (CRS) pada titik koordinat \(x\) dan \(y\) yang berasosiasi dengan proyeksi dapat digunakan untuk mentransformasi bentuk 3D permukaan bumi.

    Sejauh ini belum ada datum yang benar-benar memproyeksikan secara tepat keseluruhan bentuk permukaan bumi. Beberapa hanya dapat digunakan pada area kecil. Contoh beberapa proyeksi yaitu “Mercator”, “UTM”, “Robinson”, “Lambert”, “Sinusoidal”, and “Albers”. Pada dua gugus data spasial yang sama namun berbeda proyeksi tidak akan sejajar saat dirender bersamaan3.

  • 3 Perhatikan perbedaan bentuk dari masing-masing proyeksi. Perbedaan ini adalah hasil kalkukasi transformasi data ke peta 2D. Sumber: M. Corey, opennews.org

    • Peta US dengan berbagai proyeksi berbeda

    Notasi

    Planar CRS didefinisikan dengan proyeksi, datum, dan sejumlah parameter yang bergantung pada masing-masing proyeksinya. Parameter tersebut digunakan untuk menentukan titik pusat pada peta. Sebagai bagian dari pengembangan open source geospatial software, PROJ.44 hadir sebagai dukungan dalam proyeksi kartografik. Pada R, notasi dan definisi PROJ.4 digunakan untuk proyeksi dan basis transformasi CRS.

  • 4 Beberapa referensi spasial dapat ditemukan di alamat berikut spatialreference.org

    • CRS yang paling umum digunakan telah ditetapkan sebagai kode EPSG (European Petroleum Survey Group). Unik ID ini dapat digunakan untuk mengidentifikasi sebuah CRS. Misalnya EPSG:27561 ekivalen dengan string PROJ.4:

    +proj=lcc +lat_1=49.5 +lat_0=49.5 +lon_0=0 +k_0=0.999877341 +x_0=6 +y_0=2 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs

    Semua kode unik EPSG dapat dilihat pada fungsi built-in proj4string

    kable(head(rgdal::make_EPSG(), 10))
    code note prj4 prj_method
    3819 HD1909 +proj=longlat +ellps=bessel +no_defs +type=crs (null)
    3821 TWD67 +proj=longlat +ellps=aust_SA +no_defs +type=crs (null)
    3822 TWD97 +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs (null)
    3823 TWD97 +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
    3824 TWD97 +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
    3887 IGRS +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs (null)
    3888 IGRS +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
    3889 IGRS +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
    3906 MGI 1901 +proj=longlat +ellps=bessel +no_defs +type=crs (null)
    4000 MOLDREF99 +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs (null)

    Menetapkan CRS

    Terkadang data spasial yang dimiliki belum memilik CRS karena data tersebut dibuat dengan dari data tabular biasa sehingga perlu menetapkan CRS pada data tersebut.

    # empty the crs original first
crs(pts) <- ""

# set new crs
crs(pts) <- "+proj=longlat +datum=WGS84"
cat(crs(pts))
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

    Cara di atas pada dasarnya bukan cara yang tepat untuk menetapkan CRS pada data spasial. Oleh karena itu, perlu adanya transformasi data spasial

    Transformasi data vector

    Transformasi data dengan CRS lain dapat menggunakan fungsi project. Misalnya digunakan proyeksi Robinson.

    newcrs <- "+proj=robin +datum=WGS84"
    rob <- terra::project(pts, newcrs)
rob
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : -10366659, -9704732, 4020871, 4583692  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

    Setelah transformasi, unit geometri tidak lagi dalam satuan derajat tetapi dalam meter, begitu juga pada spatial extent yang digunakan. Transformasi balik dapat juga dilakukan

    pts2 <- terra::project(rob, "+proj=longlat +datum=WGS84")

    Transformasi data raster

    Data vector data ditransform dari koordinat lon/lat ke planar dan sebaliknya tanpa kehilangan presisi. Namun hal ini tidak untuk data raster. Raster terdiri dari cell persegi panjang yang berukuran sama sehingga tidak mungkin ditransformasi cell per cell.

    Masing-masing nilai pada cell baru perlu diduga berdasarkan nilai yang overlap pada cell yang lama. Pada kasus data kategorik, nearest neighbor dapat digunakan selainnya menggunakan interpolasi.

    Karena mengubah proyeksi raster dapat mengubah nilai suatu cell, pada banyak kasus hal ini harus dihindari namun dapat dilakukan

    r <- rast(xmin=-110, xmax=-90, ymin=40, ymax=60, ncols=40, nrows=40)
values(r) <- 1:ncell(r)
r
## class       : SpatRaster 
## dimensions  : 40, 40, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -110, -90, 40, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        : lyr.1 
## min value   :     1 
## max value   :  1600
## class       : SpatRaster
## dimensions  : 40, 40, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -110, -90, 40, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s)   : memory
## name        : lyr.1
## min value   :     1
## max value   :  1600
plot(r)

    # menggunakan proyeksi Robinson
pr1 <- terra::project(r, "+proj=robin +datum=WGS84")
crs(pr1)
## [1] "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"unknown\",\n        METHOD[\"Robinson\"],\n        PARAMETER[\"Longitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\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[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

plot(pr1)

    Cara di atas bukanlah metode yang tepat karena seharusnya parameter raster yang diproyeksikan harus sesuai atau sejajar dengan raster asal.

    Cara terbaik dalam memproyeksikan data raster adalah dengan menggunakan SpatRaster yang memiliki proyeksi yang dibutuhkan. Pada contoh ini, raster proyeksi harus dibuat terlebih dahulu.

    x <- rast(pr1)
# Set the cell size
res(x) <- 200000

pr3 <- terra::project(r, x)
pr3
## class       : SpatRaster 
## dimensions  : 10, 14, 1  (nrow, ncol, nlyr)
## resolution  : 2e+05, 2e+05  (x, y)
## extent      : -9577685, -6777685, 4283463, 6283463  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        :    lyr.1 
## min value   : 111.1541 
## max value   :  1523.58

plot(pr3)

    Analisis data raster sangat penting menggunakan proyeksi yang sama terutama saat area yang digunakan besar. Hal ini untuk memastikan grid cell berukuran sama dan dapat dikomparasi satu sama lain, khususnya saat perhitungan luasan data digunakan.