Pengenalan R untuk Data Spasial

Packages

# Install packages sp dan gstat kemudian load
library(sp) 
library(gstat)
## Warning: package 'gstat' was built under R version 4.0.5

Data : Meuse Dataset

data(meuse)
ls()
## [1] "meuse"
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 ...
# Plot data
# poins
help(meuse)
## starting httpd help server ... done
coordinates(meuse)<- ~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")

# gabung grid(daerah) dengan points(titik) dan tepinya polygon
image(meuse.grid, col="lightgrey") #grid (daerahnya)
plot(meuse.sr, col="grey", add=TRUE) #polygon di tepinya
plot(meuse, add=TRUE)

Menggunakan Data Eksternal

Data pertumbuhan penduduk

# import 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 ...
# lihat data dalam tabel/dataframe
View(datapop)

# Plot the Data points
coordinates(datapop) <- c("Longitude", "Latitude")
plot(datapop)

#memberi warna dan ketebalan titiknya
size <- datapop$PopGrowth_2000/sum(datapop$PopGrowth_2000)
plot(datapop, pch=20, col="steelblue", cex=size*100)

# Menambah Peta Dunia
# install packages and load
library(rworldmap)
## Warning: package 'rworldmap' was built under R version 4.0.3
## ### Welcome to rworldmap ###
## For a short introduction type :   vignette('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)

# Rasterize
library(raster)
## Warning: package 'raster' was built under R version 4.0.4
r <- raster(datapop)
nc <- rasterize(coordinates(datapop), r, fun=mean, background=NA)
plot(nc)
plot(countriesCoarse, add=TRUE)

Peta Pulau Jawa

1. Load package

library(rgdal)
## Warning: package 'rgdal' was built under R version 4.0.4
## 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/ASUS/OneDrive/Dokumen/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/ASUS/OneDrive/Dokumen/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/ASUS/OneDrive/Dokumen/R/win-library/4.0/rgdal/proj

2. Impor data dari direktori

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

jawa <- readOGR(dsn='C:/Users/ASUS/Downloads/Map of Jawa (original)', layer='jawa')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum D_unknown in Proj4 definition: +proj=longlat
## +ellps=GRS80 +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\ASUS\Downloads\Map of Jawa (original)", layer: "jawa"
## with 116 features
## It has 5 fields
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)
text(jawa, 'NAMA_KAB', cex=0.5)

4. Menampilkan nama kota/kabupaten

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

5. Memberi warna yang berbeda untuk setiap provinsi

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)


  1. Statistika dan Sains Data [IPB University], ↩︎