Spatial projection in R

Membaca dan merubah proyeksi atau sistem kordinat dengan program R

Istilah "Proyeksi" pasti sering Anda dengar ketika banyak bekerja dengan data yang memiliki informasi geografis. Proyeksi dan sistem kordinat adalah suatu sistem yang dapat mentranformasi data dari bentuk sphere (rupa muka bumi) yang berbentuk bulat menjadi peta dua dimensi melalui pendekatan matematis.

Pada program R, kedua informasi proyeksi dan sistem kordinat disimpan dalam format yang spesifik yakni proj4 string yang disimpan dalam pustaka PROJ.4 (http://trac.osgeo.org/proj). Berbagai sistem proyeksi berikut definisinya dijelaskan secara komprehensif pada pustaka tersebut.

Sebagai contoh, sistem kordinat geografi (tanpa proyeksi) dapat dituliskan seperti ini
"+proj=longlat +datum=wgs84"

Sementara untuk proyeksi UTM seperti pada wilayah Jawa Barat (WGS 84/UTM Zona 48S) dapat dituliskan seperti:
"+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Pada tutorial kali ini, kita akan mencoba membaca dan merubah proyeksi data spasial dengan berbagai format data. Kita akan menggunakan data yang sama pada tutorial sebelumnya yang akan menghasilkan dua data spasial yang memiliki proyeksi yang sama seperti gambar dibawah.

Data yang digunakan dalam tutorial ini dapat diunduh disini.

Gambar 1. Peta curah hujan tahunan (biru) di sekitar TNGP (merah)

Seperti biasa, kita harus menyiapkan lokasi atau folder aktif dimana hasil ataupun output akan disimpan, dengan menggunakan fungsi setwd pada R.

## Mengaktifkan folder kerja, dimana semua hasil luaran akan disimpan
setwd("D:/LUBIS_PRIVATE_DATA/R_GIS_PROJECT/04_Hasil")

Kemudian memanggil kembali data spasial baik dalam bentuk vektor (point/polyline/polygone) maupun dam bentuk raster. Kita akan membutuhkan pustaka "rgdal" untuk data vektor dan pustaka "raster"untuk data raster.

# Memanggil pustaka rgdal dan raster
library(rgdal)
library(raster) 

 # Memanggil data batas kawasan TNGP (polygone)
tngp <- readOGR(dsn = "D:/LUBIS_PRIVATE_DATA/R_GIS_PROJECT/02_Data_mentah/01_Data_vektor", layer="Batas TNGP")

# Memanggil data lokasi penelitian di TNGP (point)
lokasi <- readOGR(dsn = "D:/LUBIS_PRIVATE_DATA/R_GIS_PROJECT/02_Data_mentah/01_Data_vektor", layer="Lokasi penelitian")

# Memanggil data curah hujan tahunan yang diperoleh dari website WORLDCLIME (raster)
annprec <- raster("D:/LUBIS_PRIVATE_DATA/R_GIS_PROJECT/02_Data_mentah/02_Data_Raster/AnnPrep_TNGP.tif") 

Sistem kordinat dan proyeksi dari data spasial yang telah dipanggil tersebut dapat dilihat dengan menggunakan perintah projection seperti berikut;

# Membaca dan memeriksa proyeksi dan sistem kordinat
projection(tngp) # Proyeksi geografi
projection(lokasi) # Proyeksi UTM Zone 48S
projection(annprec) # Proyeksi geografi

Batas kawasan Taman Nasional Gede Pangrango (TNGP) dan data curah hujan tahunan yang telah dipanggil memiliki proyeksi geografis, sementara untuk data lokasi penelitian memiliki proyeksi UTM 48S. Proyeksi geografi tersebut dapat diubah menjadi UTM 48S ataupun sebaliknya. Untuk data yang memiliki format dalam bentuk vektor, kita akan menggunakan perintah spTransform dari pustaka sp, sementara untuk data raster bisa menggunkan perintah projectRaster dari pustaka raster.

# Mengubah proyeksi data vektor TNGP dari Geografis ke UTM
# CARA 1 - Menggunakan proyeksi kordinat yang diperoleh dari data "lokasi"
tngp_utm.1 <- spTransform(tngp, CRS = crs(lokasi))

# CARA 2 - Mengubah proyeksi data vektor Geografi ke UTM
# Menentukan proyeksi tujuan yang diinginkan (UTM)
utm48s <- crs("+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Lalu ubah proyeksinya
tngp_utm.2 <- spTransform(tngp, CRS = crs(utm48s))

# Periksa kembali proyeksi sebelum dan sesudah diproses
projection(tngp) # Geografi
projection(tngp_utm.1) # UTM
projection(tngp_utm.2) # UTM 

Selanjutnya merubah proyeksi data raster dari geografi ke UTM
# Merubah proyeksi data raster menggunakan informasi sebelumnya
annprec_utm <- projectRaster(annprec, crs = utm48s)

# Periksa kembali proyeksi sebelum dan sesudah
projection(annprec) # Geografi
projection(annprec_utm) # UTM

Sebaliknya, kita juga dapat merubah proyeksi UTM ke sistem kordinat geografi, sama halnya dengan sebelumnya.

Kemudian plot peta. Kali ini kita akan menggunakan pola warna yang berbeda dengan menggunakan pustaka RColorBrewer

library(RColorBrewer)
# Menampilkan pilihan warna yang ada pada RColorBrewer
display.brewer.all() # pilih kode warna yang disukai, misalnya "Blues"

# Plot data yang akan menghasilkan Gambar 1 diatas
plot(annprec_utm, col=brewer.pal(30,"Blues"), xlab="Timur (mtr)", ylab="Utara(mtr)")
plot(tngp_utm.1, border="red", lwd=2, add=TRUE)
points(lokasi, pch=18, col = "black")

Data yang telah ditransformasi dapat disimpan kembali kedalam file. Untuk data vektor dapat disimpan menggunakan perintah writeOGR, sementara untuk data raster bisa menggunakan perintah writeRaster seperti contoh dibawah.

# Menyimpan data vektor dengan format shapefile
writeOGR(tngp_utm.1, dsn = "D:/LUBIS_PRIVATE_DATA/R_GIS_PROJECT/04_Hasil",
layer="Batas_TNGP_utm",overwrite_layer = TRUE, driver="ESRI Shapefile")

# Menyimpan data raster dengan format TIF
writeRaster(annprec_utm,overwrite=TRUE, filename="annprec_utm.tif")


Selamat mencoba.....