5. vector-raster conversions, reprojection, warping

Edzer Pebesma

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

このビネットは、stars オブジェクトをベクタとラスタの表現から移動させる方法を示している。

sf ベクタオブジェクトのラスタライズ

library(stars)
## Loading required package: abind
system.file("gpkg/nc.gpkg", package = "sf") %>%
  read_sf() %>%
  st_transform(32119) -> nc
nc$dens = nc$BIR79 / units::set_units(st_area(nc), km^2)
(nc.st = st_rasterize(nc["dens"], dx = 5000, dy = 5000))
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##                    Min.  1st Qu. Median     Mean  3rd Qu.     Max. NA's
## dens [1/km^2] 0.2545128 1.225654 1.9322 3.345956 3.825793 21.24795 4808
## dimension(s):
##   from  to offset delta                 refsys point x/y
## x    1 162 123830  5000 NAD83 / North Carolina FALSE [x]
## y    1  61 318256 -5000 NAD83 / North Carolina FALSE [y]
plot(nc.st)

使用されるアルゴリズムは GDAL rasterize ユーティリティで、このユーティリティのすべてのオプションは st_rasterize に渡すことができる。最終的なラスタのジオメトリは、ターゲット境界ボックスとラスタ寸法 nx および ny、またはピクセルサイズ dx および dy パラメータを渡すことで制御することができる。

ラスタオブジェクトを sf オブジェクトにベクタライズする

stars オブジェクトは、 を使って オブジェクトに変換することができる。 これは、ピクセルがピクセル中心の点値を表すか、単一の値を持つ小さな正方形のポリゴンを表すかによって、いくつかのオプションがある。 st_as_sf sf

ランドサット-7の6バンドの画像で再度作業するが、最初のバンドを選択し、値を丸めることにする。

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)[, 1:50, 1:50, 1:2]
x[[1]] = round(x[[1]]/5)

ポリゴン化

ラスタセルが点描画で、フィールド全体をベクタ表示したい場合は、等高線を引いて等高線セットをエクスポートする (GDALのバージョンが2.4.0以上の場合のみ有効) 。

l =  st_contour(x, contour_lines = TRUE, breaks = 11:15)
plot(l[1], key.pos = 1, pal = sf.colors, lwd = 2, key.length = 0.8)

ポイントへのエクスポート

あるいは、単純にすべてのピクセルをポイントとしてエクスポートし、ポイントごとにすべてのバンドを持つワイドテーブルとして取得し、POINT ジオメトリの複製を作成しないこともできる。

st_as_sf(x, as_points = TRUE, merge = FALSE)
## Simple feature collection with 2500 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 288790.5 ymin: 9119350 xmax: 290187 ymax: 9120747
## Projected CRS: SIRGAS 2000 / UTM zone 25S
## First 10 features:
##    L7_ETMs.tif.V1 L7_ETMs.tif.V2                 geometry
## 1              14             11 POINT (288790.5 9120747)
## 2              14             11   POINT (288819 9120747)
## 3              13             10 POINT (288847.5 9120747)
## 4              12              9   POINT (288876 9120747)
## 5              12             10 POINT (288904.5 9120747)
## 6              12             10   POINT (288933 9120747)
## 7              12             10 POINT (288961.5 9120747)
## 8              12             10   POINT (288990 9120747)
## 9              13             10 POINT (289018.5 9120747)
## 10             13             10   POINT (289047 9120747)

または、1つの属性と全ポイントが複製された長いテーブルとして。

st_as_sf(x, as_points = TRUE, merge = FALSE, long = TRUE)
## Simple feature collection with 5000 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 288790.5 ymin: 9119350 xmax: 290187 ymax: 9120747
## Projected CRS: SIRGAS 2000 / UTM zone 25S
## First 10 features:
##    band L7_ETMs.tif                 geometry
## 1     1          14 POINT (288790.5 9120747)
## 2     1          14   POINT (288819 9120747)
## 3     1          13 POINT (288847.5 9120747)
## 4     1          12   POINT (288876 9120747)
## 5     1          12 POINT (288904.5 9120747)
## 6     1          12   POINT (288933 9120747)
## 7     1          12 POINT (288961.5 9120747)
## 8     1          12   POINT (288990 9120747)
## 9     1          13 POINT (289018.5 9120747)
## 10    1          13   POINT (289047 9120747)

このように、band という属性が追加され、どのバンドに関係しているかがわかるようになった。

ポリゴンへの書き出し

あるいは、ポリゴンにエクスポートして、次のようにピクセルごとに1つのポリゴンを取得することもできる。

st_as_sf(x[1], as_points = FALSE, merge = FALSE)
## Simple feature collection with 2500 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 288776.3 ymin: 9119336 xmax: 290201.3 ymax: 9120761
## Projected CRS: SIRGAS 2000 / UTM zone 25S
## First 10 features:
##    L7_ETMs.tif.V1 L7_ETMs.tif.V2                       geometry
## 1              14             11 POLYGON ((288776.3 9120761,...
## 2              14             11 POLYGON ((288804.8 9120761,...
## 3              13             10 POLYGON ((288833.3 9120761,...
## 4              12              9 POLYGON ((288861.8 9120761,...
## 5              12             10 POLYGON ((288890.3 9120761,...
## 6              12             10 POLYGON ((288918.8 9120761,...
## 7              12             10 POLYGON ((288947.3 9120761,...
## 8              12             10 POLYGON ((288975.8 9120761,...
## 9              13             10 POLYGON ((289004.3 9120761,...
## 10             13             10 POLYGON ((289032.8 9120761,...

または同じピクセル値を持つポリゴンをマージする。

p = st_as_sf(x, as_points = FALSE, merge = TRUE)

境界線付きでプロットすると、同じ画素値を持つ領域の境界線が分解されて表示される。

plot(p)

さらにオプション connect8 を TRUE に設定すると、デフォルトの4連結ではなく、8連結のアルゴリズムが使用される。 どちらの場合も、返されるポリゴンは、シンプルフィーチャの基準では無効であることが多いが、lwgeom::st_make_valid を使って有効にすることができる。

stars オブジェクトでベクタとラスタの切り替えを行う

ラスタ次元をベクタ次元に変換し、その他の次元は stars オブジェクトのままにしておくには、次のような方法がある。

x.sf = st_xy2sfc(x, as_points = TRUE)
x.sf
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##              Min. 1st Qu. Median    Mean 3rd Qu. Max.
## L7_ETMs.tif     7       9     11 11.2548      12   28
## dimension(s):
##          from   to                     refsys point
## geometry    1 2500 SIRGAS 2000 / UTM zone 25S  TRUE
## band        1    2                         NA    NA
##                                                       values
## geometry POINT (288790.5 9120747),...,POINT (290187 9119350)
## band                                                    NULL

のように as_points の引数を設定することも必要である。st_as_sf。

ラスタを再投影する

曲線ラスタもラスタであり、正方形および長方形のグリッドは曲線グリッドの特殊なケースであると考えると、ラスタの再投影はもはや「問題」ではなく、すべてのラスタセルに対して新しい座標を再計算するだけで、一般に曲線グリッド (時には正方形または長方形グリッドに戻すことができる) になる。曲線的なグリッドセルがセル中心の座標で表現される場合、グリッドセルの実際の形状は失われ、グリッドセルが大きい場合や変換がより強い非線形である場合に、より大きな影響を与える可能性がある。

上記で作成したグリッドを再投影した例を以下に示す。

nc.st %>% st_transform("+proj=laea +lat_0=34 +lon_0=-60") -> nc.curv
nc.curv
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##                    Min.  1st Qu. Median     Mean  3rd Qu.     Max. NA's
## dens [1/km^2] 0.2545128 1.225654 1.9322 3.345956 3.825793 21.24795 4808
## dimension(s):
##   from  to                       refsys point                         values
## x    1 162 +proj=laea +lat_0=34 +lon... FALSE [162x61] -2210936,...,-1371611
## y    1  61 +proj=laea +lat_0=34 +lon... FALSE    [162x61] 90645.7,...,538200
##   x/y
## x [x]
## y [y]
## curvilinear grid
plot(nc.curv, border = NA, graticule = TRUE)

同じラスタセルが新しいCRSに再プロットされたが、今度は曲線的なグリッドに変更された。

ラスタを歪ませる

ラスタのワーピングとは、別のCRSにある (通常は規則的な) グリッドをもとに、新しいCRSに新しい_規則的なグリッドを作成することを意味する。まずターゲットグリッドを作成することで、前節の変換を行うことができる。

nc %>% st_transform("+proj=laea +lat_0=34 +lon_0=-60") %>% st_bbox() %>%
    st_as_stars() -> newgrid

を作成し、古いラスタを新しいラスタにワープさせる。

nc.st %>% st_warp(newgrid) -> nc.new
nc.new 
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##                    Min.  1st Qu. Median     Mean  3rd Qu.     Max.  NA's
## dens [1/km^2] 0.2545128 1.225654 1.9322 3.344844 3.825793 21.24795 36155
## dimension(s):
##   from  to   offset    delta                       refsys x/y
## x    1 380 -2188108  2098.08 +proj=laea +lat_0=34 +lon... [x]
## y    1 171   494920 -2098.08 +proj=laea +lat_0=34 +lon... [y]
plot(nc.new)

この新しいオブジェクトは、新しいCRSで、新しいX軸とY軸に整列した規則的なグリッドを持つ。