P3c. Working with Spatial Data

Spatial Data Structure

Tipe data spasial yang paling umum digunakan adalah shapefile, adapun tipe lain yang juga cukup populer adalah KML (Keyhole Markup Language). Data spasial tidak hanya berisi baris dan kolom, namun objek geometrik, seperti titik, garis, ataupun poligon. Data shapefile sebenarnya terdiri dari beberapa file dengan beberapa extension, di antaranya adalah .shp, .shx, dan .dbf. Beberapa package yang umum digunakan untuk bekerja dengan data spasial adalah sp dan rgdal.

Reading Spatial Data

Sebagai ilustrasi, akan digunakan data yang tersedia pada laman http://rtwilson.com/downloads/SnowGIS_SHP.zip. Pastikan Anda mengekstrak folder data tersebut pada direktori yang Anda inginkan. Selanjutnya, package rgdal akan digunakan untuk membaca data SnowGIS tersebut.

install.packages("rgdal")

library(rgdal)

## Loading required package: sp

## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/BIRO PKH - BIG/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/BIRO PKH - BIG/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/BIRO PKH - BIG/Documents/R/win-library/4.0/rgdal/proj

Untuk melihat file apa saja yang ada di dalam folder shapefile tersebut, kita dapat menggunakan fungsi list.files() dan tuliskan direktori Anda masing-masing, ini diberi nama sebagai dsn.

dsn<-paste("D:/GDrive/02-Kuliah/STA581-Sains Data/Bahan/SnowGIS/SnowGIS_SHP")
list.files(dsn)

##  [1] "Cholera_Deaths.dbf"          "Cholera_Deaths.prj"         
##  [3] "Cholera_Deaths.sbn"          "Cholera_Deaths.sbx"         
##  [5] "Cholera_Deaths.shp"          "Cholera_Deaths.shx"         
##  [7] "OSMap.tfw"                   "OSMap.tif"                  
##  [9] "OSMap_Grayscale.tfw"         "OSMap_Grayscale.tif"        
## [11] "OSMap_Grayscale.tif.aux.xml" "OSMap_Grayscale.tif.ovr"    
## [13] "Pumps.dbf"                   "Pumps.prj"                  
## [15] "Pumps.sbx"                   "Pumps.shp"                  
## [17] "Pumps.shx"                   "README.txt"                 
## [19] "SnowMap.tfw"                 "SnowMap.tif"                
## [21] "SnowMap.tif.aux.xml"         "SnowMap.tif.ovr"

Terlihat pada output di atas bahwa folder tersebut memuat beberapa shapefile, di antaranya terdapat 6 file dengan nama Cholera_Deaths dan 5 file bernama Pumps .

Kedua set data tersebut dikenal sebagai layer.

ogrListLayers(dsn)

## [1] "Cholera_Deaths" "Pumps"         
## attr(,"driver")
## [1] "ESRI Shapefile"
## attr(,"nlayers")
## [1] 2

Kita dapat menggunakan fungsi ogrInfo() untuk mengetahui informasi mengenai layer tersebut.

ogrInfo(dsn, layer = "Cholera_Deaths")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum OSGB_1936 in Proj4 definition: +proj=tmerc +lat_0=49
## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs

## Source: "D:\GDrive\02-Kuliah\STA581-Sains Data\Bahan\SnowGIS\SnowGIS_SHP", layer: "Cholera_Deaths"
## Driver: ESRI Shapefile; number of rows: 250 
## Feature type: wkbPoint with 2 dimensions
## Extent: (529160.3 180857.9) - (529655.9 181306.2)
## CRS: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
## LDID: 87 
## Number of fields: 2 
##    name type length typeName
## 1    Id    0      6  Integer
## 2 Count    0      4  Integer

Importing Spatial Data

Fungsi readOGR() dapat digunakan untuk membaca data shapefile.

CholeraDeaths <- readOGR(dsn, layer = "Cholera_Deaths")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum OSGB_1936 in Proj4 definition: +proj=tmerc +lat_0=49
## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GDrive\02-Kuliah\STA581-Sains Data\Bahan\SnowGIS\SnowGIS_SHP", layer: "Cholera_Deaths"
## with 250 features
## It has 2 fields

summary(CholeraDeaths)

## Object of class SpatialPointsDataFrame
## Coordinates:
##                min      max
## coords.x1 529160.3 529655.9
## coords.x2 180857.9 181306.2
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs]
## Number of points: 250
## Data attributes:
##        Id        Count       
##  Min.   :0   Min.   : 1.000  
##  1st Qu.:0   1st Qu.: 1.000  
##  Median :0   Median : 1.000  
##  Mean   :0   Mean   : 1.956  
##  3rd Qu.:0   3rd Qu.: 2.000  
##  Max.   :0   Max.   :15.000

Selanjutnya kita dapat memeriksa class dari data CholeraDeaths tersebut.

class(CholeraDeaths)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Data tersebut merupakan SpatialPointsDataFrame yang termasuk S4 class, maka untuk mengakses data slot perlu digunakan notasi @. Catatan : Mengenai S4 Class ini, akan dipelajari pada Mata Kuliah STA561-Pemrograman Statistika pada pertemuan ke 6.

str(CholeraDeaths@data)

## 'data.frame':    250 obs. of  2 variables:
##  $ Id   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Count: int  3 2 1 1 4 2 2 2 3 2 ...

Making Map

Fungsi plot() dapat digunakan untuk membuat grafik paling sederhana dari data CholeraDeaths.

par(mfrow=c(1,2))
plot(CholeraDeaths)
plot(CholeraDeaths, pch=20, col="steelblue")

Perhatikan bahwa plot di atas hanya menunjukkan sebaran titik spasial, tanpa memberikan informasi yang jelas tentang lokasi data tersebut. Jika kita memiliki peta dalam bentuk data polygon, kita dapat mengimpor data tersebut dengan cara yang sama (seandainya datanya berupa shapefile), kemudian kita plot peta baru kemudian plot data titik seperti di atas.

Alternatif lainnya jika kita tidak ingin menggunakan peta polygon dari shapefile, kita dapat menggunakan beberapa package yang tersedia di R software, seperti ggmap, OpenStreetMap, leaflet, atau yang lain. Namun perhatikan bahwa untuk bisa menggunakan package OpenStreetMap, Anda harus memastikan bahwa jika Anda menggunakan R 64-bit maka Java yang terinstall di PC Anda juga harus sesuai, yaitu 64-bit.

Berikut ini akan ditunjukkan salah satu cara menampilkan peta dengan memanfaatkan package leaflet.

install.packages("leaflet")

library(leaflet)
map <- leaflet() %>% setView(lng =  -0.13659, lat =51.51328 , zoom = 12)
map %>% addTiles()

Sebelum kedua peta dan data titik digabungkan. Pastikan terlebih dahulu apakah koordinat yang digunakan menggunakan skala yang sama.

head(coordinates(CholeraDeaths))

##      coords.x1 coords.x2
## [1,]  529308.7  181031.4
## [2,]  529312.2  181025.2
## [3,]  529314.4  181020.3
## [4,]  529317.4  181014.3
## [5,]  529320.7  181007.9
## [6,]  529336.7  181006.0

Seperti terlihat di atas, koordinat pada data CholeraDeaths diukur pada skala yang berbeda dengan peta yang diambil dari package leaflet. Terdapat beberapa macam coordinate reference system (CRS), beberapa di antaranya yang cukup populer adalah suatu set EPSG (European Petroleum Survey Group) berikut:

EPSG:4326 juga dikenal sebagai WGS84, ukuran standard yang digunakan pada sistem GPS dan Google Earth.
EPSG:3857 digunakan pada Google Maps, Open Street Maps, dsb.
EPSG:27700 juga dikenal sebagai OSGB 1936, atau British National Grid: United Kingdom Ordnance Survey.

Berikut ini dilakukan transformasi menjasi WGS84.

cholera_latlong <- CholeraDeaths %>% 
  spTransform(CRS("+init=epsg:4326"))

Berikut ini, untuk menggabungkan peta dan titik.

leaflet(data = CholeraDeaths) %>% 
  addTiles() %>%
  addMarkers(cholera_latlong@coords[,1], cholera_latlong@coords[,2])

Dapat dilihat di atas, bahwa setelah koordinatnya disamakan, kita dapat menampilkan data CholeraDeaths pada peta yang diperoleh dari Open Street Map melalui package leaflet.

Referensi

Agafonkin, V. (n.d.). Leaflet for R - Markers. rstudio.github.io. Retrieved from https://rstudio.github.io/leaflet/markers.html

UQ SLC Digital Team. (2020, April 16). Creating maps using R. Language Technology and Data Analysis Laboratory (LADAL). Retrieved from https://slcladal.github.io/maps.html

IPB University, rahmaanisa@apps.ipb.ac.id ↩︎
IPB University, gerry_dito@apps.ipb.ac.id ↩︎
Badan Informasi Geospasial, abdul.aziz@big.go.id ↩︎