[R] starsとsfを用いた，ラスタデータとベクタデータを組み合わせたデータ抽出方法の調査

目標

この記事の目標は以下の通りです．

• {stars}のラスタデータを，{sf}のジオメトリを用いて抽出する
• ↑で抽出したラスタデータを，各ジオメトリごとに要約する
• POINTジオメトリがある位置のラスタデータ値を取得する

参考にしたページ

• https://r-spatial.github.io/stars/reference/index.html
• https://r-spatial.github.io/sf/reference/index.html

使用するデータ

ラスタデータは国土地理院が公開している日本の標高データ を使用します．
ベクタデータはesriジャパンが公開している全国市区町村界データを使用します．

(左)標高を表すラスタデータ (右)全国の市区町村の境界を表すベクタデータ

また，データに対してあらかじめ以下の処理をおこなっておきます．

#データ読み込み
jp_sf <- st_read("ここにベクタデータのパス")
jp_elev <- read_stars("ここにラスタデータのパス")

#座標系を設定
crs <- 4326
jp_sf <- st_transform(jp_sf, 4326)
jp_elev <- st_set_crs(jp_elev, 4326)

#ラスタデータの値を数値型に変更
jp_elev$el.tif <- as.numeric(jp_elev$el.tif)

#st_cropなど重ね合わせの処理をする場合，以下のコードを実行しないと動かない
sf_use_s2(FALSE)

ラスタデータの値をジオメトリを用いて抽出する

まず以下のコードで取得した，静岡県を表すポリゴンジオメトリを用いて，ポリゴン内にあるラスタデータを抽出してみます．

shizuoka_sf <- jp_sf %>% filter(KEN == "静岡県")
ggplot() + geom_sf(data=shizuoka_sf)

{sf}によって表現されるジオメトリを用いた{stars}のラスタデータの抽出は，
st_crop()によって実現できます．

shizuoka_elev <- st_crop(jp_elev, shizuoka_sf)

ggplot() + geom_stars(data=shizuoka_elev) + geom_sf(data=shizuoka_sf, alpha=0.2) +
  scale_fill_viridis_c() + theme(legend.position="none")

抽出したラスタデータを各ジオメトリごとに要約する

先の処理で取得した静岡県の標高データに関して，各自地区ごとの平均値を求めてみたいと思います．
これは，st_extract()を用いて実現できます．

shizuoka_elev_means <- st_extract(shizuoka_elev, shizuoka_sf) %>% st_as_sf()

ggplot() + geom_sf(data=shizuoka_elev_means, mapping = aes(fill=el.tif))

この処理はaggregate()でも実現できます．

shizuoka_elev_means <- aggregate(shizuoka_elev, shizuoka_sf, FUN = mean) %>% st_as_sf()

POINTジオメトリがある位置のラスタデータ値を取得する

最後に，以下の処理で得られる静岡県の各自地区の重心地点の標高を抽出してみます．

shizuoka_centroids <- st_centroid(shizuoka_sf)

ggplot() + geom_sf(data=shizuoka_sf) + geom_sf(data=shizuoka_centroids, color="red")

これは，st_extract()を用いて実現できます．

shizuoka_centroids_elev <- st_extract(shizuoka_elev, shizuoka_centroids)

print(shizuoka_centroids_elev, n=nrow(shizuoka_centroids_elev))

出力:

Simple feature collection with 43 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 137.5265 ymin: 34.64836 xmax: 139.0872 ymax: 35.36001
Geodetic CRS:  WGS 84
   el.tif                  geometry
1      90 POINT (138.2642 35.23961)
2       3 POINT (138.3823 34.95176)
3      29 POINT (138.4855 35.09715)
4       5 POINT (137.7196 34.72904)
5       2 POINT (137.7779 34.74979)
6       3 POINT (137.6368 34.72872)
7       2  POINT (137.7473 34.6808)
8       8 POINT (137.6404 34.83672)
9       5 POINT (137.7823 34.82072)
10     97  POINT (137.8849 35.0598)
11     NA POINT (138.8306 35.07731)
12     15 POINT (139.0632 35.09066)
13     14 POINT (138.9489 35.13941)
14     51 POINT (138.6088 35.30254)
15     15 POINT (139.0872 34.93564)
16     20 POINT (138.1266 34.91221)
17     24  POINT (138.6991 35.2011)
18      6 POINT (137.8505 34.74466)
19      2 POINT (138.3013 34.84019)
20      5 POINT (138.0189 34.78002)
21      6 POINT (138.2338 34.91958)
22     59  POINT (138.8791 35.2989)
23      2 POINT (137.9281 34.74177)
24     21 POINT (138.9213 34.70887)
25     38  POINT (138.8774 35.2255)
26      2 POINT (137.5265 34.72268)
27     15 POINT (138.9259 34.92231)
28      7 POINT (138.1483 34.64836)
29      6 POINT (138.1004 34.73884)
30     19  POINT (138.973 35.03799)
31      7 POINT (138.1857 34.73189)
32     41  POINT (139.0283 34.8266)
33     26 POINT (138.9526 34.78611)
34     11 POINT (138.8315 34.66438)
35     42 POINT (138.8204 34.74584)
36     45 POINT (138.8173 34.82334)
37     28 POINT (138.9941 35.11396)
38      3  POINT (138.8987 35.0975)
39     19  POINT (138.871 35.16614)
40     57 POINT (138.9194 35.36001)
41      2 POINT (138.2572 34.77345)
42    103 POINT (138.1033 35.16363)
43     15 POINT (137.9476 34.88657)

各重心点の標高値が取得できました．

今後の調査

　この記事では{stars}と{sf}を組み合わせたデータの抽出方法について調査しました．
今後は{stars},{sf}を，{spstat}や{spdep}など他の空間分析系のパッケージと併せて使用し，空間データを分析する方法について調査する予定です．