Sub project: 3.5 Plot Data dan Layer OpenStreetMap

File: Project_03_5.R

Tujuan: Menampilkan data spasial dan layer dari OpenStreetMap.

Yang dilakukan pertama adalah mempersiapkan lingkungan sistem dan mendefinisikan folder kerja.

rm(list=ls())
graphics.off()
setwd("D:/_randmm/project03")
getwd()
## [1] "D:/_randmm/project03"

Berikutnya adalah mengunggah paket GISTools, OpenStreetMap, dan data ke dalam sistem. Data spasial yang digunakan adalah batas kecamatan Kabupaten Karawang sebagaimana data sebelum pekerjaan ini. Saat mengunggah data langsung kita berikan referensi koordinat.

# Unggah paket GISTools dan OpenStreetMap
library(GISTools)
## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: TRUE
## Loading required package: RColorBrewer
## Loading required package: MASS
## Loading required package: rgeos
## rgeos version: 0.3-19, (SVN revision 524)
##  GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084 
##  Linking to sp version: 1.2-3 
##  Polygon checking: TRUE
library(OpenStreetMap)

# Unggah data kedalam object
KW_kec <-readShapePoly("karawang_admkec", proj4string = CRS("+init=EPSG:32748"))

Lakukan transformasi dari UTM ke Geographic (lat/lon). Buat obyek spasial baru untuk outline.

# Transform proyeksi UTM 48S ke Geo Lat/Lon
KW_kec.geo <- spTransform(KW_kec, CRS("+proj=longlat +datum=WGS84"))
KW_kec.outline <- gUnaryUnion(KW_kec.geo, id=NULL)

Definisikan koordinat untuk batas kiri atas dan kanan bawah.

# Definisikan batas kiri atas dan kanan bawah
ul <- as.vector(cbind(bbox(KW_kec.geo)[2,2], bbox(KW_kec.geo)[1,1]))
lr <- as.vector(cbind(bbox(KW_kec.geo)[2,1], bbox(KW_kec.geo)[1,2]))

Mengunduh peta melalui OpenStreetMap dari MapQuest. Untuk mengunduh ini kita gunakan fungsi openmap(). Lama pengunduhan tergantung pada koneksi internet. Setelah itu lalu siapkan window plot kemudian lakukan plotting layer dari MapQuest, layer batas kecamatan, dan layer outline.

# Unduh peta dari MapQuest
MyMapQuest <- openmap(ul, lr, NULL, 'mapquest')

# Definisikan window plot dan tampilkan peta
par(mar= c(0,0,0,0))
plot(MyMapQuest, removeMargin=FALSE)
plot(spTransform(KW_kec.geo, osm()), add=TRUE, lwd=2, border="blue")
plot(spTransform(KW_kec.outline, osm()), add=TRUE, lwd=2)