目標
この記事の目標は以下の通りです．
- {stars}パッケージの概要説明
- 基本的な操作を{tidyverse}ライクにおこなう方法の説明
パッケージの詳しい仕様については調査しておりません．
stars，Rのどちらに関しても初学者なため，誤った情報が書かれている可能性が高いです．
以下は参考にしたページです．
概要
{stars}は，画像データや時空間データなどのラスタデータ(多次元配列データ)を扱うための機能を提供します．
ラスタデータのイメージとして，以下の画像を参照してください．この画像は時空間データを表現した画像です．時空間データは，緯度(latitude)と経度(longitude)による空間的な軸と時間軸(time)を併せた，3次元配列のデータとして表現できます．
(引用: https://r-spatial.github.io/stars/index.html)
ラスタデータのファイル形式
ラスタデータのファイル形式は以下のようなものがあります．
- jpg
- png
- tiff
etc....
{stars}はこれらの形式のデータを読み込んで扱うことができます．
基本的な使い方
パッケージのインストール
install.packages("stars")
library(stars)
データの読み込みと表示
次のコードでラスタデータを読み込むことができます．
read_stars("読み込むファイルのパス")
国土地理院が公開している日本の標高データ を読み込んで表示してみます．
jp_elev <- read_stars("ここにデータへのパスを入力")
plot(jp_elev)
出力結果:
{ggplot2}を用いても出力ができます．
#geom_stars を用いる
jp_elev <- read_stars("ここにデータへのパスを入力")
ggplot() + geom_stars(data=jp_elev) + scale_fill_viridis_d()
出力結果:
starsのデータ構造
{stars}パッケージに付属している衛星写真を読み込んで表示してみます．
tif = system.file("tif/L7_ETMs.tif", package = "stars")
rast = read_stars(tif)
plot(rast, axes = TRUE)
出力結果:
ここでrastに格納されているデータを見てみます
rast
出力結果:
stars object with 3 dimensions and 1 attribute
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max.
L7_ETMs.tif 1 54 69 68.91242 86 255
dimension(s):
from to offset delta refsys point x/y
x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band 1 6 NA NA NA NA
「attribute(s)」と「dimension(s)」というデータがあります．
「attribute(s)」は，各マスに与えられている属性に関する情報が書かれています．今回読み込んだデータでは，各マスは「L7_ETMs.tif」という属性名の値を持つようです．この値が画像のピクセルの濃淡を表しています．出力では，その値の四分位範囲と最大値，最小値，平均値が表示されています．
「dimension(s)」はデータの各次元に関するデータがまとめられています．
dimension(s)の各要素の意味は次のとおりです．
| 列名 | 内容 |
|---|---|
| from | 次元配列の原点のインデックス |
| to | 次元配列の最後の値のインデックス |
| offset | 空間上の原点となる値 (例:基準とする位置の緯度や経度のこと) |
| delta | 1マスあたりに値がどのくらい変わるか (例:1マス変わるごとに緯度や経度がどのくらい増えるか) |
| refsys | 座標参照系 |
| point | マスの値が間隔の端の値を表すのか，間隔の中の値を表すのかを示す． (TRUE:間隔の中の値, FALSE:間隔の端の値) |
| values | その次元に定まる値 (色を表す軸であれば，色を表す値が収まる) |
starsの各要素には以下のようにアクセスできる．
rast[1, x軸のインデックス, y軸のインデックス, z軸のインデックス]
以下はrastの部分集合を取得する例である
#band軸[2:4]のx軸[1:100]におけるy軸の全ての値を取得している
sub_rast <- rast[1, 1:100, , 2:4]
sub_rast %>% plot()
出力結果:
{tidyverse}ライクな操作方法
{stars}は{tidyverse}と併用して処理できるように作られているようです．ここからは，{tidyverse}のメソッド群を用いてラスターデータを操作する方法をまとめます．
slice()を用いた「範囲選択」による部分集合の抽出
上述した部分集合を抽出する処理は，{dplyr}のslice()によって以下のように書き直すことができます．
rast %>% slice(band, 2:4) %>% slice(x, 1:100) %>% plot()
filterを用いた「値の条件」による部分集合の抽出
rast %>% filter(x > 288776, x < (288776 + 28.5*349/2) , band > 3) %>% plot()
出力結果:
なお，(288776 + 28.5*349/2)はx軸の定義域の中点を求めています．
xのoffsetが288776であり，deltaは28.5, インデックスは349個あることが下表からわかります．したがって，上述の式でx軸の中点が求まります．
stars object with 3 dimensions and 1 attribute
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max.
L7_ETMs.tif 1 54 69 68.91242 86 255
dimension(s):
from to offset delta refsys point x/y
x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band 1 6 NA NA NA NA
mutate()による属性の作成 ・ select()による属性の選択
mutate()を用いて，各マスに新たな属性値を与えられます．
rast %>% mutate(band2 = sqrt((L7_ETMs.tif-200)^2))
結果:
stars object with 3 dimensions and 2 attributes
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max.
L7_ETMs.tif 1 54 69 68.91242 86 255
band2 0 114 131 131.10913 146 199
dimension(s):
from to offset delta refsys point x/y
x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band 1 6 NA NA NA NA
新しい属性「band2」は．元の属性値に対して，200を引いて二乗している．これにより元画像を反転した濃淡情報をもった属性ができています．
select()を用いることで，属性を選択することができます．
rast %>% mutate(band2 = sqrt(L7_ETMs.tif-100)^2) %>% select(band2) %>% plot()
出力:
st_redimension()によるラスタデータの並列化(?)
以下のコードによって新たに作った属性を分離して並列化(?)することができます．
new_rast <- rast %>% mutate(band2 = sqrt(L7_ETMs.tif-100)^2)
new_rast = st_redimension(new_rast)
names(new_rast) = "value" #この行の処理いる？
new_rast
結果:
stars object with 4 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
Min. 1st Qu. Median Mean 3rd Qu. Max.
value 47 65 76 77.3419 87 255
dimension(s):
from to offset delta refsys point values x/y
x 1 349 288776 28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [x]
y 1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE NULL [y]
band 1 6 NA NA NA NA NULL
new_dim 1 2 NA NA NA NA L7_ETMs.tif, band2
#↑新たな軸である「ned_dim」ができており，この軸の「values」に属性の名前が割り当てられている
この処理は以下の画像をイメージするといいと思います．st_redimension()を実行する前は，1つのマスに複数の属性(値)が同時に存在しています．それを，st_redimension()によって，画像のように，(x,y,band)は変えずに属性の値だけを変えたラスタデータを横に並べて作っているイメージです．
(引用:https://r-spatial.github.io/stars/index.html)
ggplotで画像を見比べてみる
最後にggplotのfacetを使って画像を見比べてみます．途中でbandの値が6のものだけを表示するようにfilter()にかけています．
ggplot() + geom_stars(data=new_rast %>% filter(band==6)) + facet_grid(~new_dim)
出力結果:
関連記事
同じGIS系のパッケージである{sf}と組み合わせて，データを抽出する方法について調査しました．
→ https://qiita.com/piry/items/f9ecef9df8f06887fff9