国土地理院基盤地図情報をRを使ってGISデータに変換する

国土地理院の基盤地図情報は独自のxmlファイルで記述されています。
基本項目は国土地理院から「基盤地図情報ビューア」というコンバーターが提供されており、数値標高モデルはエコリスさんから「基盤地図情報標高DEMデータ変換ツール」というものが提供されています。普通はこの2種のコンバーターを使用して、基盤知事情報のデータをQGISなどで使用できる一般的なGISデータ（tiffやshp形式）に変換します。

しかし、これらのコンバーターはいずれもWindows専用で、MacOSやLinuxでは使えません（Wineを使用する方法もありますがうまくいかないことも多いです）。そこで、Rのfgdrパッケージを使用して基盤地図情報のデータの変換を行います。

パッケージ読み込み

最低限必要なパッケージは次のとおりです。

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

ダウンロードしたデータの処理

ダウンロードした基盤地図情報のzipファイルを解凍するとさらにzipファイルが現れます。

一つ一つ解凍するのは大変なのでunzipコマンドで一括解凍します。unzipコマンドのワイルドカードはちょっと変わっていて、クォーテーションマークで囲みます。

unzip '*zip'

解凍が終わるとxmlファイルが現れます。

数値標高モデルの変換

まずは、標高データをGeoTiffに変換します。
すでに国土地理院の数値標高モデルデータをラスタとしてRで扱うで紹介されている手法です。
あらかじめデータを格納するmapディレクトリを作成し、解凍したxmlファイルの入ったディレクトリをその直下に置いています。

# 1 ダウンロードして解答したxmlファイルのリストを作成します。
xml <- 
  list.files(path="map/PackDLMap_DEM/", 
             pattern=".*\\.xml", 
             full.names=T) 
# 2 最初のラスタを読み込む
dem <- 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(dem , dem1,tolerance = 0.5)
}
# 4  海域 = -9999に0を代入
dem[dem<0] <- 0
# 地図を表示
plot(dem)
# 5 座標系を変換
dem_12k <- projectRaster(dem , crs = CRS("+init=epsg:2454") )
# 6 tiffに出力
writeRaster(dem_12k, "map/FranoDEM_JGD2000_12k.tif")

基本項目

続いて基本項目（道路や河川）の変換です。こちらもベースは「Rで国土地理院 基盤地図情報データを扱う」を参考にしています。

こちらはpurrr::map関数によるループ処理を利用しています。

# 水域
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "WA",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/WA_JGD2000_12k.gpkg")  # ファイルを書き出し

# 河川
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "WL",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/WL_JGD2000_12k.gpkg")  # ファイルを書き出し


# 等高線
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "Contr",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/Contr_JGD2000_12k.gpkg") 

# 建物
list.files(path = "map/PackDLMap_kihon/",  
                  pattern = "BldA",
                  full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/BldA_JGD2000_12k.gpkg") 

# 鉄道
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "RailCL",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/RailCL_JGD2000_12k.gpkg") 

# 市町村域_polygon
list.files(path = "map/PackDLMap_kihon/",  
                  pattern = "AdmArea",
                  full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  group_by(name) %>%  # 市町村名でグループ化
  summarize() %>%  # 市町村名で結合
  st_write(driver = "GPKG" ,
           append = FALSE , "map/AdmArea_JGD2000_12k.gpkg") 

# 市町村域_Line
list.files(path = "map/PackDLMap_kihon/",  
           pattern = "AdmBdry",
           full.names = T) %>%  # 基盤地図情報のxmlファイルをリスト化
  map(read_fgd) %>%  # xmlファイルにread_fgdを適用
  bind_rows() %>%  # xmlファイルを結合
  st_transform(2454) %>%  # JGD2000平面直角座標系12系に変換
  st_write(driver = "GPKG" ,
           append = FALSE , "map/AdmBdry_JGD2000_12k.gpkg") 

ポリゴンを属性で結合する方法でつまずいたのですが、なんのことはなく、summarize()で良いということがわかりました。こちらを参考にしました。（Rでポリゴンを融合する（実践編））

以上の作業で出力したtiffやgpkgファイルを読み込んで、QGISで表示できました。