4. stars data model

Edzer Pebesma

star の vignette のより良いバージョンは、 https://r-spatial.github.io/stars/articles/ を参照。

この vignette は、stars オブジェクトのデータモデルを、人工的なデータセットと実際のデータセットを使って説明している。

#star 座オブジェクト {#stars-objects}

stars オブジェクトは

dimensions オブジェクトは、dimension 要素の名前付きリストで、それぞれがデータアレイの次元 (空間、時間、型など) のセマンティクスを記述している。それに加えて、dimensions オブジェクトは、クラス stars_raster の raster という属性を持ち、これは3つの要素を持つ名前付きリストである。

affine と curvilinear の値は、ラスタデータの場合にのみ関係し、dimensions で非NA値であることが示される。

dimension オブジェクトは、_single_dimension を記述するもので、名前付きの要素を持つリストである。

from と は通常1とディメンションサイズになるが、 はサブグリッドが選択された (またはトリミングされた) 場合、1より大きくなることがある。 to from

offset および は、_regularly_離散化された次元にのみ適用され、そうでない場合は 。 である場合、寸法値は フィールドに保持されるかもしれない。 直線および曲線グリッドは、 にグリッド値が必要で、次のいずれかになる。 delta NA NA values values

あるいは、values は、sfc ベクタ (「リスト列」) にエンコードされた空間ジオメトリのセットを含むことができ、この場合、 [vector data cube] (https://r-spatial.org/r/2022/09/12/vdc.html) を持つ。

グリッドタイプ

正規のグリッド \(4 \times 5\) の行列から作成された非常にシンプルなファイル付き {#regular-grids

-with-a-very-simple-file-created-from-a-\(4-\times-5\)-matrix}

suppressPackageStartupMessages(library(stars))
m = matrix(1:20, nrow = 5, ncol = 4)
dim(m) = c(x = 5, y = 4) # named dim
(s = st_as_stars(m))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##     Min. 1st Qu. Median Mean 3rd Qu. Max.
## A1     1    5.75   10.5 10.5   15.25   20
## dimension(s):
##   from to offset delta point x/y
## x    1  5      0     1 FALSE [x]
## y    1  4      0     1 FALSE [y]

見るからに

dim(s[[1]])
## x y 
## 5 4

このオブジェクトをプロットするとき、image のメソッドを使って、stars のオブジェクトをプロットする。

image(s, text_values = TRUE, axes = TRUE)

を見ると、 \((0,0)\) はグリッドの原点 (グリッドの角) 、 \(1\) はあるインデックス (行、列) から次のインデックスへ向かって増加する座標値であることがわかる。これは、連続した行列の列が、南から北へ向かうグリッド線を表していることを意味する。この方法で定義されたグリッドは 正則 であり、グリッドセルのサイズはどこでも一定である。

実際のグリッドデータセットの多くは、y座標 (グリッド行) が北から南 (上から下) に向かっている。これは、delta の負の値で実現される。 \((0,0)\) のグリッドは変化していないことがわかる。

attr(s, "dimensions")[[2]]$delta = -1
image(s, text_values = TRUE, axes = TRUE)

例えば、パッケージに含まれる GeoTIFF は、GDAL を通して読み込まれるすべてのデータソースと同様に、y -座標に負の delta を持っている。

tif = system.file("tif/L7_ETMs.tif", package = "stars")
st_dimensions(read_stars(tif))["y"]
##   from  to  offset delta                     refsys point
## y    1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE

ラスタ属性、回転グリッド、シアーグリッド

stars オブジェクトのディメンションテーブルは、raster 属性を持つ。

str(attr(st_dimensions(s), "raster"))
## List of 4
##  $ affine     : num [1:2] 0 0
##  $ dimensions : chr [1:2] "x" "y"
##  $ curvilinear: logi FALSE
##  $ blocksizes : NULL
##  - attr(*, "class")= chr "stars_raster"

を保持するリストである。

これらのフィールドは、個々の次元よりも高いレベルで配列の特性を記述するため、このレベルで必要とされる。次元のペアはラスタを形成し、affine と curvilinear は、次元ごとに行うことができない場合に、x と y がグリッドインデックス (下記参照) から as a pair をどのように派生させるかを記述している。

2つのアフィンパラメータ \(a_1\) と \(a_2\) を用いて、 \(x\) と \(y\) の座標は、 (1ベースの) グリッドインデックス \(i\) と \(j\)、グリッドオフセット値 \(o_x\) と \(o_y\)、グリッドセルサイズ \(d_x\) と \(d_y\) から、以下のように導き出される。

\[x = o_x + (i-1) d_x + (j-1) a_1\] \[y = o_y + (i-1) a_2 + (j-1) d_y\] 明らかに、 \(a_1=a_2=0\)、 \(x\)、 \(y\) は、それぞれのインデックス、オフセット、セルサイズからすべて導き出される。

なお、整数インデックスの場合、座標はグリッドセルの開始辺のものである。左上のグリッドセルの中心を得るには (負の数の場合 \(d_y\) ) 、 \(i=1.5\) と \(j=1.5\) を用いる。

\(a_1\) と \(a_2\) を 0 以外の値に設定することで、グリッドを回転させることができる。

attr(attr(s, "dimensions"), "raster")$affine = c(0.1, 0.1)
plot(st_as_sf(s, as_points = FALSE), axes = TRUE, nbreaks = 20)

回転角度を度数で表すと

atan2(0.1, 1) * 180 / pi
## [1] 5.710593

2つの回転係数 ( \(a_1\) と \(a_2\) ) が不等である場合、シェアードグリッドが得られる。

attr(attr(s, "dimensions"), "raster")$affine = c(0.1, 0.2)
plot(st_as_sf(s, as_points = FALSE), axes = TRUE, nbreaks = 20)

ここで、Y軸とX軸はそれぞれ異なる度数で回転している

atan2(c(0.1, 0.2), 1) * 180 / pi
## [1]  5.710593 11.309932

レクトリニアグリッド

Rectilinear gridsは軸が直交しているが、セルが合同 (同じ大きさ、形) ではなく、各軸がそれぞれ不規則に細分化されている。

セルの境界を指定することで、直方体のグリッドを定義することができる。つまり、各次元に対して、次元サイズより1つ多い値を指定することになる。

x = c(0, 0.5, 1, 2, 4, 5)  # 6 numbers: boundaries!
y = c(0.3, 0.5, 1, 2, 2.2) # 5 numbers: boundaries!
(r = st_as_stars(list(m = m), dimensions = st_dimensions(x = x, y = y)))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##    Min. 1st Qu. Median Mean 3rd Qu. Max.
## m     1    5.75   10.5 10.5   15.25   20
## dimension(s):
##   from to point                values x/y
## x    1  5 FALSE     [0,0.5),...,[4,5) [x]
## y    1  4 FALSE [0.3,0.5),...,[2,2.2) [y]
st_bbox(r)
## xmin ymin xmax ymax 
##  0.0  0.3  5.0  2.2
image(r, axes = TRUE, col = grey((1:20)/20))

最後の値を省略すると、stars、最後のセルの境界は、1つ前のセルの幅から導き出されるため、異なるセル境界を作成する可能性がある。

x = c(0, 0.5, 1, 2, 4)  # 5 numbers: offsets only!
y = c(0.3, 0.5, 1, 2)   # 4 numbers: offsets only!
(r = st_as_stars(list(m = m), dimensions = st_dimensions(x = x, y = y)))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##    Min. 1st Qu. Median Mean 3rd Qu. Max.
## m     1    5.75   10.5 10.5   15.25   20
## dimension(s):
##   from to point              values x/y
## x    1  5 FALSE   [0,0.5),...,[4,6) [x]
## y    1  4 FALSE [0.3,0.5),...,[2,3) [y]
st_bbox(r)
## xmin ymin xmax ymax 
##  0.0  0.3  6.0  3.0

セルの幅が一定であれば問題ないが、その場合、上限の境界が与えられているかどうかに関わらず、境界は offset と delta の値に縮小される。

x = c(0, 1, 2, 3, 4)  # 5 numbers: offsets only!
y = c(0.5, 1, 1.5, 2)   # 4 numbers: offsets only!
(r = st_as_stars(list(m = m), dimensions = st_dimensions(x = x, y = y)))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##    Min. 1st Qu. Median Mean 3rd Qu. Max.
## m     1    5.75   10.5 10.5   15.25   20
## dimension(s):
##   from to offset delta point x/y
## x    1  5      0     1 FALSE [x]
## y    1  4    0.5   0.5 FALSE [y]
st_bbox(r)
## xmin ymin xmax ymax 
##  0.0  0.5  5.0  2.5

また、st_dimensions の呼び出しに引数 cell_midpoints を指定することで、cell midpoints を設定することも可能である。

x = st_as_stars(matrix(1:9, 3, 3), 
                st_dimensions(x = c(1, 2, 3), y = c(2, 3, 10), cell_midpoints = TRUE))

寸法が規則的である場合、offset は半分の delta で後ろにシフトされ、さもなければセル中心間の距離から得られる間隔でシフトされる結果になる。 セル境界が指定されている場合、このようなことは明らかに避けなければならない。

曲線的なグリッド

曲線グリッドとは、グリッド線が直線でないグリッドのことである。曲率をパラメトリックに記述するのではなく、典型的な (HDF5 または NetCDF) ファイルでは、2つのラスタレイヤに残りのレイヤの対応するピクセルごとに経度と緯度が記述されている。

例として、Sentinel 5P のデータセットを使用する。このパッケージは、starsdata でインストールできる。

install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", type = "source")

データセットはこちらでご覧いただける。

(s5p = system.file("sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", package = "starsdata"))
## [1] ""

右のサブアレイで read_stars を呼び出すと、曲線的な stars ラスタを構築することができる。

subs = gdal_subdatasets(s5p)
subs[[6]]

この配列については、項目 GEOLOCATION の下に GDAL メタデータを見ることができる。

gdal_metadata(subs[[6]], "GEOLOCATION")

で、このデータセットのどこに経度と緯度の配列が保存されているかがわかる。

nit.c = read_stars(subs[[6]]) 
threshold = units::set_units(9e+36, mol/m^2)
nit.c[[1]][nit.c[[1]] > threshold] = NA
nit.c

曲線配列は、その次元テーブルに経度と緯度の値を読み込んだ実際の配列 (ラスタレイヤー、マトリクス) を持っている。このファイルをプロットすることができる。

plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, 
         pch = 16,  logz = TRUE, key.length = 1)
maps::map('world', add = TRUE, col = 'red')
plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, 
         border = NA, logz = TRUE, key.length = 1)
maps::map('world', add = TRUE, col = 'red')

でダウンサンプリングすることができる。

(nit.c_ds = stars:::st_downsample(nit.c, 8))
plot(nit.c_ds, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, 
         pch = 16, logz = TRUE, key.length = 1)
maps::map('world', add = TRUE, col = 'red')

というのは見栄えが悪いが、セルをポリゴンでプロットした方が見栄えが良い。

plot(nit.c_ds, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, 
         border = NA, logz = TRUE, key.length = 1)
maps::map('world', add = TRUE, col = 'red')

もう一つの方法は、曲線グリッドを規則的なグリッドにワープさせることである。

w = st_warp(nit.c, crs = 4326, cellsize = 0.25)
plot(w)