fgdrパッケージで基盤地図情報（数値標高モデル）を一気にgeotiffに変換する

データのダウンロード

基盤地図情報のダウンロードは1次メッシュ単位で選択することができます。広域を一気に取得したい場合におすすめです。

ダウンロードしたデータはPackDLMap.zipという形式になっているのでこれを解凍します。解凍してできたディレクトリの中にもさらにzipファイルがたくさん格納されています。最大で70個ぐらいのファイルが格納されています。

一気に解凍する

対象が単一のディレクトリだけならそれほど大変ではないのですが、北海道全域ともなると40以上のディレクトリが処理対象となります。再帰的にunzip処理を施します。

あらかじめ、メタデータの入ったxmlファイルを削除し、その後一気にunzip処理をします。

# 不要なxmlを削除する
find . -name "*.xml" -delete

# 配下の複数のディレクトリに再帰的にunzipを適用する場合
# その場に出力する場合は -execdir を使う　コマンドの直下のディレクトリの場合は -execを使う
find . -name "*.zip" -execdir unzip -o {} \;

fgdrパッケージを使って一気にgeotiffに変換する

ここからRでの処理になります。
二重for文を回しているだけのシンプルなコードですが、30GBぐらいのメモリがないと落ちます。

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

# ディレクトリの一覧を取得
# ただし、配下のディレクトリの中のファイルも取得する
dir <- list.dirs(path = "親ディレクトリのパス", full.names = TRUE, recursive = TRUE)[-1]

## DEMのしゅつりょｋ
#1 ダウンロードして解凍したxmlファイルのリストを作成
for(i in 1:length(dir)) {
  xml <- 
    list.files(path = dir[i], 
               pattern = ".*\\.xml", 
               full.names = T) 
  #2 最初のxmlを読み込み
  dem <- read_fgd_dem(xml[1], resolution = 10 , return_class = "raster")  # 最初のラスタ
  #3 for文でxmlファイルをラスタオブジェクトに変換して結合
  for(j in 2:length(xml) ) {
    dem1 <- read_fgd_dem(xml[j], 
                         resolution = 10 ,
                         return_class = "raster" )
    dem <- merge(dem , dem1, tolerance = 0.5)
  }
  #4  海域=-9999に0を代入
  dem[dem<0] <- 0
  #5 座標系を変換
  dem_utm54 <- projectRaster(dem , crs = CRS("+init=epsg:3100") )
  #6 tiffに出力 -> 元のディレクトリに出力
  writeRaster(dem_utm54,
              paste0(dir[i],"/","DEM_JGD2000utm54.tif"),
              overwrite = T)
}

参考サイト

開発者のUryu Shinyaさんのサイトを参考にしました。