Installing R and Rstudio

If you do not have R and RStudio on your computer, proceed as follows:

Download base R for your operating system from https://cran.r-project.org.
Install it on your system.
Download RStudio desktop version for your operating system from https://www.rstudio.com/products/RStudio/.
Install it on your system.

Installing Packages

Jalankan R

Jika terhubung ke internet:
install.packages(“namapackage”)

Contoh:

install.packages(c("sp", "gstat"))

Jika sudah download binary (*.zip):
install.packages(“drive:/namafile.zip”,repos=NULL)

Contoh:

install.packages(“C:/Program Files/R/nama file.zip”,repos=NULL)

atau melalui menu: Packages > Install packages(s) from local files…

First Interaction with R

Task 1: Compute the number of radians in one circular degree.

In this document, input and output are shown as follows:

2 * pi/360

## [1] 0.01745329

This means: when R prompts for input:

You type an expression, in this case 2*pi/360, and press the Enter
R then shows the result in the console. In this case it is a one-element vector (shown by the [1]), with the value 0.017453

Note: that pi is constant value known to R.

Task 2: Load the `sp` and `gstat` packages into the workspace.

You can load these in RStudio by checking the small checkbox next to the packages name in the Packages tab; see Figure 3.

You can also load these this from the R console with the library function. This loads a packages into the R workspace:

library(sp)
library(gstat)

Loading Meuse Dataset

Task 3: Load the `meuse` dataset into the workspace.

The data function loads a dataset. We show the contents of the workspace before and after with the ls() function:

ls()

## character(0)

data(meuse)
ls()

## [1] "meuse"

What Objects Were in the workspace before and after loading the meuse dataset?

Task 4: Examine the structure of the meuse dataset

The str() function shows the structure of an R object:

str(meuse)

## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...

Description of `meuse` dataset

help(meuse)

Meuse river data set

Description

This data set gives locations and topsoil heavy metal concentrations, along with a number of soil and landscape variables at the observation locations, collected in a flood plain of the river Meuse, near the village of Stein (NL). Heavy metal concentrations are from composite samples of an area of approximately 15 m x 15 m.

Usage

data(meuse)

Format

This data frame contains the following columns:

x
a numeric vector; Easting (m) in Rijksdriehoek (RDH) (Netherlands topographical) map coordinates

y
a numeric vector; Northing (m) in RDH coordinates

Introduction to Spatial Data Type

Types of spatial Data

Plot the data

Points

coordinates(meuse) <- c("x", "y")
plot(meuse)
title("points")

Polygons

data(meuse.riv)
meuse.1st <- list(Polygons(list(Polygon(meuse.riv)),"meuse.riv"))
meuse.sr <- SpatialPolygons(meuse.1st)
plot(meuse.sr, col = "grey")
title("polygons")

Grid

data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixels")
image(meuse.grid, col = "grey")
title("grid")

image(meuse.grid, col = "lightgrey")
plot(meuse.sr, col = "grey", add = TRUE)
plot(meuse,add = TRUE)

Data Raster

f <- system.file("external/test.grd", package="raster")
library(raster)
r <- raster(f)
plot(r)

Menggunakan Data Eksternal

Data pertumbuhan penduduk

Impor data:

datapop <- read.csv('http://bit.ly/Popgrowth2000', header=T, sep=',')

Dimensi data:

dim(datapop)

## [1] 293   5

Struktur data:

str(datapop)

## 'data.frame':    293 obs. of  5 variables:
##  $ City          : chr  "Abidjan" "Adana" "Addis Ababa" "Adelaide" ...
##  $ Country       : chr  "Cote d'Ivoire" "Turkey" "Ethiopia" "Australia" ...
##  $ Latitude      : num  5.31 37 9.02 -34.93 23.04 ...
##  $ Longitude     : num  -4.01 35.33 38.75 138.6 72.57 ...
##  $ PopGrowth_2000: int  1929000 1041000 2424000 1092000 2954000 1582000 3339000 2562000 1135000 1147000 ...

Tabel data:

View(datapop)

City	Country	Latitude	Longitude	PopGrowth_2000
Abidjan	Cote d’Ivoire	5.31	-4.013	1929000
Adana	Turkey	37	35.33	1041000
Addis Ababa	Ethiopia	9.025	38.75	2424000
Adelaide	Australia	-34.93	138.6	1092000
Ahmadabad	India	23.04	72.57	2954000
Aleppo	Syria	36.22	37.17	1582000
Alexandria	Egypt	31.2	29.92	3339000
Algiers	Algeria	36.75	3.042	2562000
Almaty	Kazakhstan	43.26	76.91	1135000
Amman	Jordan	31.96	35.95	1147000

Plot the Data points

coordinates(datapop) <- c("Longitude","Latitude")
plot(datapop)

size<-datapop$PopGrowth_2000/sum(datapop$PopGrowth_2000)
plot(datapop,pch=20, col="steelblue", cex=size*100)

Menambahkan Peta Dunia

Install Packages:

install.packages("rworldmap")
library(rworldmap)
data(package="rworldmap")
data(countriesCoarse,envir=environment(),package="rworldmap")

Struktur data:

str(countriesCoarse)

Peta dunia:

plot(countriesCoarse)

## Warning in wkt(obj): CRS object has no comment

Peta dunia + Data pertumbuhan penduduk:

plot(datapop,add=T, pch=20)

## Warning in wkt(obj): CRS object has no comment

Coba tampilkan data pertumbuhan penduduk dengan ukuran titik sesuai dengan besarnya pertumbuhan, pada peta dunia.!

Rasterize:

library(raster)
r <- raster(datapop)
res(r)<-c(5,5)
nc <- rasterize(coordinates(datapop), r, fun=mean, background=NA)
plot(nc)
plot(countriesCoarse, add=TRUE)

Cara lain menampilkan peta

install.packages("ggmap")
library(ggmap)
map <- get_map(location = 'Indonesia', zoom = 4)
plot(map)
datapop<-read.csv('http://bit.ly/Popgrowth2000', header=T,
sep=',')
ggmap(map)+geom_point(aes(x = Longitude, y = Latitude),
data = pop_ina,col="red")

Cara lain menampilkan peta

Ilustrasi di bawah ini dirujuk dari Data Technik (2019). Data tersedia pada package datasets, berisi 50 observasi dan 9 peubah.

states<-as.data.frame(state.x77)
head(states)

##            Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
## Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
## Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
## Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
## Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
## California      21198   5114        1.1    71.71   10.3    62.6    20 156361
## Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766

Pertama, kita akan coba menampilkan visualisasi data populasi di U.S menggunakan fungsi di dalam package ggplot2.

library(ggplot2)
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.0     v forcats 0.5.1
## v purrr   0.3.4

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()

states$region <- tolower(rownames(states))
states_map <- map_data("state")
fact_join <- left_join(states_map, states, by = "region")
ggplot(fact_join, aes(long, lat, group = group))+
  geom_polygon(aes(fill = Population), color = "white")+
  scale_fill_viridis_c(option = "C")+
  theme_classic()

Kita dapat pula mengganti pewarnaan grafik seperti di bawah ini.

ggplot(fact_join, aes(long, lat, group = group))+
  geom_polygon(aes(fill = Income), color = "white")+
  scale_fill_viridis_c(option = "C")+
  theme_classic()

Berikut apabila ingin menampilkan peubah lain, misalnya usia harapan hidup.

fact_join$`Life Exp` <- as.numeric(fact_join$`Life Exp`)
ggplot(fact_join, aes(long, lat, group = group))+
geom_polygon(aes(fill = `Life Exp`), color = "white")+
scale_fill_viridis_c(option = "C")+
theme_classic()

Ilustrasi: Data Cholera

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.

library(rgdal)

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 dikenal sebagai dsn.

dsn<-paste("D:/Dept.STK/Courses/Statistika Spasial (S2)/TA 2020-2021/Semester Genap/P1 - Genap 2021/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_Deathsdan 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:\Dept.STK\Courses\Statistika Spasial (S2)\TA 2020-2021\Semester Genap\P1 - Genap 2021\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

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:\Dept.STK\Courses\Statistika Spasial (S2)\TA 2020-2021\Semester Genap\P1 - Genap 2021\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 @.

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 ...

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.

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.

cholera_latlong <- CholeraDeaths %>% 
  spTransform(CRS("+init=epsg:4326"))
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.

Tugas Mandiri: Peta Pulau Jawa

1. Load package untuk mengimpor data polygon (shape file):

install.packages('rgdal')
library(rgdal)

2. Mengimpor data dari direktori

Data tersedia di http://bit.ly/ShapeFile_Jawa)

Menggunakan fungsi readOGR: readOGR(dsn=‘folder (direktori) tempat shape file disimpan’, layer=‘nama shape file (tanpa .shp extension)’)

Contoh:

jawa <- readOGR(dsn='D:/Belajar/Spasial/R/Map of Jawa (original)', layer='jawa')

Pratinjau data:

head(jawa@data)

##   MISKIN KODE_KAB    NAMA_KAB KODE_PROP  NAMA_PROP
## 0      0     3501     Pacitan        35 Jawa Timur
## 1      0     3502    Ponorogo        35 Jawa Timur
## 2      0     3503  Trenggalek        35 Jawa Timur
## 3      0     3504 Tulungagung        35 Jawa Timur
## 4      0     3508    Lumajang        35 Jawa Timur
## 5      0     3511   Bondowoso        35 Jawa Timur

3. Menampilkan peta Pulau Jawa

plot(jawa)

4. Menampilkan nama kota/kabupaten

text(jawa,'NAMA_KAB',cex=0.5)

5. Memberi warna yang berbeda utk setiap Propinsi

plot(jawa,col=jawa$KODE_PROP-30)

Mengganti referensi warna:

palette(rainbow(6))
plot(jawa, col=jawa$KODE_PROP-30)

palette(terrain.colors(6))
plot(jawa, col=jawa$KODE_PROP-30)

Pengenalan Struktur R untuk Data Spasial

rahmaanisa@apps.ipb.ac.id

Praktikum 1 | Statistika Spasial

January 28, 2022

Installing R and Rstudio

Installing Packages

First Interaction with R

Task 1: Compute the number of radians in one circular degree.

Task 2: Load the `sp` and `gstat` packages into the workspace.

Loading Meuse Dataset

Task 3: Load the `meuse` dataset into the workspace.

Task 4: Examine the structure of the meuse dataset

Description of `meuse` dataset

Meuse river data set

Introduction to Spatial Data Type

Types of spatial Data

Plot the data

Points

Polygons

Grid

Data Raster

Menggunakan Data Eksternal

Data pertumbuhan penduduk

Plot the Data points

Menambahkan Peta Dunia

Cara lain menampilkan peta

Ilustrasi: Data Cholera

Tugas Mandiri: Peta Pulau Jawa

1. Load package untuk mengimpor data polygon (shape file):

2. Mengimpor data dari direktori

3. Menampilkan peta Pulau Jawa

4. Menampilkan nama kota/kabupaten

5. Memberi warna yang berbeda utk setiap Propinsi

Pengenalan Struktur R untuk Data Spasial

rahmaanisa@apps.ipb.ac.id

Praktikum 1 | Statistika Spasial

January 28, 2022

Installing R and Rstudio

Installing Packages

First Interaction with R

Task 1: Compute the number of radians in one circular degree.

Task 2: Load the sp and gstat packages into the workspace.

Loading Meuse Dataset

Task 3: Load the meuse dataset into the workspace.

Task 4: Examine the structure of the meuse dataset

Description of meuse dataset

Meuse river data set

Introduction to Spatial Data Type

Types of spatial Data

Plot the data

Points

Polygons

Grid

Data Raster

Menggunakan Data Eksternal

Data pertumbuhan penduduk

Plot the Data points

Menambahkan Peta Dunia

Cara lain menampilkan peta

Ilustrasi: Data Cholera

Tugas Mandiri: Peta Pulau Jawa

1. Load package untuk mengimpor data polygon (shape file):

2. Mengimpor data dari direktori

3. Menampilkan peta Pulau Jawa

4. Menampilkan nama kota/kabupaten

5. Memberi warna yang berbeda utk setiap Propinsi

Task 2: Load the `sp` and `gstat` packages into the workspace.

Task 3: Load the `meuse` dataset into the workspace.

Description of `meuse` dataset