国土地理院基盤地図情報数値標高モデルをマージする

複数の数値標高モデルを読み込んでマージする際の手順です。terra:rast()関数で一気にマージしたいのですが、ちょっとうまくいきませんでした。

# パッケージ読み込み
library(terra)
library(tidyverse)

resolutionエラーが出る

以下のコードでマージしてやれば良いと思ったのですが、解像度が一致しないというエラーメッセージが出てしまいます。

# GeoTIFFファイルのリスト作成
tif_list <- 
  list.files(path="rast/DEM/", 
             pattern=".*\\.tiff", 
             full.names=T) 
# purrr::map関数でSpatRaster形式でインポート
dem <-
  purrr::map(tif_list, rast)
# terra::merge()関数でラスタを結合する
dem_merge <- dem[[1]]
for(i in 2:length(dem) ){
  dem_merge <- 
    terra::merge(dem_merge , dem[[i]])
}
エラー: [merge] resolution does not match

ちなみに東西方向の結合はうまくいきます。

terra::merge(dem[[1]], dem[[2]]) %>% 
  plot()

しかし、南北方向だとうまくいかずresolutionエラーが出ます。

terra::merge(dem[[1]], dem[[3]]) %>% 
  plot()
h(simpleError(msg, call)) でエラー: 
  引数 'x' の評価中にエラーが起きました (関数 'plot' に対するメソッドの選択時): [merge] resolution does not match 

raster::merge()を使う

ちなみに、raster::merge()でも同じようなエラーが出ます。

library(raster)
dem <-
  purrr:map(tif_list, raster)
raster::merge(dem[[1]], dem[[3]]) %>% 
  plot() 
h(simpleError(msg, call)) でエラー: 
  引数 'x' の評価中にエラーが起きました (関数 'plot' に対するメソッドの選択時): different origin 

ただし、raster::merge()の場合は、tolerance引数を0.2などにしてやるとうまく結合できるようです。

dem_merge <- dem[[1]]
for(i in 2:length(dem) ){
  dem_merge <- 
    raster::merge(dem_merge , dem[[i]], tolerance = 0.2)
}

raster形式をSpatRasterに変換してあげれば、無事にSpatRaster形式になります。
これが解決法の一つ。ただし、図のようになぜかうまくいかないことがしばしばあります。この場合は、左上がおかしなことになっています。

dem_merge_spraster <- rast(dem_merge)
plot(dem_merge_spraster)

fgdr::read_fgd_dem()関数でresolutionを指定して結合する

こちらが一番確実な方法です。地理院のxmlをGeoTIFFに変換する手間も不要です。

### 地理院xml読み込み
#1 ダウンロードして解凍したxmlファイルを一つのフォルダに格納しておく。
xml <- 
  list.files(path="rast/gml/xml/", 
             pattern=".*\\.xml", 
             full.names=T) 
#2 xmlをrasterオブジェクトに変換
dem_merge <- 
  read_fgd_dem(xml[1], resolution = 10 , return_class = "raster")  # 最初のラスタ
#3 for文でxmlファイルをラスタオブジェクトに変換して結合
for(i in 2:length(xml) ) {
  dem1 <- read_fgd_dem(xml[i], 
                       resolution = 10 ,
                       return_class = "raster" )
  dem_merge <- merge(dem_merge , dem1, tolerance = 0.5)
}
#4  海域=-9999に0を代入
dem_merge[dem_merge < 0] <- 0
# SpatRasterに変換
dem_merge_spraster <-
  terra::rast(dem_merge)
#5 描画
plot(dem_merge_spraster)