やりたいこと
- 北海道の9つの2次メッシュごとに作成したラスタデータをリストにして一括で読み込みたい
- それぞれの2次メッシュにはDEM、傾斜量、日照量など複数のラスタがある
- 2次メッシュごとのラスタはstackに格納したい
全コード
# パッケージ読み込み
library(tidyverse)
library(terra)
# ラスタファイルが含まれるディレクトリのパスを指定
dir_paths <-
list.dirs(path = "data/rast")
dir_paths <- dir_paths[-1] # リストの先頭は親ディレクトリなので削除する
# 各ディレクトリ内のラスタファイルのファイルパスをリストに読み込む関数を定義
get_raster_files <- function(dir_path) {
list.files(dir_path, pattern = ".tif$", full.names = TRUE)
}
# ディレクトリごとにラスタファイルのファイルパスをリストに分割する関数を定義
split_raster_files <- function(raster_files) {
map(raster_files, rast)
}
# 各ディレクトリ内のラスタファイルをリストに読み込む
raster_list <- dir_paths %>%
map(get_raster_files) %>%
map(split_raster_files)
# 領域（extent）が不揃いなのでresample関数でDEMの領域に統一する
raster_list_ex <- raster_list # リストを新規作成
# for文でリストごとにextentを統一する
for(i in 1:length(raster_list)){
for(j in 1:length(raster_list[[i]])){
terra::resample(
x = raster_list[[i]][[j]],
y = raster_list[[i]][[1]],
method = 'bilinear'
) -> raster_list_ex[[i]][[j]]
}
}
# ラスタをスタックにまとめる
raster_list_stack <-
raster_list_ex %>%
map(rast)
ラスタファイルが格納されているディレクトリのリストを作る
# ラスタファイルが含まれるディレクトリのパスを指定
dir_paths <-
list.dirs(path = "data/rast")
dir_pathsに格納されているのはこんな感じのディレクトリです。
それぞれのディレクトリには複数のラスタファイルが格納されています。
dir_pathsは次のようなリストになっています。
# dir_pathsの中身
dir_paths
[1] "data/rast" "data/rast/624065" "data/rast/644143" "data/rast/644263" "data/rast/644322" "data/rast/644454" "data/rast/654252" "data/rast/654367" "data/rast/664243"
[10] "data/rast/664320"
この中で最初の"data/rast"は親ディレクトリなので、削除します。
dir_paths <- dir_paths[-1] # リストの先頭は親ディレクトリなので削除する
ディレクトリ内のディレクトリに格納されたファイルパスを取得する
dir_pathsに格納されたディレクトリ内にラスタデータが格納されています。このラスタデータのファイルパスを取得する関数を定義します。
# 各ディレクトリ内のラスタファイルのファイルパスをリストに読み込む関数
get_raster_files <- function(dir_path) {
list.files(dir_path, pattern = ".tif$", full.names = TRUE)
}
この関数の働きは、list.files関数で第一引数で受け取ったオブジェクトのうちtif形式のファイルのリストを作成します。
ディレクトリごとにラスタファイルを読み込む
# ディレクトリごとにラスタファイルのファイルパスをリストに分割する関数
split_raster_files <- function(raster_files) {
map(raster_files, terra::rast)
}
この関数の働きは、第一引数で受け取ったファイルパスのリストにpurrr::map関数を適用して、ディテクトリごとのリスト内にラスタファイルをSpatRaster形式のデータで格納します。
関数を実行して、ラスタファイルをリストに格納する
# 各ディレクトリ内のラスタファイルをリストに読み込む
raster_list <- dir_paths %>%
purrr::map(get_raster_files) %>% # ラスタファイル名を取得する
purrr::map(split_raster_files) # SpatRaster形式でファイルを読み込む
dir_pathsオブジェクトに格納されたディレクトリのリストを、先ほど作成したget_raster_files関数をpurrr::map関数を通して適用することで、ディレクトリごとにリスト化してファイルパスを読み込みます。
リスト化したファイルパスにsplit_raster_files関数をpurrr::map関数を通して適用することで、リスト内にラスタファイルを読み込みます。
raster_listに格納されたのはこのようなデータです。
raster_list
[[1]]
[[1]][[1]]
class : SpatRaster
dimensions : 904, 1012, 1 (nrow, ncol, nlyr)
resolution : 10.28001, 10.28001 (x, y)
extent : 468862.5, 479265.8, 4631299, 4640592 (xmin, xmax, ymin, ymax)
coord. ref. : JGD2000 / UTM zone 54N (EPSG:3100)
source : aspect.tif
name : aspect
min value : 1
max value : 4
[[1]][[2]]
class : SpatRaster
dimensions : 904, 1012, 1 (nrow, ncol, nlyr)
resolution : 10.28001, 10.28001 (x, y)
extent : 468862.5, 479265.8, 4631299, 4640592 (xmin, xmax, ymin, ymax)
coord. ref. : JGD2000 / UTM zone 54N (EPSG:3100)
source : DEM.tif
name : DEM
min value : 3.841128
max value : 691.401672
入れ子になったリスト内にSpatRaster形式のデータが格納されています。たとえば、リストの一番最後に格納されているラスタをビジュアル化すると次のようになります。
plot(raster_list[[9]][[7]])
領域（extent）を統一する
このままでは同じ領域であるはずのDEM、傾斜量などのラスタの領域が統一されていないのでこれをterra::resample関数を用いて、DEMの領域に統一します。
# 領域（extent）が不揃いなのでresample関数DEMの領域に統一する
raster_list_ex <- raster_list # リストを新規作成
# for文でリストごとにextentを統一する
for(i in 1:length(raster_list)){
for(j in 1:length(raster_list[[i]])){
terra::resample(
x = raster_list[[i]][[j]], # X = 統一される側のラスタ
y = raster_list[[i]][[1]], # Y = 基準となるラスタ（DEM）
method = 'bilinear'
) -> raster_list_ex[[i]][[j]]
}
}
この関数を適用することで、同じ範囲のラスタの領域がきっちりと統一されました。extentの統一がなされないと、これ以降の処理がうまく回りません。
ラスタをStackにまとめる
Stackはrasterパッケージやterraパッケージで用いられるマルチレイヤのデータ形式です。同一範囲、同一解像度である必要があります。
# ラスタをスタックにまとめる
raster_list_stack <-
raster_list_ex %>%
purrr::map(rast)
RasterStackに格納されたデータは次のようになります。
raster_list_stack
[[1]]
class : SpatRaster
dimensions : 904, 1012, 7 (nrow, ncol, nlyr)
resolution : 10.28001, 10.28001 (x, y)
extent : 468862.5, 479265.8, 4631299, 4640592 (xmin, xmax, ymin, ymax)
coord. ref. : JGD2000 / UTM zone 54N (EPSG:3100)
source(s) : memory
varnames : aspect
aspect
aspect
...
names : aspect, DEM, irradiation, sky_v~actor, slope, TRI, ...
min values : 1, 3.605107, 3970.149, 0.5063446, 0.00000, 0.00000, ...
max values : 4, 711.388428, 7830.668, 0.9999604, 53.57681, 34.55105, ...
names:に記載されるaspect,DEM,irradiation...が同一範囲の各ラスタです。各Stackに格納されたデータは次のように可視化できます。
plot(raster_list_stack[[1]])
まとめ
複数のディレクトリに分けて保存されたラスタを読み込んで、リスト形式のRasterStackにかくのうしました。ラスタデータをこのような形に整形することで多数のラスタデータに同一の処理を行うことが可能となります。この後は、このラスタを使用して北海道の遺跡データに地形指標を付値する作業を行います。