Eksplorasi Data Spasial

Author

Alfa Nugraha

Manipulasi Data Vektor

Pada bagian ini ditunjukkan beberapa cara dalam memanipulasi data vector

depok_vect <- vect("kota_depok_2020.shp")
depok_vect
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 8  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  source      : kota_depok_2020.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area    ADM3_EN ADM3_PCODE    ADM2_EN
##  type        :      <num>      <num>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207       Beji  ID3276050 Kota Depok
##                    0.2806   0.001583 Bojongsari  ID3276011 Kota Depok
##                    0.2442   0.001271   Cilodong  ID3276031 Kota Depok
##  ADM2_PCODE    ADM1_EN ADM1_PCODE
##       <chr>      <chr>      <chr>
##      ID3276 Jawa Barat       ID32
##      ID3276 Jawa Barat       ID32
##      ID3276 Jawa Barat       ID32

plot(depok_vect, "ADM3_EN")

Dasar

Atribut dan geometri

Untuk menampilkan informasi atribut dari SpatVector, gunakan:

depok_df <- as.data.frame(depok_vect)
depok_df

Shape_Leng	Shape_Area	ADM3_EN	ADM3_PCODE	ADM2_EN	ADM2_PCODE	ADM1_EN	ADM1_PCODE
0.1950645	0.0012065	Beji	ID3276050	Kota Depok	ID3276	Jawa Barat	ID32
0.2805737	0.0015832	Bojongsari	ID3276011	Kota Depok	ID3276	Jawa Barat	ID32
0.2442310	0.0012712	Cilodong	ID3276031	Kota Depok	ID3276	Jawa Barat	ID32
0.3565971	0.0017641	Cimanggis	ID3276040	Kota Depok	ID3276	Jawa Barat	ID32
0.1552971	0.0008562	Cinere	ID3276061	Kota Depok	ID3276	Jawa Barat	ID32
0.1708723	0.0008929	Cipayung	ID3276021	Kota Depok	ID3276	Jawa Barat	ID32

Ekstrak geometri sebagai matriks juga dapat dilakukan (tapi jarang sekali diperlukan)

depok_mtx <- geom(depok_vect)
head(depok_mtx)
##      geom part        x         y hole
## [1,]    1    1 106.8353 -6.356507    0
## [2,]    1    1 106.8352 -6.357034    0
## [3,]    1    1 106.8352 -6.357333    0
## [4,]    1    1 106.8352 -6.357709    0
## [5,]    1    1 106.8352 -6.358073    0
## [6,]    1    1 106.8352 -6.358186    0

Well-Known Text (WKT)

Skema penulisan bentuk representasi dari suatu geometri secara sederhana diatur ke dalam sebuah kumpulan string menjadi format teks baku yang dikenal dengan Well-Known Text (WKT)¹.

¹ Lott, Roger. 2015. “Geographic Information-Well-Known Text Representation of Coordinate Reference Systems.” Open Geospatial Consortium. http://docs.opengeospatial.org/is/12-063r5/12-063r5.html

Bentuk geometri yang umum digunakan ditulis sebagai berikut:


tipe	deskripsi
`POINT`	single point geometry
`MULTIPOINT`	set of points
`LINESTRING`	single linestring (two or more points connected by straight lines)
`MULTILINESTRING`	set of linestrings
`POLYGON`	exterior ring with zero or more inner rings, denoting holes
`MULTIPOLYGON`	set of polygons
`GEOMETRYCOLLECTION`	set of the geometries above

Sebagai contoh berikut ini adalah bentuk 7 geometri yang umum diketahui:

Sehingga WKT dalam bentuk geometri pada gambar di atas:

POINT (0 1)
MULTIPOINT ((1 1), (2 2), (4 1), (2 3), (1 4))
LINESTRING (1 1, 5 5, 5 6, 4 6, 3 4, 2 3)
MULTILINESTRING ((1 1, 5 5, 5 6, 4 6, 3 4, 2 3), (3 0, 4 1, 2 1))
POLYGON ((2 1, 3 1, 5 2, 6 3, 5 3, 4 4, 3 4, 1 3, 2 1),
    (2 2, 3 3, 4 3, 4 2, 2 2))
MULTIPOLYGON (((2 1, 3 1, 5 2, 6 3, 5 3, 4 4, 3 4, 1 3, 2 1),
    (2 2, 3 3, 4 3, 4 2, 2 2)), ((3 7, 4 7, 5 8, 3 9, 2 8, 3 7)))
GEOMETRYCOLLECTION (
    POLYGON ((2 1, 3 1, 5 2, 6 3, 5 3, 4 4, 3 4, 1 3, 2 1),
      (2 2 , 3 3, 4 3, 4 2, 2 2)),
    LINESTRING (1 6, 5 10, 5 11, 4 11, 3 9, 2 8),
    POINT (2 5),
    POINT (5 4)
)

Data tersebut dapat dilihat dengan sintaks berikut

depok_wkt <- geom(depok_vect, wkt=TRUE)
substr(depok_wkt, 1, 50)
##  [1] "POLYGON ((106.835259 -6.356507, 106.835175 -6.3570"
##  [2] "POLYGON ((106.730734 -6.3605, 106.73082 -6.360508,"
##  [3] "POLYGON ((106.859915 -6.398481, 106.86 -6.39849, 1"
##  [4] "POLYGON ((106.850464 -6.336821, 106.850632 -6.3368"
##  [5] "POLYGON ((106.807165 -6.314414, 106.80893 -6.31451"
##  [6] "POLYGON ((106.809376 -6.413101, 106.809749 -6.4137"
##  [7] "POLYGON ((106.79044 -6.3476, 106.79045 -6.34844, 1"
##  [8] "POLYGON ((106.800798 -6.383479, 106.800775 -6.3837"
##  [9] "POLYGON ((106.769799 -6.356591, 106.770083 -6.3567"
## [10] "POLYGON ((106.86012 -6.372697, 106.86038 -6.372718"
## [11] "POLYGON ((106.88813 -6.38563, 106.888646 -6.38578,"

Variabel

Ekstrak variabel tertentu dapat dilakukan seperti pada data frame

depok_vect$ADM3_EN
##  [1] "Beji"         "Bojongsari"   "Cilodong"     "Cimanggis"    "Cinere"      
##  [6] "Cipayung"     "Limo"         "Pancoran Mas" "Sawangan"     "Sukma Jaya"  
## [11] "Tapos"

Berbeda dengan cara di atas, ekstrak satu atau lebih variabel juga dapat dilakukan beserta bentuk geometrinya

depok_vect[, "ADM3_EN"]
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 1  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  source      : kota_depok_2020.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :    ADM3_EN
##  type        :      <chr>
##  values      :       Beji
##                Bojongsari
##                  Cilodong

Variabel baru juga dapat ditambahkan ke SpatVector seperti halnya pada data frame

depok_vect$lets <- sample(letters, nrow(depok_vect))
depok_vect
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 9  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  source      : kota_depok_2020.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area    ADM3_EN ADM3_PCODE    ADM2_EN
##  type        :      <num>      <num>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207       Beji  ID3276050 Kota Depok
##                    0.2806   0.001583 Bojongsari  ID3276011 Kota Depok
##                    0.2442   0.001271   Cilodong  ID3276031 Kota Depok
##  ADM2_PCODE    ADM1_EN ADM1_PCODE  lets
##       <chr>      <chr>      <chr> <chr>
##      ID3276 Jawa Barat       ID32     o
##      ID3276 Jawa Barat       ID32     d
##      ID3276 Jawa Barat       ID32     p

Menghapus variabel

p$lets <- NULL

Merge

Menambahkan atribut baru berdasarkan kolom unik dapat dilakukan dengan perintah merge

df_val <- data.frame(Kecamatan = depok_vect$ADM3_EN, 
                     Value = round(runif(length(depok_vect), 100, 1000)))

depok_merge <- merge(depok_vect, df_val, by.x="ADM3_EN", by.y="Kecamatan") 
depok_merge
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 10  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :    ADM3_EN Shape_Leng Shape_Area ADM3_PCODE    ADM2_EN
##  type        :      <chr>      <num>      <num>      <chr>      <chr>
##  values      :       Beji     0.1951   0.001207  ID3276050 Kota Depok
##                Bojongsari     0.2806   0.001583  ID3276011 Kota Depok
##                  Cilodong     0.2442   0.001271  ID3276031 Kota Depok
##  ADM2_PCODE    ADM1_EN ADM1_PCODE  lets Value
##       <chr>      <chr>      <chr> <chr> <num>
##      ID3276 Jawa Barat       ID32     o   123
##      ID3276 Jawa Barat       ID32     d   938
##      ID3276 Jawa Barat       ID32     p   275

Baris/record

Memilih baris/record tertentu pada SpatVector

beji <- depok_vect[which(depok_vect$ADM3_EN == 'Beji'), ]
beji
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 9  (geometries, attributes)
##  extent      : 106.7943, 106.8419, -6.390085, -6.354528  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area ADM3_EN ADM3_PCODE    ADM2_EN ADM2_PCODE
##  type        :      <num>      <num>   <chr>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207    Beji  ID3276050 Kota Depok     ID3276
##     ADM1_EN ADM1_PCODE  lets
##       <chr>      <chr> <chr>
##  Jawa Barat       ID32     o

Append

depok_rast <- rast('kota_depok.tif')
dim(depok_rast) <- c(2, 2)
values(depok_rast) <- 1:4
names(depok_rast) <- 'zona'

depok_zone <- as.polygons(depok_rast)
depok_zone
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 4, 1  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  zona
##  type        : <int>
##  values      :     1
##                    2
##                    3

depok_zone2 <- depok_zone[2,]
plot(depok_vect)
plot(depok_zone, add=TRUE, border='blue', lwd=5)
plot(depok_zone2, add=TRUE, border='red', lwd=2, col='red')

depok_append <- rbind(depok_vect, depok_zone)
depok_append
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 15, 10  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  source      : kota_depok_2020.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area    ADM3_EN ADM3_PCODE    ADM2_EN
##  type        :      <num>      <num>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207       Beji  ID3276050 Kota Depok
##                    0.2806   0.001583 Bojongsari  ID3276011 Kota Depok
##                    0.2442   0.001271   Cilodong  ID3276031 Kota Depok
##  ADM2_PCODE    ADM1_EN ADM1_PCODE  lets  zona
##       <chr>      <chr>      <chr> <chr> <int>
##      ID3276 Jawa Barat       ID32     o    NA
##      ID3276 Jawa Barat       ID32     d    NA
##      ID3276 Jawa Barat       ID32     p    NA

Aggregate

Agregasi juga dapat dilakukan pada data spasial yang dilakukan pada atribut dengan nilai yang sama, teknik ini dikenal pula dengan istilah dissolve

depok_aggr <- aggregate(depok_vect, by='ADM2_EN')
plot(depok_aggr)

Jika ditambahkan dengan depok_zone pada data sebelumnya

dz <- aggregate(depok_zone)
plot(dz, col='light gray', border='light gray', lwd=5)
plot(depok_aggr, add=TRUE, col=rainbow(3), lwd=3, border='white')

Overlay

Hapus

depok_erase <- erase(depok_vect, depok_zone2)
plot(depok_erase)

Intersect

depok_intersect <- intersect(depok_vect, depok_zone2)
plot(depok_intersect)

Crop

# xmin, xmax, ymin, ymax
crop_box <- ext(106.77, 106.82, -6.46, -6.40) 
depok_crop <- crop(depok_vect, crop_box)
plot(depok_crop)

plot(depok_vect)
plot(crop_box, add=TRUE, lwd=3, col="red")
plot(depok_crop, col='light blue', add=TRUE)
plot(crop_box, add=TRUE, lwd=3, border="blue")

Union

depok_union <- union(depok_vect, depok_zone)
depok_union
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 27, 10  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area ADM3_EN ADM3_PCODE    ADM2_EN ADM2_PCODE
##  type        :      <num>      <num>   <chr>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207    Beji  ID3276050 Kota Depok     ID3276
##                    0.1951   0.001207    Beji  ID3276050 Kota Depok     ID3276
##                    0.1951   0.001207    Beji  ID3276050 Kota Depok     ID3276
##     ADM1_EN ADM1_PCODE  lets  zona
##       <chr>      <chr> <chr> <int>
##  Jawa Barat       ID32     o     1
##  Jawa Barat       ID32     o     2
##  Jawa Barat       ID32     o     3

Dapat dilihat bahwa terdapat banyak poligon yang berisi semua kombinasi poligon dan atribut

plot(depok_union, col=sample(rainbow(length(depok_union))))

Cover

cover adalah kombinasi dari intersect dan union. intersect mengembalikan geometri baru dengan atribut yang telah terpotong dari dua gugus data. union menggabungkan kombinasi kedua data. cover mengembalikan hasil pemotongan dan penggabungan geometri dan atribut dari dua gugus data

depok_cover <- cover(depok_vect, depok_zone[c(1, 4), ])
depok_cover
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 11, 10  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       : Shape_Leng Shape_Area    ADM3_EN ADM3_PCODE    ADM2_EN
##  type        :      <num>      <num>      <chr>      <chr>      <chr>
##  values      :     0.1951   0.001207       Beji  ID3276050 Kota Depok
##                    0.2806   0.001583 Bojongsari  ID3276011 Kota Depok
##                    0.2442   0.001271   Cilodong  ID3276031 Kota Depok
##  ADM2_PCODE    ADM1_EN ADM1_PCODE  lets  zona
##       <chr>      <chr>      <chr> <chr> <int>
##      ID3276 Jawa Barat       ID32     o    NA
##      ID3276 Jawa Barat       ID32     d    NA
##      ID3276 Jawa Barat       ID32     p    NA

plot(depok_cover)

Difference

Perbedaan simetris dari dua SpatVector

depok_dif <- symdif(depok_zone, depok_vect)
depok_dif
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 4, 1  (geometries, attributes)
##  extent      : 106.7167, 106.9188, -6.461283, -6.314414  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  zona
##  type        : <int>
##  values      :     1
##                    2
##                    3

plot(depok_dif, col=rainbow(length(depok_dif)))

Spatial query

Berbeda dengan query pada umumnya, query pada data spasial dalam hal ini poligon dapat dilakukan dengan titik (point-in-polygon query)

# menentukan titik
pts <- matrix(c(106.75, 106.80, 106.83, 106.90, -6.43, -6.38, -6.4, -6.41), ncol=2)
# membuat data vector titik
spts <- vect(pts, crs=crs(depok_vect))

plot(depok_zone, col='light blue', lwd=2)
points(spts, col='light gray', pch=20, cex=6)
text(spts, 1:nrow(pts), col='red', font=2, cex=1.5)
lines(depok_vect, col='blue', lwd=2)

Kemudian untuk mendapatkan informasi pada titik tersebut, gunakan

extract(spts, depok_vect)
##      id.y id.x
## [1,]    1   NA

extract(spts, depok_zone)
##      id.y id.x
## [1,]    1   NA

Manipulasi Data Raster

Terdapat beberapa format data raster yang umum digunakan sebagai bagian dari data penginderaan jauh (remote sensing) yaitu GeoTiff, ESRI, ENVI, dan ERDAS. Semua gugus data raster tersebut berukuran sangat besar jika dimuat ke dalam RAM sehingga penggunaannya dalam bentuk file.

Pada bahasan ini akan melibatkan fungsi-fungsi pengolahan data raster yang ada pada terra, salah satunya adalah raster aljabar. Raster yang digunakan dibuat di dalam fungsi rast yang disimpan ke dalam memory (RAM) secara temporary.

Membuat objek `SpatRaster`

Pengaturan bawaan fungsi rast adalah membuat struktur data raster dengan longitude/latitude dan 1-by-1 degree cell. Pengaturan ini dapat diubah dengan memasukkan sejumlah argumen. Proyeksi bisa diatur pada data raster namun untuk mengubahnya harus menggunakan fungsi projection.

# rast dengan argumen default
rst <- rast()
rst
## class       : SpatRaster 
## dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84

# rast dengan input argumen 
rst <- rast(ncol=36, nrow=18, xmin=-1000, xmax=1000, ymin=-100, ymax=900)

# mengubah resolusi 
res(rst)
## [1] 55.55556 55.55556
res(rst) <- 100
res(rst)
## [1] 100 100

Mengubah jumlah kolom (berefek pada resolusi)

ncol(rst)
## [1] 20
ncol(rst) <- 18
ncol(rst)
## [1] 18
res(rst)
## [1] 111.1111 100.0000

Mengatur CRS (mendefinisikan proyeksi)

crs(rst) <- "+proj=utm +zone=48 +datum=WGS84"
rst
## class       : SpatRaster 
## dimensions  : 10, 18, 1  (nrow, ncol, nlyr)
## resolution  : 111.1111, 100  (x, y)
## extent      : -1000, 1000, -100, 900  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs

Contoh rast di atas hanya membuat geometri pada data raster dengan mendefinisikan baris dan kolom, kemudian posisi pada ruang geografisnya, namun belum ada nilai pada cell tersebut.

rst <- rast(ncol=10, nrow=10)
ncell(rst)
## [1] 100
hasValues(rst)
## [1] FALSE

Gunakan fungsi values

values(rst) <- runif(ncell(rst))
hasValues(rst)
## [1] TRUE
values(rst)[1:10]
##  [1] 0.373774545 0.793159979 0.876560054 0.019921316 0.656087097 0.944766912
##  [7] 0.969056357 0.007533308 0.954169897 0.626943715
plot(rst, main='Raster dengan 100 cells')

Pada beberapa kasus jika merubah besaran kolom dan baris maka nilai yang berasosiasi pada cell tersebut akan hilang. Kebanyakan studi kasus, jika resolusi diubah akan berdampak pada jumlah baris dan kolom. Nilai cell tidak akan hilang saat merubah extent karena perubahan ini menyesuaikan resolusi, namun tidak mengubah jumlah baris atau kolom.

Hubungan antara resolusi, spatial extent, dan jumlah piksel

Raster terdiri atas sekumpulan piksel dengan dimensi dan ukuran yang sama mewakili luas area yang menutupi permukaan bumi. Ukuran luas permukaan tersebut dikenal sebagai resolusi citra spasial. Misalnya, suatu gambar memiliki resolusi spasial 1m artinya masing-masing piksel pada citra tersebut mewakili luas 1m x 1m area².

² Referensi tekait https://www.neonscience.org/resources/learning-hub/tutorials/raster-res-extent-pixels-r

Raster pada extent yang sama tetapi memiliki lebih banyak piksel akan memiliki resolusi spasial yang semakin tinggi.

Spatial extent pada raster melambangkan koordinat x dan y pada tepi ruang geografis. Informasi ini, dalam penambahan ukuran cell memberikan perintah bagaimana melakukan render setiap piksel pada ruang dimensi yang ada.

Extent dan resolusi spasial sangat berkaitan. Untuk menghitung luas suatu raster, dibutuhkan koordinat kiri bawah dari raster terebut. Dalam sistem koordinat UTM dengan satuan meter, dapat ditambahkan jumlah kolom dan baris pada tepi koordinat lokasi raster dikalikan dengan ukuran piksel (resolusi) dari raster tersebut.

Raster aljabar

Banyak fungsi generic yang memungkinkan operasi aljabar sederhana pada objek raster seperti operator +, -, *, /, operator logika >, <, ==, !, dan fungsi seperti sqrt, max, min, dan sebagainya.

r <- rast(nrow=10, ncol=10)
values(r) <- 1:ncell(r)

Lakukan operasi aljabar

s <- r + 10
s <- sqrt(s)
s <- s * r + 5
values(r) <- runif(ncell(r))
r <- round(r)
r <- r == 1

Apabila menggunakan multi-layer objek dengan nomor atau layer yang berbeda, objek yang paling sedikit akan dioperasikan secara berulang. Sebagai contoh, perkalian 4-layer objek (a1, a2, a3, a4) dengan 2-layer objek (b1, b2), hasilnya menjadi (a1*b1, a2*b2, a3*b1, a4*b2).

r <- rast(ncol=5, nrow=5)
values(r) <- 1
s <- c(r, r+1)
q <- c(r, r+2, r+4, r+6)
x <- r + s + q
x
## class       : SpatRaster 
## dimensions  : 5, 5, 4  (nrow, ncol, nlyr)
## resolution  : 72, 36  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## names       : lyr1, lyr2, lyr3, lyr4 
## min values  :    3,    6,    7,   10 
## max values  :    3,    6,    7,   10

Fungsi summary (min, max, mean, prod, sum, median, cv, range, any, all) juga mengembalikan objek SpatRaster

a <- mean(r,s,10)
b <- sum(r,s)
st <- c(r, s, a, b)
sst <- sum(st)
sst
## class       : SpatRaster 
## dimensions  : 5, 5, 1  (nrow, ncol, nlyr)
## resolution  : 72, 36  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source      : memory 
## name        :      sum 
## min value   : 17.33333 
## max value   : 17.33333

Fungsi tingkat lanjut

Aggregate dan dissaggregate

r <- rast()
values(r) <- 1:ncell(r)
ra <- aggregate(r, 20)
rd <- disagg(ra, 20)

Crop dan merge

r1 <- crop(r, ext(-50,0,0,30))
r2 <- crop(r, ext(-10,50,-20, 10))
m <- merge(r1, r2)
plot(m)

Classify

classify dapat digunakan untuk mengganti rentang nilai tertentu dengan nilai tunggal, atau mengganti nilai tunggal dengan nilai lain

r <- rast(ncol=3, nrow=2)
values(r) <- 1:ncell(r)
values(r)
##      lyr.1
## [1,]     1
## [2,]     2
## [3,]     3
## [4,]     4
## [5,]     5
## [6,]     6

plot(r)

Ubah semua nilai di atas 4 menjadi NA

s <- app(r, fun=function(x){ x[x < 4] <- NA; return(x)} )
as.matrix(s)
##      lyr.1
## [1,]    NA
## [2,]    NA
## [3,]    NA
## [4,]     4
## [5,]     5
## [6,]     6

plot(s)

Bagi raster pertama dengan 2 kali akar kuadrat raster kedua lalu ditambah 5

rs <- c(r, s)
w <- lapp(rs, fun=function(x, y){ x / (2 * sqrt(y)) + 5 } )
as.matrix(w)
##          lyr1
## [1,]       NA
## [2,]       NA
## [3,]       NA
## [4,] 6.000000
## [5,] 6.118034
## [6,] 6.224745

Hapus semua nilai pada r yang berisi NA pada w

u <- mask(r, w)
as.matrix(u)
##      lyr.1
## [1,]    NA
## [2,]    NA
## [3,]    NA
## [4,]     4
## [5,]     5
## [6,]     6

Mengidentifikasi nilai cell di u yang sama isinya dengan s

v <- u==s
as.matrix(v)
##      lyr.1
## [1,]   NaN
## [2,]   NaN
## [3,]   NaN
## [4,]     1
## [5,]     1
## [6,]     1

Klasifikasi nilai antara 0 dan 2 menjadi 1, dan seterusnya

x <- classify(w, rbind(c(0,2,1),  c(2,5,2), c(4,10,3)))
as.matrix(x)
##      lyr1
## [1,]  NaN
## [2,]  NaN
## [3,]  NaN
## [4,]    3
## [5,]    3
## [6,]    3

Mengganti 2 dengan 40 dan 3 dengan 50

y <- classify(x, cbind(id=c(2,3), v=c(40,50)))
as.matrix(y)
##      lyr1
## [1,]  NaN
## [2,]  NaN
## [3,]  NaN
## [4,]   50
## [5,]   50
## [6,]   50

Fungsi summarizing

Statisik deskriptif secara normal dapat menggunakan global untuk summary seluruh cell pada raster objek.

r <- rast(ncol=36, nrow=18)
values(r) <- runif(ncell(r))
global(r, mean)
##            mean
## lyr.1 0.4888329

zonal untuk summary berdasarkan zona atau kelas tertentu.

s <- r
values(s) <- round(runif(ncell(r)) * 5)
zonal(r, s, 'mean')
##   lyr.1     lyr.1
## 1     0 0.5108268
## 2     1 0.5249918
## 3     2 0.4981636
## 4     3 0.4841889
## 5     4 0.4483667
## 6     5 0.4533527

freq digunakan untuk menghitung area suatu kelas cell

freq(s)
##      layer value count
## [1,]     1     0    75
## [2,]     1     1   133
## [3,]     1     2   138
## [4,]     1     3   116
## [5,]     1     4   122
## [6,]     1     5    64
freq(s, value=3)
##      layer value count
## [1,]     1     3   116

# contoh lain
depok_rast <- rast("kota_depok.tif")
freq(depok_rast)
##    layer        value count
## 1      1         Beji   303
## 2      1   Bojongsari   404
## 3      1     Cilodong   320
## 4      1    Cimanggis   447
## 5      1       Cinere   216
## 6      1     Cipayung   230
## 7      1         Limo   241
## 8      1 Pancoran Mas   373
## 9      1     Sawangan   545
## 10     1   Sukma Jaya   359
## 11     1        Tapos   690

crosstab untuk membuat tabulasi silang dua objek raster

ctb <- crosstab(c(r*3, s))
head(ctb)
##      lyr.1.1
## lyr.1  0  1  2  3  4  5
##     0 12 18 29 23 25 14
##     1 22 44 40 39 43 22
##     2 31 44 39 37 39 22
##     3 10 27 30 17 15  6

Fungsi helper

Nomor cell adalah bagian yang sangat penting. Data raster dapat dianggap sebagai matriks, namun dalam SpatRaster lebih umum diperlakukan sebagai vektor. Cell diberi nomor dari cell kiri atas ke cell kanan atas lalu dilanjutkan di sisi kiri baris berikutnya, dan seterusnya hingga cell terakhir di sisi kanan bawah raster. Terdapat beberapa fungsi helper untuk menentukan nomor kolom atau baris suatu cell dan sebaliknya, serta menentukan nomor cell untuk koordinat x, y dan sebaliknya.

r <- rast(ncol=36, nrow=18)
ncol(r)
## [1] 36
nrow(r)
## [1] 18
ncell(r)
## [1] 648
rowFromCell(r, 100)
## [1] 3
colFromCell(r, 100)
## [1] 28
cellFromRowCol(r,5,5)
## [1] 149
xyFromCell(r, 100)
##       x  y
## [1,] 95 65
cellFromXY(r, cbind(0,0))
## [1] 343
colFromX(r, 0)
## [1] 19
rowFromY(r, 0)
## [1] 10

Mengakses nilai cell

Nilai cell dapat diakses dengan beberapa metode. values dapat digunakan untuk menampilkan semua nilai atau subset dalam baris atau blok tertentu

r <- rast(system.file("ex/elev.tif", package = "terra"))
v <- values(r)
v[3075:3080,]
## [1] 324 288 342 313 311 291
values(r, row=33, nrow=1, col=35, ncol=6)
##      elevation
## [1,]       324
## [2,]       288
## [3,]       342
## [4,]       313
## [5,]       311
## [6,]       291

Menampilkan nilai menggunakan nomor cell atau koordinat dengan extract

cells <- cellFromRowCol(r, 33, 35:40)
cells
## [1] 3075 3076 3077 3078 3079 3080
r[cells]
##   elevation
## 1       324
## 2       288
## 3       342
## 4       313
## 5       311
## 6       291
xy <- xyFromCell(r, cells)
xy
##             x        y
## [1,] 6.029167 49.92083
## [2,] 6.037500 49.92083
## [3,] 6.045833 49.92083
## [4,] 6.054167 49.92083
## [5,] 6.062500 49.92083
## [6,] 6.070833 49.92083
extract(r, xy)
##   elevation
## 1       324
## 2       288
## 3       342
## 4       313
## 5       311
## 6       291